Assembly

Everything this package assembles comes out of one data structure. A SplineQuadrature tabulates a basis and its derivatives at the Gauß-Legendre points of every cell, and every matrix of the form

\[A_{kl} = \int_\Omega f(x) \, D^a \phi_k \, D^b \phi_l \, \mathrm{d}x\]

is then a single weighted contraction of that table,

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

using SimpleSplines
using LinearAlgebra
using SparseArrays

Construction

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

Any of the three bases serves. 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 are identical.

b = BSplineBasis(UniformMesh(16, 0 .. 2π), 3, Periodic())
q = SplineQuadrature(b)
nbasis(q), degree(q), ncells(q), length(quadrature_nodes(q))

(16, 3, 16, 80)

• nq is the number of quadrature points per cell. The default, quadrature_order(p) = ⌈3p/2⌉, integrates degree $3p-1$ exactly — see Choosing nq below.
• dmax is the highest derivative order tabulated. Three is what a third-order operator needs after one integration by parts. Asking basis_values for more is an error naming the fix.

quadrature_order.(0:5)

6-element Vector{Int64}:
 0
 2
 3
 5
 6
 8

try
    basis_values(SplineQuadrature(b; dmax = 1), 2)
catch err
    println(err.msg)
end

derivatives up to order 1 were tabulated, but order 2 was requested; rebuild the quadrature with dmax = 2

The table

sum(quadrature_weights(q)) ≈ domainlength(b)

true

The global index of the r-th point of cell k is (k-1)*nq + r, so a hand-written assembly loop reads

Φ = basis_values(q, 0)
w = quadrature_weights(q)
nq = length(quadrature_nodes(q)) ÷ ncells(q)
total = sum(w[(k - 1) * nq + r] * Φ[i, (k - 1) * nq + r]
for k in 1:ncells(q) for i in local_indices(basis(q), k) for r in 1:nq)
total ≈ domainlength(b)              # ∑ᵢ ∫ φᵢ = ∫ 1 = L

true

The table is sparse, and that is the point

A basis function is supported on p+1 cells, so only (p+1)*nq of the n*nq entries in its row can be nonzero, and only those are stored.

nnz(Φ), ncells(q) * (degree(q) + 1) * nq, nnz(Φ) / length(Φ)

(320, 320, 0.25)

Storing $\Phi$ densely would make every contraction $\Phi \, \mathrm{diag}(f \odot w) \, \Phi^{\mathsf T}$ cost $O(N^2 n n_q)$ instead of $O(N p^2 n_q)$.

The matrices

stiffness_matrix and derivative_matrix are mixed_matrix(q, 1, 1) and mixed_matrix(q, 0, 1) — the same objects, not merely equal ones:

stiffness_matrix(q) === mixed_matrix(q, 1, 1),
derivative_matrix(q) === mixed_matrix(q, 0, 1)

(true, true)

mixed_matrix(q, 0, 2) === mixed_matrix(q, 0, 2)

true

The same holds for mass_matrix and basis_integrals, which are assembled when the quadrature is built and returned by reference.

Structural identities

On a periodic mesh derivative_matrix is already antisymmetric, with no skew-symmetrisation needed, because integration by parts over a torus leaves no boundary term:

S = derivative_matrix(q)
maximum(abs, S + S')

6.661338147750939e-16

The stiffness matrix is symmetric positive semi-definite with the constants in its kernel:

K = stiffness_matrix(q)
(maximum(abs, K - K'),
    minimum(eigvals(Symmetric(Matrix(K)))),
    maximum(abs, K * ones(nbasis(q))))

(1.1102230246251565e-16, 6.657001339060997e-17, 8.881784197001252e-16)

And integration by parts holds between the mixed matrices, which is what lets a third-order operator be carried by a $\mathcal{C}^2$ cubic basis:

maximum(abs, mixed_matrix(q, 1, 2) + mixed_matrix(q, 0, 3))

2.955968803064479e-14

A variable coefficient costs no more than a constant one

weighted_matrix takes either a function of the coordinate or a vector already sampled at quadrature_nodes — the second form is what a coefficient that is itself a spline expansion becomes, and it avoids resampling a field already in hand.

A1 = weighted_matrix(q, sin, 0, 1)
A2 = weighted_matrix(q, sin.(quadrature_nodes(q)), 0, 1)
A1 == A2, maximum(abs, weighted_matrix(q, one, 0, 0) - mass_matrix(q))

(true, 1.3877787807814457e-17)

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

eltype(weighted_matrix(q, cis, 0, 0)), eltype(mixed_matrix(q, 0, 0))

(ComplexF64, Float64)

Projection

l2_projection(q, f)            # f a function, or a vector sampled at quadrature_nodes(q)
l2_projection!(û, q, f)        # in place
l2_projection!(s::Spline, q, f)

$\hat{u} = \mathbb{M}^{-1} \int_\Omega f \phi_i$, which is the best approximation of f in the space in the $L^2$ norm.

û = l2_projection(q, sin)
maximum(abs(evaluate(b, û, x) - sin(x)) for x in range(0, 2π; length = 401))

3.4258634479700234e-5

v̂ = similar(û)
l2_projection!(v̂, q, sin)
v̂ == û

true

The in-place form allocates nothing on every representation but one: the $f \odot w$ product goes into a buffer the quadrature holds, 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 — the FactorizedMass case — where CHOLMOD has no in-place solve.

fq = sin.(quadrature_nodes(q))
l2_projection!(v̂, q, fq)                       # warm up
@allocated l2_projection!(v̂, q, fq)

0

A sample whose element type is wider than the quadrature's — a complex f — gets its own product rather than being narrowed into that buffer, so the in-place form accepts exactly what the allocating form accepts and merely stops being allocation-free there.

Mass operators

Solves against the mass matrix go through a MassOperator, and which representation is built is decided by the basis, not by the mesh. A uniform mesh is necessary for circulance but not sufficient: a clamped basis on one has $p$ boundary functions at each end that are not translates of anything.

for (name, bb) in (("periodic, uniform", BSplineBasis(UniformMesh(16, 0 .. 2π), 3, Periodic())),
    ("periodic, graded", BSplineBasis(GradedMesh(16, 0 .. 2π), 3, Periodic())),
    ("periodic, random", BSplineBasis(RandomMesh(16, 0 .. 2π), 3, Periodic())),
    ("clamped, uniform", BSplineBasis(UniformMesh(16, 0 .. 2π), 3)),
    ("Dirichlet, random", BSplineBasis(RandomMesh(16, 0 .. 2π), 3, Dirichlet())))
    println(rpad(name, 20), nameof(typeof(mass_operator(SplineQuadrature(bb)))))
end

periodic, uniform   CirculantMass
periodic, graded    FactorizedMass
periodic, random    FactorizedMass
clamped, uniform    BandedMass
Dirichlet, random   BandedMass

The circulant construction verifies circulance rather than assuming it, to rtol = sqrt(eps(T)) by default — 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. A matrix that is banded but not circulant would still produce plausible numbers through the transform, and the failure would surface much later as a wrong conservation law.

qg = SplineQuadrature(BSplineBasis(GradedMesh(16, 0 .. 2π), 3, Periodic()))
try
    CirculantMass(mass_matrix(qg), nbasis(qg))
catch err
    println(err.msg)
end

the mass matrix is not circulant to a relative tolerance of 1.4901161193847656e-8; a CirculantMass is only valid on a uniform mesh

A singular assembly

The mass matrix is positive definite, but the same three representations carry the other assemblies of a quadrature, and stiffness_matrix on a periodic basis is singular: $-\phi'' = \rho$ says nothing about the constants, and a periodic basis represents them.

bp = BSplineBasis(UniformMesh(64, 0 .. 2π), 3, Periodic())
qp = SplineQuadrature(bp)
Sp = stiffness_matrix(qp)
N = nbasis(bp)
maximum(abs, Sp * ones(N))

1.0658141036401503e-14

By default mass_operator refuses it, since a singular assembly is usually too coarse a quadrature rather than an intended one:

try
    mass_operator(Sp, bp)
catch err
    println(err.msg)
end

the circulant mass matrix has a zero eigenvalue and is not invertible; pass kernel = :project to solve in the complement of the kernel instead

kernel = :project asserts that the kernel is the constants, and asks for the solution that has no constant component. Each representation deflates that subspace in the way its own structure allows — CirculantMass gives the constant mode a zero factor, which is the Moore–Penrose pseudoinverse, and FactorizedMass factorises the minor that drops the first degree of freedom and takes the mean out afterwards. Both verify the assertion, and neither changes the matrix:

op = mass_operator(Sp, bp; kernel = :project)
rhs = sin.(2π .* (0:(N - 1)) ./ N)     # one whole period, so its mean is already zero
φ = op \ rhs
maximum(abs, Sp * φ - rhs), sum(φ)     # it solves the system, in the mean-free gauge

(5.928590951498336e-14, 8.43769498715119e-15)

The textbook cure is the rank-one shift by the mean projector, $\mathbb{S} + \mathbb{1}\mathbb{1}^T/N$, which is invertible and gives the same answer. It is also structurally full, so it costs $O(N^2)$ to store and to factorise where the assembly itself costs $O(N)$:

shifted = mass_operator(Matrix(Sp) .+ inv(N), bp)
maximum(abs, φ - shifted \ rhs),
Base.summarysize(Sp), Base.summarysize(Matrix(Sp) .+ inv(N))

(8.852404527578515e-18, 7848, 32816)

Solving

mass_solve!(y, op, x)          # in place; y and x may alias
ldiv!(y, op, x)   ldiv!(op, x)
op \ x                         # x a vector or a matrix of right-hand sides
mass_matrix(op)   Matrix(op)   size(op)   eltype(op)

op = mass_operator(q)
rhs = basis_integrals(q)
x = op \ rhs
maximum(abs, x .- 1)             # M⁻¹ ∫φᵢ = 1, by the partition of unity

3.9968028886505635e-15

y = similar(rhs)
mass_solve!(y, op, rhs)
y == x, (@allocated mass_solve!(y, op, rhs))

(true, 0)

mass_factorization(q) is the same operator, named for the case where the factorisation rather than the matrix is what is wanted.

Choosing nq

An nq-point Gauß-Legendre rule is exact to degree $2n_q - 1$. Three thresholds matter:

Below the second threshold the antisymmetry does not degrade gracefully — it fails outright, and only on a non-uniform mesh, because a uniform periodic mesh gives circulant assemblies for which it holds at any nq. That is exactly the case where a check run on a uniform mesh confirms the property for the wrong reason:

p = 3
for nq in 2:6
    qu = SplineQuadrature(BSplineBasis(UniformMesh(16, 0 .. 2π), p, Periodic()); nq = nq)
    qr = SplineQuadrature(BSplineBasis(RandomMesh(16, 0 .. 2π), p, Periodic()); nq = nq)
    du = maximum(abs, derivative_matrix(qu) + derivative_matrix(qu)')
    dr = maximum(abs, derivative_matrix(qr) + derivative_matrix(qr)')
    println("nq = ", nq, "   uniform ", rpad(du, 24), "   random ", dr)
end

nq = 2   uniform 2.4147350785597155e-15     random 0.009260061147332663
nq = 3   uniform 8.326672684688674e-16      random 1.1032841307212493e-15
nq = 4   uniform 4.996003610813204e-16      random 1.5221485339659564e-15
nq = 5   uniform 6.661338147750939e-16      random 6.106226635438361e-16
nq = 6   uniform 4.440892098500626e-16      random 6.308373111082252e-16

Too coarse an nq makes the mass matrix singular rather than merely inexact, and the constructor says so by name:

try
    SplineQuadrature(BSplineBasis(UniformMesh(8, 0 .. 1), 3); nq = 1)
catch err
    println(err.msg)
end

the mass matrix assembled with nq = 1 points per cell is not usable for a degree-3 basis on 8 cells: the mass matrix is not positive definite; the quadrature is too coarse to resolve the basis. Use nq ≥ 4 for an exact mass matrix, or at least nq ≥ 2

Two traps