Boundary Conditions

The seven conditions, how to write them, and what the recombined basis they produce does. For why recombination works and what it costs, see Boundary Conditions in the Theory section.

using SimpleSplines
using LinearAlgebra
using SparseArrays

The seven types

type	condition	order $m$	constructor
`Free`	none — the plain clamped basis	$-1$	`Free()`
`Periodic`	the periodic closure of the whole axis	$-1$	`Periodic()`
`Dirichlet`	$u = 0$	$0$	`Dirichlet()`
`Neumann`	$u' = 0$	$1$	`Neumann()`
`Natural`	$u'' = 0$	$2$	`Natural()`
`Robin`	$\alpha u + \beta u' = 0$	$1$, or $0$ if $\beta = 0$	`Robin(α, β)`
`Constraint`	$\sum_k c_{k+1} D^k u = 0$	`findlast(!iszero, c) - 1`	`Constraint(c...)`

Robin promotes its arguments, so Robin(1, 2.0) works, and rejects $\alpha = \beta = 0$. Constraint takes the coefficient of $u$ first, then of $u'$, and so on; it rejects an all-zero and an empty coefficient list. The named types are preferable where they apply: they say what is meant, and Dirichlet in addition takes the cheaper elimination path.

Robin(1, 2.0), Constraint(0, 0, 0, 1), constraint_order(Constraint(1, 0))

(Robin(1.0, 2.0), Constraint(0, 0, 0, 1), 0)

Three traits, one meaning

Every condition is given its meaning in exactly one place, so that the recombination has a single implementation:

for bc in (Free(), Periodic(), Dirichlet(), Neumann(), Natural(),
    Robin(2.0, 3.0), Robin(1.0, 0.0), Constraint(1, 0, -2))
    println(rpad(repr(bc), 20),
        "  coefficients = ", rpad(repr(constraint_coefficients(bc)), 12),
        "  order = ", rpad(constraint_order(bc), 3),
        "  costs = ", nconstraints(bc))
end

Free()                coefficients = nothing       order = -1   costs = 0
Periodic()            coefficients = nothing       order = -1   costs = 0
Dirichlet()           coefficients = (1,)          order = 0    costs = 1
Neumann()             coefficients = (0, 1)        order = 1    costs = 1
Natural()             coefficients = (0, 0, 1)     order = 2    costs = 1
Robin(2.0, 3.0)       coefficients = (2.0, 3.0)    order = 1    costs = 1
Robin(1.0, 0.0)       coefficients = (1.0, 0.0)    order = 0    costs = 1
Constraint(1, 0, -2)  coefficients = (1, 0, -2)    order = 2    costs = 1

constraint_coefficients — the $(c_0, c_1, \dots)$, or nothing for the two conditions that impose none.
constraint_order — the highest derivative appearing, or -1. Note constraint_order(Robin(1.0, 0.0)) == 0: the trailing zero does not count.
nconstraints — the degrees of freedom removed: 0 or 1, never more, whatever the order.

Writing a specification

boundary_conditions normalises whatever you pass:

you write	you get
a single condition	`(bc, bc)` — both ends
a two-tuple	`(left, right)`
`Periodic()`	`Periodic()`, bare — not a pair
a `Symbol`, or a tuple of them	the same, with each symbol resolved

(boundary_conditions(Dirichlet()),
    boundary_conditions((Dirichlet(), Neumann())),
    boundary_conditions(:periodic),
    boundary_conditions((:free, Natural())))

((Dirichlet(), Dirichlet()), (Dirichlet(), Neumann()), Periodic(), (Free(), Natural()))

Periodic comes back unpaired because it is a condition on the axis, not on its ends: there is no such thing as being periodic at the left end only, so pairing it with anything is rejected.

try
    boundary_conditions((Periodic(), Dirichlet()))
catch err
    println(err.msg)
end

Periodic identifies the two ends of the domain, so it cannot be given for one end only; pass `Periodic()` for the whole axis

Symbol sugar

:periodic, :dirichlet, :neumann, :natural, :free. Lowercase only, and deliberately not case-insensitive: :Natural, :nothing, :Dirichlet and :Periodic are rejected with a message saying what to write instead, and any other symbol with the list of the five that are accepted.

for s in (:Natural, :nothing, :Dirichlet, :quasiperiodic)
    try
        BoundaryCondition(s)
    catch err
        println(err.msg, "\n")
    end
end

boundary condition :Natural is not accepted: `:Natural` used to mean the unconstrained clamped basis, which is now `Free()`. `:natural` now means the natural condition u'' = 0, i.e. `Natural()`. Say which you mean.

boundary condition :nothing is not accepted: `:nothing` used to select the unconstrained clamped basis by falling through a branch. Write `Free()`.

boundary condition :Dirichlet is not accepted: boundary conditions are given in lower case or as types: `:dirichlet` or `Dirichlet()`.

unknown boundary condition :quasiperiodic; the accepted symbols are :dirichlet, :free, :natural, :neumann, :periodic, or pass a BoundaryCondition directly

Which basis a condition selects

m = UniformMesh(8, 0 .. 1)
for bc in (Free(), Periodic(), Dirichlet(), Neumann(), Natural(),
    Robin(1.0, 2.0), (Dirichlet(), Free()), (Natural(), Neumann()))
    bb = BSplineBasis(m, 3, bc)
    println(rpad(repr(bc), 24), rpad(nameof(typeof(bb)), 24),
        "N = ", rpad(nbasis(bb), 4), "reproduces ", polynomial_reproduction(bb))
end

Free()                  BSplineBasis            N = 11  reproduces 3
Periodic()              PeriodicBSplineBasis    N = 8   reproduces 0
Dirichlet()             RecombinedBSplineBasis  N = 9   reproduces -1
Neumann()               RecombinedBSplineBasis  N = 9   reproduces 0
Natural()               RecombinedBSplineBasis  N = 9   reproduces 1
Robin(1.0, 2.0)         RecombinedBSplineBasis  N = 9   reproduces -1
(Dirichlet(), Free())   RecombinedBSplineBasis  N = 10  reproduces -1
(Natural(), Neumann())  RecombinedBSplineBasis  N = 9   reproduces 0

Read it back with boundary:

(boundary(BSplineBasis(m, 3)),
    boundary(BSplineBasis(m, 3, Periodic())),
    boundary(BSplineBasis(m, 3, (Dirichlet(), Neumann()))))

((Free(), Free()), Periodic(), (Dirichlet(), Neumann()))

Working with a recombined basis

A RecombinedBSplineBasis answers the whole AbstractBSplineBasis interface, so most code needs to know nothing about it. Four things are specific to it.

Its parent is reachable. parent(b) is the clamped basis the recombination was built from, and recombination_matrix is the sparse $R$ with $\psi_j = \sum_i R_{ij}\varphi_i$.

b = BSplineBasis(m, 3, Natural())
bp = parent(b)
R = recombination_matrix(b)
nbasis(bp), nbasis(b), size(R), nnz(R)

(11, 9, (11, 9), 13)

Every assembly is the parent's conjugated by $R$, which is worth knowing when comparing against a hand-built reference:

q, qp = SplineQuadrature(b), SplineQuadrature(bp)
maximum(abs, Matrix(mass_matrix(q)) - Matrix(R' * mass_matrix(qp) * R))

1.1102230246251565e-16

local_width may exceed $p+1$. A recombined function near an end spans the union of two parent supports, so the local block is wider there. Any buffer for evaluate_all! must have exactly that many entries, and how many depends on the condition rather than on the degree alone:

[(bc, degree(BSplineBasis(m, 3, bc)) + 1, local_width(BSplineBasis(m, 3, bc)))
 for bc in (Free(), Dirichlet(), Neumann(), Natural(), Constraint(0, 0, 0, 1))]

5-element Vector{Tuple{BoundaryCondition, Int64, Int64}}:
 (Free(), 4, 4)
 (Dirichlet(), 4, 4)
 (Neumann(), 4, 4)
 (Natural(), 4, 5)
 (Constraint(0, 0, 0, 1), 4, 6)

degree(b) + 1, local_width(b),
[length(collect(local_indices(b, k))) for k in 1:ncells(b)]

(4, 5, [3, 4, 5, 4, 4, 5, 4, 3])

try
    evaluate_all!(zeros(degree(b) + 1), b, 0.3)
catch err
    println(err.msg)
end

evaluate_all! needs a buffer of 5 entries for this basis, got 4; see `local_width`

So local_width(b), never degree(b) + 1, is what a buffer should be sized by.

nodes is not a Schoenberg-Whitney set. For a recombined basis nodes returns the Greville abscissa of the parent function each column carries with unit coefficient. Those are N distinct increasing points inside the domain, interlaced with the supports — which is what a plotting or collocation grid wants — but the collocation matrix $\psi_j(\xi_i)$ is not guaranteed invertible, and that is not claimed.

nodes(b)

9-element Vector{Float64}:
 0.0
 0.041666666666666664
 0.25
 0.375
 0.5
 0.625
 0.75
 0.9583333333333334
 1.0

The dimension is n + p less one per constrained end, whatever the order of the condition. Natural() reaches three parent functions and removes one, exactly as Dirichlet() reaches one and removes one.

[(bc, nbasis(BSplineBasis(m, 4, bc)))
 for bc in (Dirichlet(), Neumann(), Natural(), Constraint(0, 0, 0, 1))]

4-element Vector{Tuple{BoundaryCondition, Int64}}:
 (Dirichlet(), 10)
 (Neumann(), 10)
 (Natural(), 10)
 (Constraint(0, 0, 0, 1), 10)

Traps

A condition is imposed on the space, not on a solution. There is no "apply the boundary condition" step after assembling: the basis has one function fewer per constrained end, so a stiffness matrix built on it is already the matrix of the constrained problem. Do not also delete rows.
Homogeneous only. $u(a) = g \ne 0$ is not a boundary condition here. Write $u = u_0 + w$ with $u_0$ any function taking the required values and $w$ in the homogeneous space, and solve for $w$.
Free is not Natural. Free() is the absence of a condition; the natural boundary condition is $u'' = 0$, which is Natural. They span different spaces and have different dimensions. The confusion is easy because a clamped basis is sometimes loosely called a "natural" spline basis, and it is why :Natural is rejected as a symbol.
A Dirichlet basis does not reproduce the constants, so $\int_\Omega u_h \, \mathrm{d}x$ is not conserved by anything built on it and basis_integrals!= mass_matrix(q) * ones(N). Check polynomial_reproduction before assuming a conservation law survives.
An ill-scaled Robin degrades silently. Robin(1.0, 1e-20) makes the recombination anchor numerically zero; the mass matrix comes out finite with condition number already Inf, and the Cholesky reports success. Keep the two coefficients within a few orders of magnitude of each other, and use the named types where they apply.

bd = BSplineBasis(m, 3, Dirichlet())
qd = SplineQuadrature(bd)
maximum(abs, basis_integrals(qd) - mass_matrix(qd) * ones(nbasis(bd)))

0.010937499999999989

That is the size of the discrepancy the partition of unity would have made zero.