Library

Mesh{T}

A subdivision of a closed interval $\Omega = [a,b]$ into n cells by the n+1 breakpoints

\[a = y_1 < y_2 < \dots < y_{n+1} = b ,\]

so that cell k is $[y_k, y_{k+1}]$.

A mesh carries no boundary condition. It is the geometry alone, and the same mesh serves a clamped basis, a recombined one and a periodic one — what differs between them is how the knot vector is closed at the two ends, which is the basis's business and not the mesh's. See BoundaryCondition.

Four families are provided. The difference between the first three matters for testing rather than for use:

• UniformMesh — equally spaced. The assembled matrices of a periodic basis are circulant, which makes several quantities vanish identically that do not vanish in general.
• GradedMesh — the image of a uniform mesh under a fixed smooth map. Refining n therefore gives a genuine family of meshes, and rates of convergence measured across it are meaningful.
• RandomMesh — randomly perturbed cell widths. A different mesh for every n, so convergence rates across it mean nothing; its use is to check the properties that must hold exactly on any mesh.
• GeneralMesh — an explicit list of breakpoints, for a subdivision that none of the three families describes.

UniformMesh(n, domain)
UniformMesh{T}(n, domain)

The uniform mesh of n cells on domain, $y_i = a + (i-1) (b-a) / n$.

domain may be a ClosedInterval such as -1 .. 1, a tuple (a, b), or a single number L standing for $[0, L]$. UniformMesh(n) is $[0, 2\pi]$.

julia> breakpoints(UniformMesh(4, 0 .. 1))
5-element Vector{Float64}:
 0.0
 0.25
 0.5
 0.75
 1.0

julia> breakpoints(UniformMesh(2, -1 .. 1))
3-element Vector{Float64}:
 -1.0
  0.0
  1.0

On a uniform mesh the B-spline basis functions of a periodic basis are translates of a single cardinal spline, so the mass, stiffness and derivative matrices assembled from it are circulant. Several identities of the discrete brackets hold on a uniform mesh and nowhere else; a test that means to check a property valid on any mesh should use RandomMesh instead.

GradedMesh(n, domain; amplitude = 0.12)
GradedMesh{T}(n, domain; amplitude = 0.12)

The image of a uniform mesh under the fixed smooth map

\[s \mapsto s + \frac{a}{2\pi} \sin 2 \pi s , \qquad s_i = \frac{i-1}{n} ,\]

with amplitude $= a$, affinely rescaled onto domain. The map does not depend on n, so refining n gives a genuine family of meshes and the rates measured across it are meaningful — which is exactly what RandomMesh does not give.

The map is a bijection of $[0,1]$ with positive derivative for $|a| < 1$, and fixes both endpoints, so the breakpoints still run from $a$ to $b$. The cell widths vary by a factor of $(1+a)/(1-a)$ across the domain.

julia> m = GradedMesh(8, 0 .. 1);

julia> issorted(breakpoints(m)) && breakpoints(m)[begin] == 0 && breakpoints(m)[end] == 1
true

RandomMesh(n, domain; seed = 1, spread = 0.6)
RandomMesh{T}(n, domain; seed = 1, spread = 0.6)

A mesh of n cells whose widths are drawn as $1 + \sigma \, r_i$ with $r_i$ uniform on $[0,1)$ and spread $= \sigma$, then rescaled to sum to the width of domain. The stream is seeded, so the mesh is reproducible.

The widths vary by up to a factor of $1 + \sigma$, and a different n gives an unrelated mesh rather than a refinement of this one. Convergence rates measured across a family of RandomMeshes are therefore meaningless; use GradedMesh for those. What this mesh is for is the properties that must hold exactly on any mesh — the antisymmetry of the discrete brackets, for instance, which a uniform mesh would confirm for the wrong reason, its assemblies being circulant.

GeneralMesh(breakpoints)
GeneralMesh{T}(breakpoints)

The mesh with the given explicit breakpoints, which must be strictly increasing. The domain is first(breakpoints) .. last(breakpoints).

This is the escape hatch for a subdivision none of the three families describes. The case it exists for is a mesh that is uniform in the interior but carries one oversized cell at each end — a device for keeping particles that stray outside the resolved region inside the support of the basis, used in particle discretisations of kinetic equations:

julia> m = GeneralMesh([-20.0; range(-10, 10; length = 5); 20.0]);

julia> ncells(m), domain(m)
(6, -20.0 .. 20.0)

julia> diff(breakpoints(m))
6-element Vector{Float64}:
 10.0
  5.0
  5.0
  5.0
  5.0
 10.0

A GeneralMesh is never circulant even when its breakpoints happen to be equally spaced, in the sense that no assembly built on it takes the CirculantMass path — that is decided by the mesh type, so an equally spaced GeneralMesh is a legitimate way to force the general path in a test.

breakpoints(m::Mesh)

The n+1 breakpoints of m, in increasing order, running from $a$ to $b$ inclusive.

These are the cell boundaries: cell k is breakpoints(m)[k] .. breakpoints(m)[k+1]. They are to be distinguished from the knot vector of a basis built on m, whose entries may repeat.

breakpoints(b::AbstractBSplineBasis)

The n+1 cell boundaries of the basis, held by the basis rather than asked of the mesh on every call.

This is what keeps evaluate_all! allocation-free, findcell being on the innermost loop of a particle deposition where one allocation per particle per step is the whole cost of the routine. UniformMesh builds its vector on demand, its breakpoints being a closed form and its own findcell a division that needs none of them; the other mesh families hold theirs, and for those the basis and the mesh share one array.

As for a mesh, the result must not be mutated: it is the basis's own storage, and where the mesh stores its breakpoints too it is the mesh's as well.

ncells(m::Mesh)

The number of cells n, one fewer than the number of breakpoints.

domain(m::Mesh)

The closed interval $[a,b]$ the mesh subdivides, as a ClosedInterval.

domain(B::TensorProductBasis)

The product domain $\Omega_1 \times \dots \times \Omega_D$, as a DomainSets ProductDomain, so that x ∈ domain(B) answers for a D-vector x.

domainlength(m::Mesh)

The width $b - a$ of the domain.

For a PeriodicBSplineBasis built on m this is the period $L$, which is why the name is not width.

meshwidth(m::Mesh)

The largest cell width, the $h$ that convergence rates are measured against.

findcell(b::AbstractBSplineBasis, x)
findcell(m::Mesh, x)

The index of the cell containing x, in 1:ncells.

The right endpoint b belongs to the last cell, not to a cell of its own — the cells are half-open except for the last, which is closed. A point outside the domain is clamped to the nearest cell for a bounded mesh; on a PeriodicBSplineBasis reduce x onto the domain first.

It allocates nothing: a basis and a non-uniform mesh both hand over a stored breakpoint vector, and on a UniformMesh the cell index is a division that does not read the breakpoints at all.

BoundaryCondition

How a B-spline basis is closed at the ends of its Mesh.

Two quite different things go by this name, and separating them is what makes the general case tractable:

How the knot vector is closed. This fixes the spline space and its dimension. Periodic continues the breakpoints periodically and gives $N = n$; Free repeats each end knot $p+1$ times — the clamped or open knot vector — and gives $N = n + p$. These are different constructions rather than variants of one, and periodicity is not a per-end setting: it couples the two ends, so it applies to a whole axis or not at all.

Homogeneous linear constraints at one end. Dirichlet $u = 0$, Neumann $u' = 0$, Robin $\alpha u + \beta u' = 0$, Natural $u'' = 0$, or a Constraint of arbitrary order. These are constraints on the clamped space, imposed by recombination — see RecombinedBSplineBasis — and each one costs one degree of freedom at the end it applies to.

Specifying them

A single condition applies to both ends; a two-tuple gives the left and the right:

julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3);                       # Free, both ends

julia> nbasis(b)
11

julia> nbasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet()))       # one per end
9

julia> nbasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3, (Dirichlet(), Free())))
10

julia> nbasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3, Periodic()))
8

Lowercase symbols are accepted as sugar and normalised on construction — :periodic, :dirichlet, :neumann, :natural, :free.

Which one conserves what

For a particle or finite-element discretisation the choice is not free: a conservation law survives the discretisation only if the conserved density lies in the span of the basis. The clamped basis reproduces every polynomial of degree $\le p$ exactly, a periodic basis reproduces only the constants, and a Dirichlet-recombined basis reproduces none — not even the constants, since every one of its functions vanishes at the ends.

polynomial_reproduction reports this as a number, so that a scheme whose conservation proof needs $1, v, v^2$ in the span can assert it rather than assume it.

See also Free, Periodic, Dirichlet, Neumann, Robin, Natural, Constraint.

Free()

No constraint: the plain clamped basis, of dimension $n + p$ on the end it applies to.

The knot vector repeats the end knot $p+1$ times, so the basis is interpolatory there — $\varphi_1(a) = 1$ and every other function vanishes at $a$. The full spline space $\mathcal{S}^p$ of $\mathcal{C}^{p-1}$ piecewise polynomials is represented, and every polynomial of degree $\le p$ is reproduced exactly.

Periodic()

The periodic closure: the breakpoints are continued periodically, $y_{i+n} = y_i + L$, and the basis is wrapped onto the torus.

Applies to a whole axis rather than to one end, since it identifies the two. Gives $N = n$, and there are no boundary functions at all: every basis function spans $p+1$ cells and the basis is $\mathcal{C}^{p-1}$ across the seam as well as inside, so integration by parts leaves no boundary terms. See PeriodicBSplineBasis.

Dirichlet()

The homogeneous Dirichlet condition $u = 0$.

On a clamped basis $\varphi_1$ is the only function that does not vanish at the end, so the recombination reduces to dropping it — the textbook elimination — and the dimension falls by one per end it applies to.

Neumann()

The homogeneous Neumann condition $u' = 0$.

Unlike Dirichlet this is a genuine recombination: both $\varphi_1$ and $\varphi_2$ have nonzero derivative at the end, and the surviving function is the combination of the two that the derivative annihilates.

Natural()

The natural boundary condition $u'' = 0$.

Requires p ≥ 2; for a lower degree the second derivative of every basis function vanishes identically and the condition constrains nothing, which is reported as an error rather than silently accepted.

Robin(α, β)

The homogeneous Robin condition $\alpha u + \beta u' = 0$.

Reduces to Dirichlet at $\beta = 0$ and to Neumann at $\alpha = 0$, though those types are preferable where they apply — they say what is meant, and Dirichlet takes the cheaper elimination path. At least one of the two coefficients must be nonzero.

The sign convention is the one written above, with both terms on the same side and the outward direction playing no role: Robin(1, 1) is $u + u' = 0$ at both ends if given as a single condition, not $u \pm u' = 0$. Pass a two-tuple where the ends differ.

Constraint(c...)

The general homogeneous local condition $\sum_{k} c_{k+1} \, D^k u = 0$, with c[1] the coefficient of $u$, c[2] of $u'$, and so on.

Constraint(1) is Dirichlet, Constraint(0, 1) is Neumann, Constraint(α, β) is Robin and Constraint(0, 0, 1) is Natural; the named types are preferable where they apply. What this adds is the conditions with no standard name — $u''' = 0$, or $u - 2u'' = 0$.

The highest derivative appearing must be at most p, since the $p$-th derivative of a degree-$p$ spline is piecewise constant and everything above it vanishes identically.

boundary_conditions(bc)

Normalise a user-supplied boundary specification to either Periodic or a pair (left, right) of BoundaryConditions.

Accepts a single condition, meaning both ends; a two-tuple, giving left and right; and a symbol or pair of symbols in place of either.

julia> boundary_conditions(Dirichlet())
(Dirichlet(), Dirichlet())

julia> boundary_conditions((:dirichlet, :neumann))
(Dirichlet(), Neumann())

julia> boundary_conditions(:periodic)
Periodic()

Periodic is returned bare rather than as a pair, because it is a condition on the axis and not on its ends: there is no such thing as being periodic at the left end only, and pairing it with anything else is rejected.

constraint_coefficients(bc::BoundaryCondition)

The coefficients $(c_0, c_1, \dots)$ of the local condition $\sum_k c_k \, D^k u = 0$ that bc imposes, as a tuple, or nothing for the conditions that impose none — Free and Periodic.

This is the one place the named conditions are given their meaning, so that the recombination has a single implementation rather than one method per condition.

constraint_order(bc::BoundaryCondition)

The highest derivative order appearing in the condition, or -1 if it imposes none.

Used to check the condition against the degree of the basis: a condition involving $D^k$ with k > p constrains nothing, because the k-th derivative of every degree-p spline vanishes identically.

nconstraints(bc::BoundaryCondition)

The number of degrees of freedom the condition removes at the end it applies to: 0 for Free and Periodic, 1 for every local condition.

AbstractBSplineBasis{T} <: Basis{T}

A B-spline basis of some degree on a Mesh, closed at the ends according to a BoundaryCondition.

Three concrete types implement it, and which one a construction produces is decided by the boundary condition rather than chosen by name — see BSplineBasis:

• BSplineBasis — the clamped basis, $N = n + p$, no constraint at either end.
• PeriodicBSplineBasis — the basis on the torus, $N = n$.
• RecombinedBSplineBasis — a clamped basis with one homogeneous constraint per constrained end, $N = n + p - \#\text{constraints}$.

All three answer nbasis, degree, order, nodes, evaluate, evaluate_all, getindex and axes, so the assembly in SplineQuadrature is written once against this interface.

BSplineBasis(mesh, p)
BSplineBasis(mesh, p, bc)
BSplineBasis{T}(mesh, p)

The B-spline basis of degree p on the Mesh mesh, closed at the ends according to the BoundaryCondition bc.

With no bc, or with Free, this is the clamped (or open) basis: the end knots are repeated $p+1$ times, the dimension is $N = n + p$, and the basis is interpolatory at both ends, $\phi_1(a) = \phi_N(b) = 1$.

julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3);

julia> nbasis(b), degree(b), order(b)
(11, 3, 4)

julia> sum(b[0.3, j] for j in eachindex(b)) ≈ 1          # partition of unity
true

julia> b[0.0, 1], b[1.0, nbasis(b)]                      # interpolatory at both ends
(1.0, 1.0)

julia> typeof(BSplineBasis(UniformMesh(8, 0 .. 1), 3, Periodic())).name.name
:PeriodicBSplineBasis

julia> typeof(BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet())).name.name
:RecombinedBSplineBasis

Knot vector

$[\,\underbrace{a, \dots, a}_{p+1}, y_2, \dots, y_n, \underbrace{b, \dots, b}_{p+1}\,]$, of length $n + 2p + 1$, so that basis function j is supported on knotvector(b)[j] .. knotvector(b)[j+p+1] and there are $N = n + p$ of them.

Repeating the end knot $p+1$ times drops the continuity there to $\mathcal{C}^{-1}$, which is what lets the basis take a nonzero value at the endpoint at all; inside, it is $\mathcal{C}^{p-1}$.

Requirements

p ≥ 0, and n ≥ 1. Unlike the periodic case there is no lower bound on the number of cells in terms of the degree: a single cell carries the full polynomial space of degree p.

See also polynomial_reproduction, which is p here and is what conservation depends on, and evaluate_all for the local evaluation a particle loop needs.

PeriodicBSplineBasis(mesh, p)
PeriodicBSplineBasis{T}(mesh, p)
PeriodicBSplineBasis(n, p; L = 2π)

The periodic B-spline basis of degree p on the Mesh mesh, spanning the spline space $\mathcal{S}^p_n$ of $\mathcal{C}^{p-1}$ piecewise polynomials of degree p on the torus obtained by identifying the two ends of $[a,b]$.

julia> b = PeriodicBSplineBasis(UniformMesh(8, 0 .. 1), 3);

julia> nbasis(b), degree(b), order(b)
(8, 3, 4)

julia> sum(b[0.3, j] for j in eachindex(b)) ≈ 1      # partition of unity
true

Dimension

The number of degrees of freedom is $N = n$, the number of cells — not $n + p$, as it is on a bounded interval. The construction differs from the clamped case only in how the knot vector is closed up: instead of repeating the end knots to clamp the basis at the two ends, the breakpoints are continued periodically, $y_{i+n} = y_i + L$, and the recursion is applied to the resulting bi-infinite sequence. The splines it produces satisfy $\phi_{j+n}^p (x) = \phi_j^p (x - L)$, so only n of them are distinct on $\Omega$, and the periodic basis is obtained by wrapping these onto the torus,

\[\phi_j^p \big\vert_\Omega (x) = \sum_{r \in \mathbb{Z}} \phi_j^p (x + rL) , \qquad 1 \le j \le n ,\]

a finite sum, since $\phi_j^p$ is supported on the p+1 cells $[y_j, y_{j+p+1}]$. The p extra functions of the clamped case are precisely those the clamping introduces at the two ends, and the wrapping identifies them in pairs.

Why this matters

There are no boundary functions at all: every basis function spans p+1 cells, and the basis is $\mathcal{C}^{p-1}$ across the seam $a \equiv b$ as well as inside. Integration by parts over $\Omega$ therefore leaves no boundary terms, which is what makes the discrete brackets assembled from this basis exactly antisymmetric.

On a UniformMesh the basis functions are in addition translates of a single cardinal B-spline, and the mass, stiffness and derivative matrices are circulant.

What is lost is polynomial reproduction beyond the constants: $v$ and $v^2$ are not periodic, so a scheme that conserves momentum or energy because they lie in the span of the basis does not do so here. See polynomial_reproduction.

Requirements

p ≥ 0 and n > p. The second is what makes the wrapping well defined: a basis function spans p+1 cells, so with n ≤ p it would wrap onto itself and the sum above would not terminate.

See also SplineQuadrature for the assembly built on this basis, and evaluate for derivatives of arbitrary order.

RecombinedBSplineBasis(parent, left, right)

A clamped BSplineBasis with a homogeneous boundary condition imposed at one or both ends by recombination: the basis functions that violate the condition are replaced by the combinations of them that satisfy it.

Usually reached through BSplineBasis rather than constructed by name:

julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet());

julia> b isa RecombinedBSplineBasis, nbasis(b)
(true, 9)

julia> û = randn(nbasis(b));

julia> abs(evaluate(b, û, 0.0)) < 1e-14, abs(evaluate(b, û, 1.0)) < 1e-14
(true, true)

julia> c = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Neumann());

julia> v̂ = randn(nbasis(c));

julia> abs(evaluate(c, v̂, 0.0, 1)) < 1e-12          # u'(a) = 0 for any coefficients
true

How the recombination works

Write the condition at the left end as $L u (a) = 0$ with $L = \sum_k c_k D^k$ of order $m$, and put $a_i = (L \phi_i)(a)$. A clamped knot vector has $D^k \phi_i (a) = 0$ for $i > k+1$, so only the first $m+1$ functions violate the condition, and only $\phi_{m+1}$ contributes to $a_{m+1} = c_m D^m \phi_{m+1}(a)$, which is nonzero whenever the leading coefficient $c_m$ is. The $m$ combinations

\[\psi_t = \phi_t - \frac{a_t}{a_{m+1}} \, \phi_{m+1} , \qquad t = 1, \dots, m ,\]

are then annihilated by $L$ by construction and span the constrained part of the space, and $\phi_{m+2}, \dots$ pass through unchanged. Counting: $m+1$ functions in, $m$ out, so each constrained end costs exactly one degree of freedom whatever the order of its condition. The right end is the mirror image, anchored on $\phi_{N_p - m}$.

For Dirichlet, $m = 0$: there are no combinations and $\phi_1$ is simply dropped, which is the textbook elimination. For Neumann, $m = 1$: one genuine combination of $\phi_1$ and $\phi_2$. For Natural, $m = 2$: $\phi_1$ and $\phi_2$ already have vanishing second derivative at $a$ and pass through, and $\phi_3$ is the one that goes.

What is lost

The recombined basis is not a partition of unity and its functions are not non-negative — $\psi_t$ is a difference. The mass matrix stays symmetric positive definite and banded, so every assembly and every solve is unaffected, but a quantity whose conservation rested on $\sum_j \phi_j \equiv 1$ no longer has it. See polynomial_reproduction, which is -1 for a Dirichlet-recombined basis: not even the constants survive.

Requirements

The two end blocks must not overlap: nbasis(parent) > m_left + m_right + 1. Enough functions must survive to span something: nbasis(parent) - nconstraints(left) - nconstraints(right) ≥ 1, which is a separate requirement because a constrained end costs one degree of freedom whatever the order of its condition. The order of each condition must be at most the degree. Periodic is a condition on the axis rather than on an end and is rejected here, as it is by boundary_conditions.

recombination_matrix(b::RecombinedBSplineBasis)

The sparse matrix R with ψ_j = Σ_i R[i,j] φ_i, of size nbasis(parent) × nbasis(b).

Every assembly of a recombined basis is the parent's assembly conjugated by this matrix — M̃ = R' * M * R for the mass matrix, Φ̃ = R' * Φ for a tabulation — which is why the quadrature needs no separate implementation for the recombined case.

recombination_matrix(B::PolarSplineBasis)

The sparse matrix R with Ψ_k = Σ_I R[I,k] Φ_I, of size nbasis(parent(B)) × nbasis(B).

Its columns 1:3 are the pole triangle and carry 2 * Nθ nonzeros each; every other column is a single one. As for a RecombinedBSplineBasis, every assembly of the polar basis is the parent's conjugated by this matrix — M̃ = R' * M * R, Φ̃ = R' * Φ — which is why PolarSplineQuadrature needs no quadrature rule of its own.

nbasis(b::AbstractBSplineBasis)

The dimension $N$ of the spline space — the number of basis functions, and the length a coefficient vector must have.

It is decided by the boundary condition, not by the degree alone:

The generic belongs to CompactBasisFunctions and is extended here rather than redefined.

degree(b::AbstractBSplineBasis)

The polynomial degree $p$ of the basis: every basis function is piecewise of degree $p$ and the space is $\mathcal{C}^{p-1}$ across an interior breakpoint.

The generic belongs to GeometricBase and is extended here rather than redefined. See order for $k = p+1$, the other convention, and why both exist.

order(b::AbstractBSplineBasis)

The order $k = p + 1$ of the spline basis b.

Note that this is the spline meaning of the word — a B-spline of order k is piecewise of degree k-1, and the knot vector of a basis of N functions has N + k entries. 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.

mesh(b::AbstractBSplineBasis)

The Mesh the basis is built on.

knotvector(b::AbstractBSplineBasis)

The knot sequence the Cox-de Boor recursion runs on.

nodes(b::BSplineBasis)

The Greville abscissae of b,

\[\xi_j = \frac{1}{p} \sum_{i=1}^{p} x_{j+i} ,\]

the averages of the p interior knots of each basis function.

These are the points a spline basis is naturally interpolated at: $\xi_j$ lies in the support of $\phi_j$ and the collocation matrix $\phi_j(\xi_i)$ is invertible (Schoenberg-Whitney). For a clamped basis $\xi_1 = a$ and $\xi_N = b$.

nodes(b::PeriodicBSplineBasis)

The Greville abscissae of b, reduced onto $[a,b)$,

\[\xi_j = \frac{1}{p} \sum_{i=1}^{p} x_{j+i} ,\]

the averages of the p interior knots of each basis function.

nodes(B::TensorProductBasis)

The per-axis node vectors, as a tuple. The tensor-product grid is their Iterators.product; it is not materialised, since at three dimensions it is the whole coefficient array.

nnodes(b::AbstractBSplineBasis)

The number of nodes, which for a spline basis is nbasis(b): there is one Greville abscissa per basis function.

The generic belongs to GeometricBase and is extended here rather than redefined.

nnodes(B::TensorProductBasis)

The number of points of the tensor-product node grid, $\prod_d N_d$ — a scalar, unlike nodes(B), which returns the per-axis vectors whose product that grid is.

grid(b::AbstractBSplineBasis)

The nodes of b, under the name ContinuumArrays uses for the points a quasi-array is sampled at. The generic belongs to that package and is extended here rather than redefined.

boundary(b::AbstractBSplineBasis)

The boundary condition of b: Periodic for a periodic basis, otherwise the pair (left, right).

basis_index(b::AbstractBSplineBasis, j::Integer)

The index of basis function j, wrapped onto 1:nbasis(b) where the basis is periodic and returned unchanged otherwise.

evaluate_all reports the first index of its local block before wrapping, because the block is contiguous only before it. Pass each index through this function rather than writing mod1 at the call site: whether the wrap is needed is a property of the basis, and getting it wrong on a non-periodic basis silently folds the two ends of the domain together.

local_width(b::AbstractBSplineBasis)

The number of entries evaluate_all! writes, i.e. the largest number of basis functions that are nonzero at a single point.

p+1 for a BSplineBasis or a PeriodicBSplineBasis. For a RecombinedBSplineBasis it can be larger, because a recombined function near an end spans the union of two parent supports; the block is still contiguous, and entries beyond the nonzeros of a particular cell are filled with zeros.

local_indices(b::AbstractBSplineBasis, cell)

The indices of the basis functions that are nonzero on cell, as an iterable of Int.

cell:(cell+p) for a clamped BSplineBasis; the same block wrapped onto 1:N for a PeriodicBSplineBasis; and the precomputed column range for a RecombinedBSplineBasis, which may be wider than p+1 near an end.

The indices are distinct, which is what lets an assembly push one entry per pair without having to worry about a duplicate being summed into place.

polynomial_reproduction(b::AbstractBSplineBasis)

The largest $m$ such that every polynomial of degree $\le m$ lies in the span of b, or -1 if not even the constants do.

This is the predicate that decides which moments an $L^2$ projection preserves. Projecting a distribution onto b reproduces $\int \pi(v) f \, dv$ for a polynomial $\pi$ exactly when $\pi$ lies in the span, so a scheme whose mass, momentum and energy diagnostics are computed from the projected distribution — $1$, $v$ and $v^2$ — needs polynomial_reproduction(b) ≥ 2.

julia> m = UniformMesh(16, -10 .. 10);

julia> polynomial_reproduction(BSplineBasis(m, 3))
3

julia> polynomial_reproduction(BSplineBasis(m, 3, Periodic()))
0

julia> polynomial_reproduction(BSplineBasis(m, 3, Dirichlet()))
-1

julia> polynomial_reproduction(BSplineBasis(m, 3, Natural()))
1

evaluate(b::AbstractBSplineBasis, j, x, d = 0)
evaluate(b::AbstractBSplineBasis, û::AbstractVector, x, d = 0)

The d-th derivative of the j-th basis function of b at x, or of the spline $u_h = \sum_j \hat{u}_j \phi_j$ with coefficients û.

x may also be a vector, in which case a vector is returned.

The explicit derivative order is what the lazy b' products cannot give: the discrete brackets need d up to 3, and stacking Derivative that deep does not compose.

julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3);

julia> evaluate(b, 1, 0.1) ≈ b[0.1, 1]
true

julia> evaluate(b, nbasis(b), 1.0)         # interpolatory at the right endpoint
1.0

On a PeriodicBSplineBasis any real x is accepted and reduced onto the domain first, so the result is the periodic extension. On the other two, x outside the domain gives zero — which is the mathematically correct value for a compactly supported basis, and is worth being aware of in a particle method, where a particle that leaves the domain silently stops contributing rather than raising an error.

This is the reference path, one basis function at a time, and it runs the recursion as written rather than through a faster equivalent: one value costs $O(2^p)$. A loop that needs every nonzero function at a point — a particle deposition, or a matrix assembly — should use evaluate_all instead, which is $O(p^2)$ for the whole block of p+1 together. The two are within a third of each other at $p = 3$ and a factor of 150 apart at $p = 12$; see scripts/evaluate_cost_scaling.jl.

evaluate(B::TensorProductBasis, I, x, d = ntuple(_ -> 0, D))
evaluate(B::TensorProductBasis, û::AbstractArray, x, d = ntuple(_ -> 0, D))

The mixed derivative $\partial_1^{d_1} \cdots \partial_D^{d_D}$ at x of the basis function indexed by I, or of the spline with coefficient array û.

I may be a CartesianIndex, a tuple, or a linear index into LinearIndices(B). x is any indexable D-vector. d is a per-axis tuple of derivative orders; the gradient component $\partial_k$ is d = ntuple(i -> i == k ? 1 : 0, D).

julia> B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(8, 0 .. 1), 2);

julia> evaluate(B, (3, 4), (0.3, 0.4)) ≈
       evaluate(bases(B)[1], 3, 0.3) * evaluate(bases(B)[2], 4, 0.4)
true

julia> û = zeros(size(B)...); û[3, 4] = 1.0;

julia> evaluate(B, û, (0.3, 0.4)) ≈ evaluate(B, (3, 4), (0.3, 0.4))
true

Evaluating the spline this way costs $O(N)$; evaluating it at a point where only $\prod_d (p_d+1)$ terms are nonzero costs that much instead, and is what evaluate_all! is for.

evaluate(B::PolarSplineBasis, k::Integer, x, d = (0, 0))
evaluate(B::PolarSplineBasis, û::AbstractVector, x, d = (0, 0))

The mixed derivative $\partial_s^{d_1} \partial_\theta^{d_2}$ at x = (s, θ) of the k-th basis function, or of the spline with coefficient vector û.

julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3),
                            PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3));

julia> sum(evaluate(B, k, (0.3, 1.1)) for k in eachindex(B)) ≈ 1     # partition of unity
true

julia> maximum(abs, [evaluate(B, 1, (0.0, θ)) - 1/3 for θ in range(0, 2π, 17)]) < 1e-15
true

The k-th basis function is evaluated through the $k$-th column of recombination_matrix, so a pole function costs $O(N_\theta)$ and every other one costs the same as the tensor-product function it is. The spline is not: it goes through parent_coefficients and then the parent's local block, which is $O(\prod_d (p_d+1))$ per point — but it forms R * û on every call, so a sweep over many points should pass the whole vector of them, which lifts that out of the loop.

evaluate(s::Spline, x, d = 0)

The d-th derivative of s at x. For a tensor product d is a per-axis tuple.

evaluate_all(b::AbstractBSplineBasis, x, d = 0)

The allocating form of evaluate_all!: returns (j₀, values) with a freshly allocated buffer of local_width entries.

Convenient at a call site that runs once. A particle loop should hold its own buffer and call evaluate_all!, which allocates nothing.

evaluate_all(B::PolarSplineBasis, x, d = (0, 0))

The basis functions that are nonzero at x = (s, θ), as a pair (idx, values) of the global indices and the matching values of $\partial_s^{d_1} \partial_\theta^{d_2} \Psi$.

Unlike the tensor-product and one-dimensional forms, which return the first index of a contiguous block, this returns the indices themselves. The polar index set is not a product and the block is not contiguous in it: in the first two radial cells the three pole functions are nonzero alongside the outer rows, and they sit at the front of the index set rather than beside them. An offset cannot say that, so the indices are given.

julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3),
                            PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3));

julia> idx, vals = evaluate_all(B, (0.02, 1.1));

julia> idx[1:3], sum(vals) ≈ 1
([1, 2, 3], true)

julia> idx, vals = evaluate_all(B, (0.9, 1.1));      # away from the pole: 4 × 4 functions

julia> length(idx), sum(vals) ≈ 1
(16, true)

evaluate_all(b::AbstractBSplineBasis, x, d = 0)
evaluate_all!(values, b::AbstractBSplineBasis, x, d = 0)

The d-th derivatives at x of the basis functions that do not vanish there, returned as (j₀, values) — or, for the in-place form, written into values with j₀ returned.

values must have local_width(b) entries — p+1 for a clamped or periodic basis, and possibly more for a recombined one; values[t] is the derivative of basis function basis_index(b, j₀ + t - 1).

julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3);

julia> j₀, v = evaluate_all(b, 0.3);

julia> j₀, length(v)
(3, 4)

julia> all(v[t] ≈ evaluate(b, j₀ + t - 1, 0.3) for t in eachindex(v))
true

julia> sum(v) ≈ 1                     # the other N - p - 1 functions vanish at 0.3
true

This is the routine a particle method needs. Depositing $N_p$ particles onto the basis costs $O(N_p \, p^2)$ through this and $O(N_p \, N \, 2^p)$ through evaluate one index at a time — the difference between a loop over the p+1 functions that actually overlap the particle and a loop over the whole basis, and, within each term of it, between the triangular scheme and the unmemoised recursion.

Outside the domain of a bounded basis values is filled with zeros, matching evaluate one index at a time, so the deposition above adds nothing for a particle that has left the domain rather than depositing an extrapolated polynomial. On a PeriodicBSplineBasis every real x is reduced onto the domain first and no point is outside.

Method

De Boor's triangular scheme evaluates the q+1 nonzero splines of degree $q = p - d$ on the span containing x in $O(q^2)$, and the derivative recursion

\[D^{m+1} \phi_j^r = r \left( \frac{D^m \phi_j^{r-1}}{x_{j+r} - x_j} - \frac{D^m \phi_{j+1}^{r-1}}{x_{j+r+1} - x_{j+1}} \right)\]

then lifts the block from degree $p-d$ to degree $p$, one order at a time, the block growing by one function at each step. The internal _bspline computes the same numbers from the recursion written out as it stands, and the test suite checks the two against each other at every degree and derivative order.

BSplineDerivative

The type of Derivative(axes(b,1)) * b for an AbstractBSplineBasis b, equivalently of b'.

A lazy product: it stores the basis and runs the derivative recursion on indexing. For derivatives of order higher than one use evaluate with an explicit d.

PeriodicBSplineDerivative

The BSplineDerivative of a PeriodicBSplineBasis specifically.

A compatibility alias. Dispatch on BSplineDerivative instead, which covers all three bases.

SplineQuadrature(basis; nq = quadrature_order(degree(basis)), dmax = 3)

The assembly table of an AbstractBSplineBasis: its basis functions and their derivatives tabulated at the global Gauß-Legendre quadrature points, together with the quadrature weights and the mass matrix.

Any of the three bases serves — clamped, periodic or recombined. What differs between them is which functions are nonzero on a cell (local_indices) and which representation the mass matrix takes (mass_operator); the contractions below are the same in all three cases.

This is the one data structure every assembly here is built from. With Φ[d+1][i,q] the d-th derivative of $\phi_i$ at the q-th quadrature point and w the weight vector, every matrix of the form $\int_\Omega f(x) \, D^a \phi_k \, D^b \phi_l \, dx$ is a single weighted contraction,

\[A = \Phi_a \, \mathrm{diag}(f \odot w) \, \Phi_b^T ,\]

which is what mixed_matrix and weighted_matrix do. Writing the assembly this way rather than element-by-element is what makes a variable coefficient — the $\int u_h \phi_k \phi_l'$ of the second KdV bracket, say — cost no more than a constant one.

Arguments

• nq: quadrature points per cell. The default integrates degree $3p-1$ exactly; see quadrature_order.
• dmax: the highest derivative order tabulated. Three is what a third-order operator needs after one integration by parts.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3));

julia> S = derivative_matrix(q);

julia> maximum(abs, S + S') < 1e-14        # antisymmetric on a periodic mesh
true

Storage

Φ is a SparseMatrixCSC, N by n * nq, one per derivative order. Only the entries inside each basis function's support of p+1 cells are ever nonzero, and only those are stored: the structurally nonzero count per row is (p+1) * nq rather than n * nq, which is what keeps the contraction $\Phi \, \mathrm{diag}(f \odot w) \, \Phi^T$ proportional to N p instead of N².

quadrature_order(p)

The number of Gauß-Legendre points per cell that integrates degree $3p-1$ exactly, $n_q = \lceil 3p/2 \rceil$.

An nq-point Gauß-Legendre rule is exact to degree $2 n_q - 1$, so this is the smallest nq with $2 n_q - 1 \ge 3p - 1$. Degree $3p-1$ is what the consistency of a Galerkin discretisation of a quadratic nonlinearity needs — an identity such as $\int 6 u_h u_{h,x} \phi_i = -3 \int u_h^2 \phi_i'$ pairs three basis functions and one derivative. The mass matrix alone would need only $2p-1$.

What this rule is not needed for is antisymmetry. A bracket assembled in explicitly skew-symmetrised form is antisymmetric to the last bit at any nq and on any mesh; it is accuracy, not structure, that the quadrature buys.

julia> quadrature_order.(1:4)
4-element Vector{Int64}:
 2
 3
 5
 6

quadrature_nodes(q::SplineQuadrature)

The global quadrature points, nq Gauß-Legendre points in each of the n cells, concatenated in cell order.

quadrature_nodes(q::TensorProductQuadrature)

The per-axis node vectors, as a tuple. The D-dimensional grid is their product; a point of it is (x[1][q1], x[2][q2], …).

quadrature_nodes(q::PolarSplineQuadrature)

The per-axis node vectors, as a tuple (s, θ), exactly as for the parent TensorProductQuadrature. The two-dimensional grid is their product and the flattening runs the radial axis fastest.

quadrature_weights(q::SplineQuadrature)

The global quadrature weights, scaled by the cell widths, so that sum(quadrature_weights(q)) == L.

quadrature_weights(q::TensorProductQuadrature)

The per-axis weight vectors, as a tuple. The weight of a D-dimensional point is the product of the per-axis weights, which is what quadrature_sample applies.

quadrature_weights(q::PolarSplineQuadrature)

The flattened weight vector of the two-dimensional grid, kron(w_θ, w_s), of length prod(quadrature_grid_size(parent(q))).

This is a vector where the parent's is a tuple of per-axis vectors: the polar tabulation has one column per grid point, so the weight that goes with it is one number per grid point too.

basis_values(q::SplineQuadrature, d = 0)

The table Φ[i,r] of the d-th derivative of $\phi_i$ at the r-th quadrature point.

basis_values(q::PolarSplineQuadrature, d::NTuple{2,Int})
basis_values(q::PolarSplineQuadrature, d::Integer = 0)

The table $\Phi_d[k,r] = \partial_s^{d_1} \partial_\theta^{d_2} \Psi_k(x_r)$ over the flattened quadrature grid, sparse, formed on first use and memoised.

d is a per-axis multi-index: (1,0) is $\partial_s$, (0,1) is $\partial_\theta$. The scalar form is accepted only for d = 0, where the multi-index is unambiguous; a scalar d ≥ 1 names no derivative on a two-dimensional space and is rejected rather than resolved to one of the axes.

basis(q::SplineQuadrature)

The AbstractBSplineBasis the quadrature was built for.

basis(q::TensorProductQuadrature)

The TensorProductBasis the quadrature was built for.

basis(q::PolarSplineQuadrature)

The PolarSplineBasis the quadrature was built for.

basis(s::Spline)

The basis the spline is expanded in.

mass_matrix(op::MassOperator)

The assembled mass matrix behind the operator.

mass_matrix(q::SplineQuadrature)

The mass matrix $\mathbb{M}_{kl} = \int_\Omega \phi_k \phi_l \, dx$, symmetric positive definite.

Assembled once when the SplineQuadrature is built and returned by reference; see mass_factorization for the Cholesky factor that goes with it.

mass_matrix(op::KroneckerMass)

The assembled Kronecker product, kron(M_D, …, M_1).

Formed on demand and not stored: it is what the representation exists to avoid, and at three dimensions it will not fit. Provided so that a test can check the factored solve against the dense one at a size where both are possible.

mass_matrix(q::TensorProductQuadrature)

The assembled Kronecker product kron(M_D, …, M_1), formed on demand — see mass_matrix(op::KroneckerMass) for why it is not stored.

mass_matrix(q::PolarSplineQuadrature)

The mass matrix $\mathbb{M}_{kl} = \int \Psi_k \Psi_l \, ds \, d\theta$, symmetric positive definite, assembled when the quadrature is built and returned by reference.

mass_factorization(q::SplineQuadrature)

The Cholesky factorization of mass_matrix, for solving with the mass matrix without refactorizing.

stiffness_matrix(q::SplineQuadrature)

The stiffness matrix $\mathbb{K}^1_{kl} = \int_\Omega \phi_k' \phi_l' \, dx$, symmetric positive semi-definite with the constants in its kernel.

stiffness_matrix(q::PolarSplineQuadrature)

The matrix $\int \nabla \Psi_k \cdot \nabla \Psi_l \, ds \, d\theta$, the sum over the two axes of the parameter square, symmetric and positive semi-definite with the constants in its kernel on a free space; a homogeneous-Dirichlet rim removes the constant from the space and the matrix is positive definite instead.

This is the gradient of the parameter square, not of the mapped domain: the metric of the map belongs in the weight, as it does for mixed_matrix.

derivative_matrix(q::SplineQuadrature)

The matrix $S_{kl} = \int_\Omega \phi_k \phi_l' \, dx$.

On a periodic mesh this is already antisymmetric, with no skew-symmetrisation needed: $\int_\Omega \partial_x (\phi_k \phi_l) \, dx = 0$ because the basis is $\mathcal{C}^{p-1}$ across the seam and there are no boundary terms. Its skew part is therefore S itself and not S/2, which is where the factor of one half in the first discrete bracket comes from.

The proviso is that the quadrature integrate that total derivative, of degree $2p-1$, exactly — which means nq ≥ p, one point per cell more than quadrature_order would suggest is needed at low degree. Below it the identity fails outright rather than gracefully: at p = 3 and nq = 2 the defect is $9 \times 10^{-3}$ on a RandomMesh. On a UniformMesh it holds at any nq, the assemblies being circulant, so a check run only there confirms it for the wrong reason.

mixed_matrix(q::SplineQuadrature, a, b)

The matrix $\int_\Omega D^a \phi_k \, D^b \phi_l \, dx$.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3));

julia> K0 = mixed_matrix(q, 1, 2);      # ∫ φ_k' φ_l''  =  -∫ φ_k φ_l'''

julia> maximum(abs, K0 + mixed_matrix(q, 0, 3)) < 1e-10
true

mixed_matrix(q::PolarSplineQuadrature, a::NTuple{2,Int}, b::NTuple{2,Int})

The matrix $\int D^a \Psi_k \, D^b \Psi_l \, ds \, d\theta$ over the parameter square, memoised.

The integral is against the parameter measure. A polar space is used for a mapped domain, where the measure carries the Jacobian of the map; that is weighted_matrix with the Jacobian as its coefficient, and it is the caller's to supply because the map is.

weighted_matrix(q::SplineQuadrature, f, a, b)

The matrix $\int_\Omega f(x) \, D^a \phi_k \, D^b \phi_l \, dx$.

f may be a function of the coordinate, or a vector already sampled at quadrature_nodes — the second form is what a variable coefficient given as a spline expansion becomes, and it avoids resampling a field that is already in hand.

The element type is the promotion of the quadrature's with f's, so a complex coefficient field gives a complex matrix — unlike mixed_matrix, which carries no weight and is always eltype(q).

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(16, 2π), 3));

julia> A = weighted_matrix(q, sin, 0, 1);   # ∫ sin(x) φ_k φ_l'

julia> size(A)
(16, 16)

weighted_matrix(q::PolarSplineQuadrature, f, a::NTuple{2,Int}, b::NTuple{2,Int})

The matrix $\int f(x) \, D^a \Psi_k \, D^b \Psi_l \, ds \, d\theta$, with f either a function of the coordinate pair (s, θ) or a vector already sampled on the flattened quadrature grid.

Not memoised, unlike mixed_matrix: the coefficient of a mapped assembly or of a metric bracket depends on the state and changes at every Newton iteration.

basis_integrals(q::SplineQuadrature)

The vector $\int_\Omega \phi_i \, dx$.

This is the gradient of the total mass $C_0 = \int_\Omega u \, dx$ with respect to the degrees of freedom, and it equals $\mathbb{M} \mathbf{1}$ because the basis is a partition of unity. It spans the kernel of the first discrete bracket, which is why the mass is a Casimir there.

Assembled once when the SplineQuadrature is built and returned by reference, as mass_matrix is — do not mutate the result.

basis_integrals(q::PolarSplineQuadrature)

The vector $\int \Psi_k \, ds \, d\theta$.

On a free space, equal to $\mathbb{M} \mathbf{1}$ because the polar basis is a partition of unity there — the pole triangle included, which is what the choice of vertex radius buys. A homogeneous-Dirichlet rim breaks the partition of unity (see polynomial_reproduction), so the two disagree there. Assembled once and returned by reference; do not mutate the result.

l2_projection(q::SplineQuadrature, f)

The coefficients of the $L^2$ projection of f onto the spline space, $\hat{u} = \mathbb{M}^{-1} \int_\Omega f \phi_i \, dx$.

f may be a function or a vector of values at quadrature_nodes.

julia> q = SplineQuadrature(PeriodicBSplineBasis(UniformMesh(32, 2π), 3));

julia> û = l2_projection(q, sin);

julia> abs(evaluate(basis(q), û, 1.0) - sin(1.0)) < 1e-5
true

l2_projection(q::TensorProductQuadrature, f)
l2_projection!(û, q::TensorProductQuadrature, f)

The coefficient array of the $L^2$ projection of f onto the tensor-product spline space, $\hat{u} = \mathbb{M}^{-1} \int_\Omega f \, \Phi_I \, dx$.

f may be a function of a D-tuple of coordinates, or an array already sampled on the quadrature grid — see quadrature_sample. The result has size size(basis(q)).

julia> B = BSplineBasis(UniformMesh(16, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(16, 0 .. 1), 3);

julia> q = TensorProductQuadrature(B);

julia> û = l2_projection(q, x -> x[1]^2 * x[2]);

julia> abs(evaluate(B, û, (0.37, 0.62)) - 0.37^2 * 0.62) < 1e-12
true

The last check is exact to round-off rather than approximate, because a clamped basis of degree 3 reproduces $x^2 y$ exactly — which is polynomial_reproduction, and the property the conservation of mass, momentum and energy rests on.

Both the load assembly and the solve exploit the Kronecker structure: the first is D sparse contractions (contract), the second D one-dimensional solves (KroneckerMass). Neither forms the $N \times N$ mass matrix.

l2_projection(q::PolarSplineQuadrature, f)
l2_projection!(û, q::PolarSplineQuadrature, f)

The coefficient vector of the $L^2$ projection of f onto the polar spline space, $\hat{u} = \mathbb{M}^{-1} \int f \, \Psi_k \, ds \, d\theta$.

f may be a function of the coordinate pair (s, θ), or a vector already sampled on the flattened quadrature grid.

julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(16, 0 .. 1), 3),
                            PeriodicBSplineBasis(UniformMesh(32, 0 .. 2π), 3));

julia> q = PolarSplineQuadrature(B);

julia> û = l2_projection(q, x -> 1.0);         # the constants are in the space exactly

julia> abs(evaluate(B, û, (0.0, 0.7)) - 1) < 1e-12
true

The load is one sparse matrix-vector product against the tabulation, and the solve is the sparse Cholesky of mass_operator. Neither has the Kronecker shortcut of a tensor-product projection, and that is the cost of the pole.

l2_projection!(û, q::SplineQuadrature, f)

In-place l2_projection, writing the coefficients into û.

This allocates nothing on every basis but one: the $f \odot w$ product goes into a buffer held by the quadrature, the load vector is formed with mul! straight into û, and both the CirculantMass and the BandedMass solve are themselves allocation-free. The exception is a periodic basis on a non-uniform mesh, where the FactorizedMass solve still allocates a CHOLMOD temporary, as mass_solve! notes.

A sample whose element type is wider than the quadrature's — a complex f — gets its own product instead of being narrowed into that buffer, so this method accepts exactly what l2_projection accepts and merely stops being allocation-free there.

The shared buffer is one of the two reasons a SplineQuadrature may not be used from two threads at once; see the warning there.

l2_projection!(s::Spline, q, f)

Project f onto the space of s, writing the coefficients into the array s already holds so that every reference to s sees the new function.

MassOperator{T}

The mass matrix in the form the solves actually want: something that can be applied and, above all, inverted, without forming $\mathbb{M}^{-1}$.

Three representations are provided, and which one is built is decided by the basis:

• CirculantMass for a periodic basis on a UniformMesh, where the basis functions are translates of a single cardinal spline and $\mathbb{M}$ is circulant, hence diagonalised by the discrete Fourier transform. A solve is two transforms and a pointwise multiplication, $O(N \log N)$, with the transforms planned once.
• BandedMass for a bounded basis. The overlaps are contiguous — there is no seam — so the matrix is banded outright and a banded Cholesky solves it in $O(Np)$ with no allocation at all.
• FactorizedMass for a periodic basis on a GradedMesh or a RandomMesh, where the matrix is banded modulo $N$ but not circulant — the basis functions are no longer translates of one another — so there is nothing for a Fourier transform to diagonalise, and the wrap-around entries put it outside the banded representation too. A sparse Cholesky factorisation is what is left.

All three answer \, ldiv! and Matrix.

A mass matrix is positive definite, but the same three representations carry the other assemblies of a SplineQuadrature, and a stiffness matrix on a periodic basis is singular: the constants the basis represents lie in its kernel. The kernel keyword says what to do about that. :reject is the default and throws, since a singular assembly is usually too coarse a quadrature rather than an intended one. :project states that the kernel is the constants and asks for the solution that has no constant component — the mean-free solution — which is what $-\phi'' = \rho$ on a periodic domain asks for. Both periodic representations implement it, each in the way its own structure allows, and both verify the assertion rather than trusting it. Neither forms the rank-one shift $\mathbb{M} + \mathbb{1}\mathbb{1}^T/N$ that would remove the singularity by filling the matrix in completely.

CirculantMass(M, n; kernel = :reject, rtol = sqrt(eps(T)))
CirculantMass(c::AbstractVector, n; kernel = :reject)

The mass matrix of a uniform periodic mesh, represented by the eigenvalues of its circulant structure and a pair of planned real transforms.

On a uniform mesh every basis function is a translate of one cardinal spline, so $\mathbb{M}_{kl}$ depends only on $k - l \bmod N$ and

\[\mathbb{M} = F^{*} \operatorname{diag}(\hat{c}) F , \qquad \hat{c} = \mathcal{F}(\mathbb{M}_{:,1}) ,\]

with $F$ the discrete Fourier transform. A solve is therefore a forward transform, a pointwise multiplication by the reciprocals of $\hat{c}$, and an inverse transform.

The plans are created once, at construction, and the spectral buffer is preallocated, so a solve allocates nothing beyond its result. \ allocates the result; mass_solve! does not.

The first column is read from the assembled matrix rather than recomputed, and the construction checks that the matrix really is circulant — a silent mismatch here would give wrong answers on every mesh that is uniform by accident rather than by construction. rtol is the tolerance of that check, relative to the largest entry of the first column, so that the verdict survives a rescaling of the assembly and holds in every element type.

The second form takes the first column itself, wrapped in a Circulant. That is the whole matrix, because a circulant matrix is its first column, and it is what a caller that can describe its operator in $O(n)$ should pass: the matrix form makes it build $n^2$ entries for this constructor to read $n$ of them back out, and then keeps them for the life of the operator. The circulance check has nothing to verify on it, so the path is $O(n)$ in construction as well as in storage, and rtol is accepted only by the matrix form because only there is there anything to compare. mass_matrix of such an operator returns the Circulant; Matrix of it materialises.

kernel = :project accepts a matrix that is singular with the constants in its kernel, and gives the constant mode a zero factor instead of an infinite one. That is the Moore–Penrose pseudoinverse: the solve drops the constant component of the right-hand side and returns the solution that has none of it. The transform does the whole of the work, so nothing is added to the matrix and nothing is taken out of it. This is what makes a periodic stiffness matrix usable here, and it is the rule an FFT Poisson solver applies when it sets the $k = 0$ factor to zero rather than dividing by it.

The deflation lives entirely in the stored reciprocal eigenvalues, so the type carries no kernel-mode parameter and there is one mass_solve! for both modes. FactorizedMass carries such a parameter because there the two modes are two different solves, and the parameter is what selects between them.

Circulant(c)

The $n \times n$ circulant matrix whose first column is c, stored as that column alone: $C_{ij} = c_{(i - j) \bmod n + 1}$.

A representation, not an arithmetic type. It exists so that a caller holding a circulant matrix it can describe in $O(n)$ — a periodic assembly, or such an assembly plus a rank-one shift — can hand CirculantMass that description instead of materialising $n^2$ entries for it to read one column back out of. getindex, size and Matrix are what it provides; Matrix materialises, and is the way out when something downstream really does want the whole matrix.

The column is copied, so a later mutation of the caller's vector cannot leave the operator's stored eigenvalues describing a different matrix from the one mass_matrix reports.

julia> C = SimpleSplines.Circulant([2.0, 1.0, 0.0, 1.0]);

julia> size(C), C[1, 1], C[2, 1], C[1, 2]
((4, 4), 2.0, 1.0, 1.0)

julia> Matrix(C)
4×4 Matrix{Float64}:
 2.0  1.0  0.0  1.0
 1.0  2.0  1.0  0.0
 0.0  1.0  2.0  1.0
 1.0  0.0  1.0  2.0

BandedMass(M)

The banded Cholesky factorisation of the mass matrix, for the bases whose mass matrix is genuinely banded rather than banded modulo $N$.

A basis function overlaps only the $2p+1$ others whose supports meet its own, and on a bounded basis those are contiguous — there is no seam, so no corner entries. The factorisation is therefore $O(Np^2)$ with no fill-in analysis and no reordering, and, unlike CHOLMOD, the factor answers ldiv! in place: a solve allocates nothing at all. That is what makes l2_projection! allocation-free on every bounded basis rather than only on the uniform periodic one.

A non-contiguous argument, which LAPACK cannot address, is staged through a buffer the operator owns, so one operator is not to be shared between threads.

FactorizedMass remains for the periodic non-uniform case, where the matrix wraps.

FactorizedMass(M; kernel = :reject)

The sparse Cholesky factorisation of the mass matrix, for meshes on which it is not circulant.

$\mathbb{M}$ is symmetric positive definite and banded modulo $N$ — a basis function overlaps only the $2p+1$ others whose supports meet its own — so the factorisation is cheap and the solve is $O(Np)$.

With kernel = :project the matrix may be singular instead, with the constants in its kernel. What is factorised is then the minor that drops the first degree of freedom, which is positive definite exactly when the kernel is the constants and nothing more: a vector supported away from the first index lies in the span of $\mathbb{1}$ only if it is zero. The dropped degree of freedom is the gauge, and mass_solve! fixes it afterwards by taking the mean out. Deleting a row and a column preserves the sparsity; the rank-one shift $\mathbb{M} + \mathbb{1}\mathbb{1}^T/N$, which also removes the singularity, makes the matrix structurally full and the factorisation $O(N^2)$.

KroneckerMass(ops...)

The mass operator of a TensorProductBasis, $\mathbb{M} = \mathbb{M}^{(D)} \otimes \dots \otimes \mathbb{M}^{(1)}$, represented by its D one-dimensional factors and never assembled.

A solve applies each factor's inverse along its own axis:

\[\mathbb{M}^{-1} = \left( \mathbb{M}^{(D)} \right)^{-1} \otimes \dots \otimes \left( \mathbb{M}^{(1)} \right)^{-1} ,\]

which is exact — this is an identity, not an approximate splitting — and costs $\sum_d (N/N_d)$ one-dimensional solves. Each factor keeps whatever representation it had: a CirculantMass on a periodic uniform axis, a FactorizedMass otherwise.

\ and ldiv! take and return arrays of size size(B). A vector of length length(B) is also accepted and is reshaped, on the convention that the first axis varies fastest — the same convention LinearIndices(B) uses, and the one that makes Matrix(op) equal kron(M_D, …, M_1) rather than the reverse.

julia> B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(6, 0 .. 1), 2);

julia> q = TensorProductQuadrature(B);

julia> op = mass_operator(q);

julia> b = randn(size(B)...);

julia> maximum(abs, Matrix(op) * vec(b) - vec(op * b)) < 1e-10
true

julia> maximum(abs, op \ (op * b) - b) < 1e-10
true

mass_operator(M, basis; kernel = :reject)

Build the MassOperator appropriate to basis: a CirculantMass for a periodic basis on a UniformMesh, a FactorizedMass for a periodic basis on any other mesh, and a BandedMass for a bounded basis — clamped or recombined — whose mass matrix has no seam to wrap across.

kernel = :project passes the deflation on to whichever of the two periodic representations is chosen, so that a caller with a singular assembly says what it wants rather than which representation implements it. See MassOperator for what the deflation means and each representation for how it is done.

mass_operator(c::AbstractVector, b::PeriodicBSplineBasis; kernel = :reject)

Build the CirculantMass of the circulant matrix whose first column is c.

This is the $O(n)$ entry point for a caller that can describe its operator without materialising it. A periodic stiffness matrix shifted by the rank-one mean projector is the case it exists for: $S + \mathbb{1}\mathbb{1}^T/n$ is circulant like $S$, with first column S[:, 1] .+ inv(n), but forming it fills an $O(n)$ sparse matrix into an $O(n^2)$ dense one.

S = stiffness_matrix(SplineQuadrature(b))
mass_operator(S[:, 1] .+ inv(size(S, 1)), b)      # O(n); the matrix form is O(n²)

Only a UniformMesh is accepted. A first column describes the whole matrix only when the matrix is circulant, and on any other mesh it is not — so this raises rather than falling through to FactorizedMass, which would read the vector as something it is not.

mass_operator(q::SplineQuadrature)

The MassOperator of the quadrature, chosen by the basis rather than by the mesh: a CirculantMass for a periodic basis on a UniformMesh, a FactorizedMass for a periodic basis on any other mesh, and a BandedMass for a bounded basis — clamped or recombined — whose mass matrix has no seam to wrap across.

mass_operator(q::TensorProductQuadrature)

The KroneckerMass built from the per-axis mass operators, never assembled.

mass_operator(q::PolarSplineQuadrature)
mass_factorization(q::PolarSplineQuadrature)

The FactorizedMass of the polar mass matrix — a sparse Cholesky, not a KroneckerMass; see PolarSplineQuadrature for why the Kronecker structure is not available.

mass_factors(op::KroneckerMass)

The tuple of one-dimensional MassOperators, axis order.

mass_solve!(y, op::MassOperator, x)

Solve $\mathbb{M} y = x$ in place. y and x may alias.

Attached to the function rather than to either method, so that the two representations of MassOperator share one piece of documentation and a cross-reference to the name resolves.

TensorProductBasis(bases...)
b₁ ⊗ b₂ ⊗ …

The tensor-product basis $\Phi_{i_1 \dots i_D}(x) = \prod_{d=1}^{D} \phi^{(d)}_{i_d}(x_d)$ built from D one-dimensional bases.

Each factor is an independent AbstractBSplineBasis, so degree, mesh, domain and boundary condition are per-axis. Nothing is shared between the axes, and nothing needs to be:

julia> B = TensorProductBasis(
           BSplineBasis(UniformMesh(32, 0 .. 2π), 3, Periodic()),      # x: cubic, periodic
           BSplineBasis(UniformMesh(41, -10 .. 10), 4, Free()),        # v: quartic, clamped
       );

julia> ndims(B), size(B), nbasis(B)
(2, (32, 45), 1440)

julia> degree(B)
(3, 4)

julia> domain(B)
(0.0 .. 6.28319) × (-10.0 .. 10.0)

julia> leftendpoint(domain(bases(B)[2])), rightendpoint(domain(bases(B)[2]))
(-10.0, 10.0)

Coefficients are an array, not a vector

A spline in this basis is $u_h = \sum_I \hat{u}_I \Phi_I$ with û a D-dimensional array of size size(B). That is the natural shape: the Kronecker structure of every operator is visible in it, and LinearIndices/CartesianIndices do the index arithmetic that would otherwise be written by hand — one of the places a hand-rolled tensor-product spline reliably goes wrong, since the flattening convention has to agree between the evaluation, the mass matrix and the projection.

The Kronecker structure is used, not just noted

The mass matrix of a tensor-product basis is $\mathbb{M} = \mathbb{M}^{(D)} \otimes \dots \otimes \mathbb{M}^{(1)}$, and it is never formed. A solve is D one-dimensional solves applied along each axis in turn — see KroneckerMass — which is $O(N \sum_d p_d)$ against the $O(N^3)$ of a dense factorisation of the Kronecker product, and needs $O(\sum_d N_d^2)$ storage rather than $O(N^2)$. At the 41-element cubic basis of a two-dimensional velocity space that is a $1681 \times 1681$ dense Cholesky avoided; in three dimensions it is a $68921^2$ one, which does not fit.

The load vector of an $L^2$ projection factorises the same way, as D successive sparse contractions — see l2_projection. The integrand itself need not be separable; only the basis is, and that is enough.

See also TensorProductQuadrature for the assembly, evaluate_all! for the local evaluation a particle loop needs, and polynomial_reproduction, which is the minimum over the axes.

⊗(b₁, b₂)

Tensor product of two bases, or of a tensor product and a basis — b₁ ⊗ b₂ ⊗ b₃ associates to a single flat TensorProductBasis of three factors rather than a nest of two.

bases(B::TensorProductBasis)

The tuple of one-dimensional bases the product is built from.

TensorProductQuadrature(B; nq = quadrature_order.(degree(B)), dmax = 3)

The assembly table of a TensorProductBasis: D one-dimensional SplineQuadratures, one per axis, together with the KroneckerMass built from their mass operators.

nq and dmax may be given per axis as a tuple, or once for all axes.

There is no D-dimensional tabulation. Everything a Galerkin assembly on a tensor-product basis needs is a sequence of contractions with the one-dimensional tabulations, and holding those is $O(\sum_d N_d n_d n_{q,d})$ rather than the $O(N \prod_d n_d n_{q,d})$ a D-dimensional table would cost — which at three dimensions is the difference between kilobytes and not fitting.

julia> B = BSplineBasis(UniformMesh(16, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(12, 0 .. 2), 2);

julia> q = TensorProductQuadrature(B);

julia> size(quadrature_nodes(q)[1]), size(quadrature_nodes(q)[2])
((80,), (36,))

julia> û = l2_projection(q, x -> sin(π * x[1]) * x[2]);

julia> abs(evaluate(B, û, (0.3, 1.1)) - sin(π * 0.3) * 1.1) < 1e-3
true

quadratures(q::TensorProductQuadrature)

The tuple of one-dimensional SplineQuadratures, one per axis, in axis order.

This is how the per-axis tables are reached: basis_values(quadratures(q)[k], d) is the tabulation of axis k, and there is no D-dimensional table to ask for.

quadrature_grid_size(q::TensorProductQuadrature)

The size of the D-dimensional quadrature grid, (n_d * nq_d) per axis.

quadrature_grid_size(q::PolarSplineQuadrature)

The per-axis size of the quadrature grid, (n_d * nq_d) per axis.

quadrature_sample(q::TensorProductQuadrature, f)

Sample f on the D-dimensional quadrature grid, returning an array of size quadrature_grid_size.

f is called with a D-tuple of coordinates. This is the array l2_projection contracts, and it is exposed because the integrand of interest is often not a plain function of position — the $\mathbb{L}_k$ of a metriplectic collision operator is $\int \varphi_i \, (1 + \log f_s)$, whose integrand is built from a spline that is already in hand:

fs = [evaluate(B, f̂, x) for x in Iterators.product(quadrature_nodes(q)...)]
L̂  = l2_projection(q, 1 .+ log.(fs))

contract(q::TensorProductQuadrature, F, d = ntuple(_ -> 0, D))

Contract the array F, given on the quadrature grid, against the one-dimensional basis tabulations, giving the load array

\[L_{i_1 \dots i_D} = \sum_{q_1 \dots q_D} F_{q_1 \dots q_D} \prod_{k=1}^{D} D^{d_k} \phi^{(k)}_{i_k}(x_{q_k}) \, w_{q_k} .\]

F must already include whatever the integrand is; the quadrature weights are applied here.

Method

One axis at a time. With $\Phi_k$ the sparse $N_k \times Q_k$ tabulation of axis k,

\[L = \Phi_1 \times_1 \Phi_2 \times_2 \dots \times_D \, (F \odot w) ,\]

evaluated by reshaping the array so that the axis being contracted is first, multiplying by $\Phi_k$, and cycling that axis to the back. After D steps the axes are back in order. Each step is one sparse matrix-matrix product, so the whole contraction costs $O\!\left(\sum_k N_k p_k \prod_{j \ne k} \cdot\right)$ — linear in the grid size, against the $O(N \prod_k n_k n_{q,k})$ of forming each entry by its own quadrature loop.

PolarSplineBasis(radial, angular)
PolarSplineBasis(B::TensorProductBasis)

The polar spline space on the parameter square $[a,b] \times [c,c+L)$, whose left radial endpoint $s = a$ is a pole: one point of the physical domain, reached from every angle.

The first two rows of the radial basis — $2 N_\theta$ tensor-product functions — are replaced by three, the pole triangle, which together span the constants and the two linear functions of the pseudo-Cartesian chart at the pole. Every other row passes through unchanged, so

\[N = 3 + (N_s - 2) \, N_\theta .\]

julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3),
                            PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3));

julia> nbasis(B), nbasis(parent(B))
(147, 176)

julia> û = zeros(nbasis(B)); û[1:3] .= 1;      # the three pole functions, added up

julia> abs(evaluate(B, û, (0.0, 0.3)) - evaluate(B, û, (0.0, 2.9))) < 1e-15
true

Why this is not a boundary condition

A BoundaryCondition constrains one axis at one end, and BSplineBasis(mesh, p, bc) builds the one-dimensional basis that satisfies it. A pole is not of that kind: single-valuedness at $s = a$ says that the $\theta$-dependence there is constant, which couples the two axes. There is no radial basis whose tensor product with anything is a polar space, so this is a two-dimensional basis type beside TensorProductBasis rather than a fourth boundary condition — and the independent-per-axis promise of the tensor-product layer stays true, because the polar space is not one.

The construction

Write a tensor-product spline as $u = \sum_{ij} \hat{u}_{ij} N_i(s) M_j(\theta)$ with the radial basis clamped of degree $p \ge 2$ and the angular basis periodic. A clamped basis has $N_1(a) = 1$ and $N_i(a) = 0$ for $i > 1$, and, at degree two or more, only $N_1$ and $N_2$ have a nonzero derivative there. So the value and the radial derivative at the pole read

\[u(a,\theta) = \sum_j \hat{u}_{1j} M_j(\theta) , \qquad \partial_s u(a,\theta) = \sum_j \bigl( \hat{u}_{1j} N_1'(a) + \hat{u}_{2j} N_2'(a) \bigr) M_j(\theta) ,\]

and nothing else in the basis can affect either. Imposing

\[\hat{u}_{1j} = c , \qquad \hat{u}_{2j} = c + \frac{\alpha \cos \theta_j + \beta \sin \theta_j}{N_2'(a)} ,\]

with $\theta_j$ the Greville abscissae of the angular basis (nodes), makes $u(a,\theta) = c$ a constant — the partition of unity of the angular basis — and

\[\partial_s u(a,\theta) = \alpha \, C(\theta) + \beta \, S(\theta) , \qquad C = \sum_j \cos\theta_j \, M_j , \quad S = \sum_j \sin\theta_j \, M_j .\]

That is a linear function of the pseudo-Cartesian coordinates

\[\tilde{x} = (s-a) \, C(\theta) , \qquad \tilde{y} = (s-a) \, S(\theta)\]

— see pseudo_cartesian — so $u$ is $C^1$ at the pole in that chart, by construction and not by refinement. The constrained set is three-dimensional, parametrised by $(c, \alpha, \beta)$.

The pole triangle

The three functions that span it are chosen to be the barycentric coordinates of a triangle in the $(\tilde{x}, \tilde{y})$ plane, with vertices

\[v_k = \frac{2}{N_2'(a)} \bigl( \cos \psi_k, \, \sin \psi_k \bigr) , \qquad \psi_k = \frac{2\pi (k-1)}{3} ,\]

so that $\Psi_k$ has the value $1/3$ at the pole and the gradient $\nabla \ell_k = \tfrac{2}{3}\|v\|^{-1}(\cos\psi_k, \sin\psi_k)$. Written out, the coefficients of $\Psi_k$ on the first two rows are

\[\lambda^{(1)}_{kj} = \tfrac{1}{3} , \qquad \lambda^{(2)}_{kj} = \tfrac{1}{3} \bigl( 1 + \cos(\psi_k - \theta_j) \bigr) ,\]

which is why that vertex radius and no other: it is the smallest for which $\lambda^{(2)} \ge 0$, so the basis stays non-negative, and the three vertices sum to zero, so $\sum_k \Psi_k = \sum_j (N_1 + N_2) M_j$ and the whole basis is still a partition of unity. Any other radius spans the same three-dimensional space and is only a worse-conditioned basis of it. See pole_triangle.

Indexing and coefficients

The index set is not a product, so a spline in this basis carries a coefficient vector of length nbasis(B), not an array: 1:3 are the pole functions and 3 + (i-2) + (j-1)*(N_s-2) is the outer function $N_i M_j$, the radial index running fastest as it does in LinearIndices(parent(B)).

Every basis function is a combination of the parent's, $\Psi_k = \sum_I R_{Ik} \Phi_I$, and every assembly is the parent's conjugated by that matrix — exactly as for a RecombinedBSplineBasis one dimension down. See recombination_matrix and parent_coefficients.

Requirements

The radial basis must be clamped and unconstrained at the pole end, of degree at least two; degree one has no $C^1$ to impose, and a periodic radial axis has no pole. A plain BSplineBasis qualifies, and so does a RecombinedBSplineBasis carrying Free at the pole end — see The rim below. A rim condition of order $m$ must also leave the first two radial functions alone, which needs $m + 3$ functions in the parent it recombines; only a condition of order two or more can fail that. At least one radial row must survive the pole triangle, which a clamped basis gives for free and a rim condition can take away.

The angular basis must be a PeriodicBSplineBasis with at least three functions, or $C$, $S$ and the constants are not independent and the triangle is degenerate.

The rim

The pole is not a boundary condition — it couples the two axes, which is why this type exists beside TensorProductBasis rather than as a fourth BSplineBasis(mesh, p, bc) method. The rim, the outer end $s = b$, is an ordinary boundary condition on the radial axis and is imposed where every other one is, by recombining that axis:

radial = RecombinedBSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3), Free(), Dirichlet())
B = PolarSplineBasis(radial, PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3))

The two compose with nothing to reconcile. A rim condition changes only rows the pole triangle does not read: the triangle is built from the first two functions and their derivatives at $s = a$, and a right-recombined basis's first two functions are the clamped parent's, unchanged — as long as the rim block stays clear of them, which is what the $m + 3$ requirement above asks for. So $R$ is built exactly as before, on the smaller $N_s$.

What this costs is the partition of unity: a homogeneous-Dirichlet rim removes the constant from the space by construction, so $\sum_k \Psi_k \equiv 1$ becomes false — near the rim, and only there. The $C^0$ and $C^1$ properties at the pole are untouched, since no function the triangle is built from has changed. polynomial_reproduction of the radial axis says which of the two regimes a basis is in.

PolarSplineQuadrature(B; nq = map(quadrature_order, degree(B)), dmax = 3)

The assembly table of a PolarSplineBasis: the parent's TensorProductQuadrature, together with the tabulations and the mass factorisation of the polar basis itself.

julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3),
                            PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3));

julia> q = PolarSplineQuadrature(B);

julia> size(basis_values(q, (0, 0)))
(147, 3200)

julia> sum(basis_integrals(q)) ≈ 2π       # ∫ Σ_k Ψ_k over the parameter square
true

Why there is no KroneckerMass

The mass matrix of a tensor-product basis is a Kronecker product and is never formed; its solve is one one-dimensional solve per axis. The pole rows destroy that structure — a pole function is a sum over the whole angular axis — so the polar mass matrix is assembled and factorised. It is sparse: the three pole columns are dense in the $2 N_\theta$ parent rows they touch, and every other column has the $(p_s+1)(p_\theta+1)$ neighbours of a tensor-product function. Only the $3 \times 3$ pole block is dense, so the factorisation is that of a banded matrix with three border rows rather than of a dense one.

Storage

The table $\Phi_d[k,r] = \partial_s^{d_1} \partial_\theta^{d_2} \Psi_k(x_r)$ over the flattened quadrature grid is R' * kron(Φ_θ, Φ_s), sparse, formed on first use and memoised. The memo is what makes weighted_matrix affordable inside a Newton loop, where the coefficient changes at every iteration and the tabulation does not. The Kronecker product is in reverse axis order because the flattening runs the radial axis fastest — the convention of parent_coefficients, of vec of a quadrature_sample array, and of the index set itself.

pole(B::PolarSplineBasis)

The radial coordinate of the pole, i.e. the left endpoint of the radial axis.

pole_triangle(B::PolarSplineBasis)

The 3 × 2 matrix of the pole-triangle vertices $v_k$ in the pseudo-Cartesian chart.

The $k$-th pole function restricted to the pole — its value and its gradient there — is the barycentric coordinate of this triangle that is one at $v_k$ and zero at the other two.

pseudo_cartesian(B::PolarSplineBasis, x)

The chart $(\tilde{x}, \tilde{y}) = (s-a) \, (C(\theta), S(\theta))$ in which the polar space is $C^1$ at the pole, with

\[C = \sum_j \cos \theta_j \, M_j , \qquad S = \sum_j \sin \theta_j \, M_j\]

the angular splines whose coefficients are the cosine and sine of the Greville abscissae.

$C$ and $S$ are splines rather than the trigonometric functions themselves, and that is the point: $\cos\theta$ is not in the angular spline space, so a chart built from it would make the $C^1$ property hold only up to the approximation error of that space. Built from $C$ and $S$ it holds exactly. On a free space the space also contains $1$, $\tilde{x}$ and $\tilde{y}$ to round-off; a homogeneous-Dirichlet rim removes them along with the constant (polynomial_reproduction says which regime a basis is in). A geometry map that is itself represented in this space — the isogeometric case — therefore carries the smoothness to the physical domain.

The price is that the unit circle of this chart is the spline through the Greville values of the cosine and the sine, which sits a little inside the true one — some 2.5% at 16 angular cells, a gap that closes at the angular order under refinement. The chart is a chart, not a measurement: distances in it are not distances on the disk.

julia> B = PolarSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3),
                            PeriodicBSplineBasis(UniformMesh(16, 0 .. 2π), 3));

julia> pseudo_cartesian(B, (0.0, 1.2))
(0.0, 0.0)

julia> pseudo_cartesian(B, (1.0, 0.7)) .≈ 2 .* pseudo_cartesian(B, (0.5, 0.7))
(true, true)

julia> x̃, ỹ = pseudo_cartesian(B, (0.5, 0.0)); abs(x̃ - 0.5) < 0.02, abs(ỹ) < 1e-15
(true, true)

parent_coefficients(B::PolarSplineBasis, û)

The coefficient array of û in the parent TensorProductBasis, of size size(parent(B)).

R * û reshaped, which is the one place the flattening convention of the polar index set meets the array shape of the parent. Every evaluation goes through it, and it is exported so that a caller with a spline in hand can reach the parent's own machinery — plotting, a tensor-product projection, evaluate_all! — without rewriting the index arithmetic.

Spline(basis, coefficients)

The spline function $u_h = \sum_I \hat{u}_I \, \phi_I$: a basis together with the coefficients of one element of its span, callable at a point.

The basis may be one-dimensional — any AbstractBSplineBasis, with a vector of coefficients — or a TensorProductBasis, with a D-dimensional array:

julia> b = BSplineBasis(UniformMesh(8, 0 .. 1), 3);

julia> s = Spline(b, l2_projection(SplineQuadrature(b), x -> x^2));

julia> s(0.3) ≈ 0.09, s(0.3, 1) ≈ 0.6            # value and first derivative
(true, true)

julia> B = b ⊗ BSplineBasis(UniformMesh(6, 0 .. 2), 2);

julia> S = Spline(B, l2_projection(TensorProductQuadrature(B), x -> x[1]^2 * x[2]));

julia> S((0.3, 1.4)) ≈ 0.09 * 1.4
true

Why this exists separately from the basis

A basis answers "what is $\phi_j$ here"; a Spline answers "what is $u_h$ here". Keeping them apart is what lets a coefficient array be updated in place — a particle deposition rewrites the coefficients every step while the basis, the quadrature and the mass factorisation stay put — and it is why coefficients returns the array itself rather than a copy. A Spline built on an array shares it, so a projection written into that array is visible through the spline with no rebuild.

The derivative is reached either as a second argument, s(x, d), or as a callable of its own through derivative. The second form is what a right-hand side wants when the derivative appears in a broadcast alongside the value.

SplineDerivative

A derivative of a Spline, as returned by derivative. Callable, so that derivative(s).(v) broadcasts over a vector of points the way s.(v) does.

julia> b = BSplineBasis(UniformMesh(16, 0 .. 1), 3);

julia> s = Spline(b, l2_projection(SplineQuadrature(b), sin));

julia> ds = derivative(s);

julia> maximum(abs, ds.([0.2, 0.5, 0.8]) - cos.([0.2, 0.5, 0.8])) < 1e-5
true

This is the shape a collision operator's right-hand side wants. An expression such as $\nu \, ( f_s'(v_\alpha) / f_s(v_\alpha) + v_\alpha )$ reads as one broadcast over the particle velocities with derivative(fs).(v) ./ fs.(v), and the derivative object holds no state of its own — it shares the coefficient array, so a reprojection is visible through it without rebuilding anything.

coefficients(s::Spline)

The coefficient array itself, not a copy: writing into it changes the spline. That is the point — a projection writes the coefficients of a spline that is already wired into a right-hand side.

derivative(s::Spline, d = 1)
derivative(s::Spline, k::Integer)

The d-th derivative of s as a SplineDerivative — a callable, so that derivative(s).(v) broadcasts over a vector of points the way s.(v) does.

On a TensorProductBasis an integer k selects the $k$-th partial derivative, derivative(s, 2) being $\partial_2 s$; a tuple gives a mixed derivative directly.