SimpleSplines.jl

B-spline finite elements on an interval, built from the Cox-de Boor recursion, with the quadrature and assembly a Galerkin discretisation needs.

The package is deliberately small: one family of bases, one assembly table, and the mass, stiffness, derivative and variable-coefficient matrices as single weighted contractions of that table. What it is not is a curve- and surface-modelling library — there are no NURBS, no knot insertion, no degree elevation and no least-squares fitting of data. What it is for is discretising a differential equation and then solving with the result.

What it provides


Bases	clamped (`BSplineBasis`), periodic (`PeriodicBSplineBasis`), and recombined (`RecombinedBSplineBasis`) — of any degree $p \ge 0$
Meshes	`UniformMesh`, `GradedMesh`, `RandomMesh`, `GeneralMesh`
Boundary conditions	`Free`, `Dirichlet`, `Neumann`, `Natural`, `Robin`, and the general `Constraint` — per end, plus `Periodic` for the whole basis
Assembly	`SplineQuadrature` and, from it, `mass_matrix`, `stiffness_matrix`, `derivative_matrix`, `mixed_matrix`, `weighted_matrix`, `basis_integrals`
Mass solves	`CirculantMass` (FFT), `BandedMass` (banded Cholesky), `FactorizedMass` (sparse Cholesky), `KroneckerMass` (factored, $D$-dimensional)
Tensor products	`TensorProductBasis` and `TensorProductQuadrature` in any number of dimensions, with degree, mesh, domain and boundary condition per axis
Functions	`Spline`, `SplineDerivative`, `l2_projection`
Particle path	`findcell` and `evaluate_all!`, allocation-free, for depositing onto the basis

Installation

using Pkg
Pkg.add(url = "https://github.com/JuliaDEC/SimpleSplines.jl")

In thirty seconds

A basis is a Mesh, a degree, and a BoundaryCondition:

using SimpleSplines

b = BSplineBasis(UniformMesh(16, 0 .. 1), 3, Dirichlet())
nbasis(b), degree(b), order(b), polynomial_reproduction(b)

(17, 3, 4, -1)

Assembly goes through a SplineQuadrature, which tabulates the basis and its derivatives at the global Gauß-Legendre points and hands back every matrix built on them:

q = SplineQuadrature(b)
M = mass_matrix(q)
K = stiffness_matrix(q)
size(M), size(K)

((17, 17), (17, 17))

Fitting a function to the space is an $L^2$ projection, and the result is callable:

u = Spline(b, l2_projection(q, x -> sin(π * x)))
u(0.5), u(0.5, 1), u(0.0)        # value, first derivative, and the imposed u(0) = 0

(1.000002084343635, 2.6645352591003757e-15, 0.0)

Solving $-u'' = f$ with $u(0) = u(1) = 0$ is then two lines, because the recombined basis has already taken care of the boundary condition:

f(x) = π^2 * sin(π * x)
rhs = basis_values(q, 0) * (quadrature_weights(q) .* f.(quadrature_nodes(q)))
û = Matrix(K) \ rhs
maximum(abs(evaluate(b, û, x) - sin(π * x)) for x in range(0, 1; length = 101))

2.0817761019475967e-6

How the pieces fit

Mesh + BoundaryCondition
        │
        ▼
AbstractBSplineBasis  ───▶  Spline        ───▶  s(x),  derivative(s)
        │
        ▼
SplineQuadrature      ───▶  mass_matrix, stiffness_matrix, derivative_matrix,
        │                   mixed_matrix, weighted_matrix, basis_integrals,
        │                   l2_projection
        ▼
MassOperator          ───▶  mass_solve!,  op \ x

The mesh is geometry alone and carries no boundary condition; the boundary condition decides which of the three bases the one constructor BSplineBasis(mesh, p, bc) returns; the basis decides which representation the mass matrix takes. A tensor product is the same picture with a tuple in each box.

Where to go next

Tutorial — one worked example from a bare using to a solved problem.
Theory — B-Splines for the recursion, the derivative formula and the dimension counts; Boundary Conditions for recombination; Tensor Products for the Kronecker structure.
Usage — Meshes, Bases, Boundary Conditions, Assembly, Tensor Products and Splines give the constructors, the accessors and the traps.
Gallery — nine solved problems with their measured errors.
Library — every exported name.