diff --git a/CHANGELOG.md b/CHANGELOG.md index 4d83635..efc0066 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,55 @@ The first release with any implementation in it. Before this the package was a `PkgTemplates` skeleton: a registered UUID, CI, docs scaffolding and a declared dependency set, with `src/SimpleSplines.jl` a bare `module … end`. +### The manual + +The documentation was one page ending in an `@autodocs` dump, and its opening two paragraphs +still described the package as periodic-only with "one basis" — which had not been true since +the boundary conditions and tensor products landed. It is now a manual of fourteen pages: + +| | | +|:--|:--| +| **Tutorial** | one narrative from a bare `using` to a solved boundary-value problem | +| **Theory** | the Cox-de Boor recursion, the derivative recursion, support, smoothness and the three dimension counts; recombination and what it costs; the tensor-product space and its flattening convention | +| **Usage** | one page per object — meshes, bases, boundary conditions, assembly, tensor products, splines — with the exact constructors, the accessor tables, what is rejected, and the traps | +| **Gallery** | eight solved problems, each printing what it measured — a convergence rate where one is meaningful, otherwise the errors themselves | +| **Library** | every export, grouped by topic, with explicit `@docs` blocks rather than `@autodocs` | +| **References** | a bibliography, via `DocumenterCitations` | + +Everything computational in it is an `@example` block, so the numbers, the rates and the +sixteen `CairoMakie` figures are produced by the build and a claim that stops being true stops +the build. Three of the gallery problems are new rather than lifted from the suite: Poisson in +one dimension with a `Dirichlet` basis and its variable-coefficient variant, which reports a +measured order of `p+1`; the Dirichlet eigenvalues against `(kπ)²`, which reports the relative +error per mode and per degree; and Poisson on a box through the separable Laplacian, which +reports the error at three refinements. + +The library page is written out by hand instead of generated, for two reasons: `@autodocs` +rendered the internal `_bspline` alongside the public API, and an explicit list is what makes +a missing docstring visible. `checkdocs = :exports` is on. + +**One documentation bug fixed.** The `@doc raw` block whose signatures read +`findcell(b, x)` / `findcell(m, x)` sat immediately above the internal `_findcell` and so +documented *that*: `?findcell` showed nothing, `@docs findcell` would have failed the build, +and `@autodocs` rendered the text under `SimpleSplines._findcell`. It is now attached to +`findcell` itself. The same class of mistake as a comment between a docstring and its +definition, and equally silent — nothing but a docs build notices. + +**One complexity claim corrected.** Writing the cost table for the manual turned up a wrong +figure in the source: `evaluate`'s and `evaluate_all!`'s docstrings gave the single-function +path as `O(p²)`, but `_bspline` runs the recursion unmemoised — it splits into two subproblems +at every level — so one value costs `O(2^p)`. Measured over 20 000 evaluations on a 40-cell +mesh, `evaluate` takes 26 ns at `p = 3` and 15.5 µs at `p = 12`, against 21 ns and 104 ns for +the *whole block* through `evaluate_all`. Nothing about the implementation changed; the +choice to run the formula as written is deliberate, and the reason `evaluate_all!` exists. +What was wrong was the number, in the one place a reader would take it from. + +Docstrings were added for the exports that had none: `quadratures`, +`PeriodicBSplineDerivative`, `nnodes`, `grid`, `basis` on a `SplineQuadrature`, and `basis`, +`mass_operator`, `mass_matrix` and `nnodes` on the tensor-product types. `Basis`, `..`, +`leftendpoint` and `rightendpoint` are re-exported unchanged and are documented by the packages +that own them, which the library page says. + ### General boundary conditions and tensor products The package was periodic-only: one basis, on `[0,L)`, on a torus. It now covers bounded @@ -491,9 +540,9 @@ Found in review of the branch, not by the suite, and each now has a regression t reaches, because only that one contributes to `a_{m+1} = c_m D^m φ_{m+1}(a)`. That anchor is nonzero whenever the leading coefficient `c_m` is — but *how* nonzero is the caller's, not the construction's. `Robin(1.0, 1e-20)` puts a small coefficient on the highest derivative, the -anchor is then numerically zero, and the recombination coefficients `-a_t/a_{m+1}` overflow: +anchor is then numerically zero, and the recombination coefficients `-a_t/a_{m+1}` become huge: the mass matrix comes out finite with a condition number already `Inf`, `cholesky(…; -check = false)` reports `issuccess` on a matrix containing `Inf`, and the projection that +check = false)` reports `issuccess` on it anyway, and the projection that follows looks plausible. A finiteness guard would catch only the most extreme case and is a symptom patch — the diff --git a/README.md b/README.md index 8072060..dd6358b 100644 --- a/README.md +++ b/README.md @@ -7,11 +7,20 @@ [![Coverage](https://codecov.io/gh/JuliaDEC/SimpleSplines.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDEC/SimpleSplines.jl) [![PkgEval](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/S/SimpleSplines.svg)](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/S/SimpleSplines.html) -Periodic 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 provides the -periodic B-spline basis of arbitrary degree on uniform, graded and random meshes, together -with an assembly table from which the mass, stiffness, derivative and variable-coefficient -matrices follow as single weighted contractions. +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 provides the clamped, periodic and recombined B-spline bases of arbitrary degree +on uniform, graded, random or arbitrary meshes; homogeneous Dirichlet, Neumann, Robin, natural +and general local boundary conditions per end; tensor products in any number of dimensions, +with degree, mesh, domain and boundary condition per axis; and an assembly table from which the +mass, stiffness, derivative and variable-coefficient matrices follow as single weighted +contractions. Mass solves go through a representation chosen by the basis — an FFT for a +periodic uniform basis, a banded Cholesky for a bounded one, and a factored Kronecker product +in several dimensions. + +The [manual](https://JuliaDEC.github.io/SimpleSplines.jl/dev/) has a tutorial, the spline +theory the package rests on, a gallery of solved problems and the full API. ## Development diff --git a/docs/Project.toml b/docs/Project.toml index 7b31ff2..c7b1cb9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,10 +1,16 @@ [deps] -CompactBasisFunctions = "a09551c4-f815-4143-809e-dd1986753ba7" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SimpleSplines = "ac480f6c-9a94-4d53-a283-9db69f109154" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [sources] SimpleSplines = {path = ".."} [compat] +CairoMakie = "0.12, 0.13, 0.14, 0.15" Documenter = "1" +DocumenterCitations = "1" diff --git a/docs/make.jl b/docs/make.jl index c6ef518..0359dd7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,19 +1,46 @@ using SimpleSplines using Documenter +using DocumenterCitations DocMeta.setdocmeta!(SimpleSplines, :DocTestSetup, :(using SimpleSplines); recursive = true) +bib = CitationBibliography(joinpath(@__DIR__, "src", "references.bib")) + makedocs(; + plugins = [bib], modules = [SimpleSplines], authors = "Michael Kraus", sitename = "SimpleSplines.jl", + checkdocs = :exports, format = Documenter.HTML(; + prettyurls = get(ENV, "CI", "false") == "true", canonical = "https://JuliaDEC.github.io/SimpleSplines.jl", edit_link = "main", + # the library page collects every docstring in the package and is legitimately + # large; the default warning threshold is aimed at pages that are large by accident + size_threshold_warn = 400 * 1024, + size_threshold = 800 * 1024, assets = String[] ), pages = [ "Home" => "index.md", + "Tutorial" => "tutorial.md", + "Theory" => [ + "B-Splines" => "theory/bsplines.md", + "Boundary Conditions" => "theory/boundary.md", + "Tensor Products" => "theory/tensorproduct.md" + ], + "Usage" => [ + "Meshes" => "usage/meshes.md", + "Bases" => "usage/bases.md", + "Boundary Conditions" => "usage/boundary.md", + "Assembly" => "usage/quadrature.md", + "Tensor Products" => "usage/tensorproduct.md", + "Splines" => "usage/splines.md" + ], + "Gallery" => "gallery.md", + "Library" => "library.md", + "References" => "references.md" ] ) diff --git a/docs/src/gallery.md b/docs/src/gallery.md new file mode 100644 index 0000000..e05de23 --- /dev/null +++ b/docs/src/gallery.md @@ -0,0 +1,467 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# Gallery + +Eight problems, solved and measured. Every number and figure on this page is produced when the +manual is built, so nothing here is a claim about what the package should do. + +```@example gal +using SimpleSplines +using CairoMakie +using LinearAlgebra +CairoMakie.activate!(type = "png") # hide +nothing # hide +``` + +## 1. ``L^2`` projection converges at order ``p+1`` + +The projection of a smooth function onto a degree-``p`` spline space is accurate to +``O(h^{p+1})``. The mesh here is a [`GradedMesh`](@ref) rather than a uniform one, deliberately: +refining a graded mesh gives a genuine *family*, so a rate measured across it means something, +whereas a uniform mesh would make every assembly circulant and could confirm the rate for the +wrong reason. + +```@example gal +target(x) = sin(x) +sample = range(0, 2π; length = 400)[1:399] + +function projection_error(p, n) + b = BSplineBasis(GradedMesh(n, 2π), p, Periodic()) + q = SplineQuadrature(b; nq = quadrature_order(p) + 3) + û = l2_projection(q, target) + (meshwidth(b), maximum(abs(evaluate(b, û, x) - target(x)) for x in sample)) +end + +ns = (8, 16, 32, 64) +for p in 1:4 + hs = Float64[] + es = Float64[] + for n in ns + h, e = projection_error(p, n) + push!(hs, h) + push!(es, e) + end + rate = log(es[end - 1] / es[end]) / log(hs[end - 1] / hs[end]) + println("p = ", p, " error at n = 64: ", es[end], + " measured rate: ", round(rate; digits = 3), " (expected ", p + 1, ")") +end +``` + +```@example gal +fig = Figure(size = (620, 420)) +ax = Axis(fig[1, 1]; xscale = log10, yscale = log10, xlabel = "h", ylabel = "max error", + title = "L² projection of sin on a graded mesh") +for p in 1:4 + data = [projection_error(p, n) for n in ns] + hh, ee = first.(data), last.(data) + scatterlines!(ax, hh, ee; label = "p = $(p)") + lines!(ax, hh, ee[1] .* (hh ./ hh[1]) .^ (p + 1); + color = (:black, 0.35), linestyle = :dash) +end +axislegend(ax; position = :rb) +fig +``` + +The dashed lines are ``h^{p+1}`` for reference. + +## 2. Which polynomials survive which boundary condition + +An ``L^2`` projection reproduces a polynomial exactly when that polynomial lies in the span of +the basis, and [`polynomial_reproduction`](@ref) says which do. This is what decides whether a +conservation law survives a discretisation, so it is worth seeing rather than trusting. + +```@example gal +function reproduces(bc, p, k) + b = BSplineBasis(UniformMesh(8, 0 .. 1), p, bc) + q = SplineQuadrature(b) + û = l2_projection(q, x -> x^k) + maximum(abs(evaluate(b, û, x) - x^k) for x in range(0, 1; length = 201)) +end + +conditions = (Free(), Periodic(), Dirichlet(), Neumann(), Natural(), + Constraint(0, 0, 0, 1), Robin(1.0, 1.0)) +println(rpad("condition", 24), rpad("reproduces", 12), + "error on 1, x, x², x³, x⁴ (p = 4)") +for bc in conditions + r = polynomial_reproduction(BSplineBasis(UniformMesh(8, 0 .. 1), 4, bc)) + errs = [reproduces(bc, 4, k) for k in 0:4] + println(rpad(repr(bc), 24), rpad(r, 12), + join((rpad(e < 1e-9 ? "exact" : string(round(e; sigdigits = 2)), 10) + for e in errs))) +end +``` + +Read the table along each row: the monomials up to degree `r` come back exactly, and the one +after does not. `Free` reproduces everything up to its degree, `Periodic` only the constants +(``x`` is not periodic), `Dirichlet` nothing at all, and `Constraint(0, 0, 0, 1)` — the +condition ``u''' = 0`` — reproduces up to ``x^2``. + +## 3. Structural identities + +Three identities that hold exactly, and one that holds only on a uniform mesh. The point of the +[`RandomMesh`](@ref) column is that a check run only on a [`UniformMesh`](@ref) confirms a +property whose proof does not need the mesh at all — the assemblies there are circulant. + +```@example gal +meshes = ("uniform" => UniformMesh(16, 2π), + "graded" => GradedMesh(16, 2π), + "random" => RandomMesh(16, 2π)) + +println(rpad("", 10), rpad("S + Sᵀ", 12), rpad("∫φ'φ'' + ∫φφ'''", 18), + rpad("K·1", 12), "λ_min(K)") +for (name, m) in meshes + q = SplineQuadrature(BSplineBasis(m, 3, Periodic())) + S = derivative_matrix(q) + K = stiffness_matrix(q) + println(rpad(name, 10), + rpad(round(maximum(abs, S + S'); sigdigits = 2), 12), + rpad(round(maximum(abs, mixed_matrix(q, 1, 2) + mixed_matrix(q, 0, 3)); + sigdigits = 2), 18), + rpad(round(maximum(abs, K * ones(nbasis(q))); sigdigits = 2), 12), + round(minimum(eigvals(Symmetric(Matrix(K)))); sigdigits = 2)) +end +``` + + - ``S + S^{\mathsf T} = 0``: integration by parts over a torus leaves no boundary term. + - ``\int \phi'\phi'' + \int \phi\phi''' = 0``: the same, one derivative higher. This is what + lets a ``\mathcal{C}^2`` cubic basis carry a third-order operator. + - ``\mathbb{K}\mathbf{1} = 0`` and ``\lambda_{\min}(\mathbb{K}) = 0``: the constants are in + the kernel of the stiffness matrix, because they are in the span of the basis. + +And the one that is *not* mesh-independent — the antisymmetry needs the quadrature to +integrate a total derivative of degree ``2p-1`` exactly, i.e. ``n_q \ge p``: + +```@example gal +for nq in 2:4 + qu = SplineQuadrature(BSplineBasis(UniformMesh(16, 2π), 3, Periodic()); nq = nq) + qr = SplineQuadrature(BSplineBasis(RandomMesh(16, 2π), 3, Periodic()); nq = nq) + println("nq = ", nq, + " uniform ", rpad(round(maximum(abs, + derivative_matrix(qu) + derivative_matrix(qu)'); sigdigits = 2), 12), + " random ", round(maximum(abs, + derivative_matrix(qr) + derivative_matrix(qr)'); sigdigits = 2)) +end +``` + +At `nq = 2` on a random mesh the defect is ``\approx 10^{-2}``, not ``10^{-15}``. On the +uniform mesh it is machine precision at every `nq`, which is why the check has to be run +somewhere else. + +## 4. A non-separable function in two dimensions + +The tensor-product basis is separable; the function projected onto it need not be. + +```@example gal +gfun(x) = sin(x[1]) * exp(-x[2]^2 / 8) + cos(2x[1]) * x[2] / 20 + +grid2 = [(2π * (i - 0.5) / 37, -10 + 20 * (j - 0.5) / 41) for i in 1:37, j in 1:41] +errs2 = Float64[] +for n in (8, 16, 32) + Bn = BSplineBasis(UniformMesh(n, 0 .. 2π), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(n, -10 .. 10), 3) + Qn = TensorProductQuadrature(Bn) + ĝ = l2_projection(Qn, gfun) + push!(errs2, maximum(abs(evaluate(Bn, ĝ, x) - gfun(x)) for x in grid2)) +end +errs2, errs2[1] / errs2[2], errs2[2] / errs2[3] +``` + +Cubics on both axes, so the error should fall by about ``2^4 = 16`` per refinement. + +```@example gal +Bshow = BSplineBasis(UniformMesh(24, 0 .. 2π), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(24, -10 .. 10), 3) +ĝshow = l2_projection(TensorProductQuadrature(Bshow), gfun) +gx = range(0, 2π; length = 121) +gy = range(-10, 10; length = 121) +Zh = [evaluate(Bshow, ĝshow, (x, y)) for x in gx, y in gy] +Ze = [gfun((x, y)) for x in gx, y in gy] + +fig = Figure(size = (780, 320)) +ax1 = Axis(fig[1, 1]; xlabel = "x₁", ylabel = "x₂", title = "projection, n = 24") +hm = heatmap!(ax1, gx, gy, Zh) +Colorbar(fig[1, 2], hm) +ax2 = Axis(fig[1, 3]; xlabel = "x₁", ylabel = "x₂", title = "error") +he = heatmap!(ax2, gx, gy, Zh .- Ze; colormap = :balance) +Colorbar(fig[1, 4], he) +fig +``` + +## 5. Poisson in one dimension + +```math +-u'' = f \ \text{ on } (0,1) , \qquad u(0) = u(1) = 0 , +``` + +discretised by Galerkin on a [`Dirichlet`](@ref)-recombined basis: +``\mathbb{K}\hat{u} = \mathbf{b}`` with ``\mathbb{K} = \int \psi_k'\psi_l'`` the +[`stiffness_matrix`](@ref) and ``b_i = \int f \psi_i``. There is no separate step that imposes +the boundary condition — the basis already has one function fewer per end, so `K` *is* the +matrix of the constrained problem. + +```@example gal +exact(x) = sinpi(x) +rhsfun(x) = π^2 * sinpi(x) + +function poisson1d(p, n) + b = BSplineBasis(UniformMesh(n, 0 .. 1), p, Dirichlet()) + q = SplineQuadrature(b) + K = stiffness_matrix(q) + load = basis_values(q, 0) * (quadrature_weights(q) .* rhsfun.(quadrature_nodes(q))) + û = Matrix(K) \ load + (b, û, maximum(abs(evaluate(b, û, x) - exact(x)) + for x in range(0, 1; length = 401))) +end + +for p in 2:4 + es = [poisson1d(p, n)[3] for n in (8, 16, 32, 64)] + println("p = ", p, " errors ", round.(es; sigdigits = 3), + " rates ", round.([log2(es[i] / es[i + 1]) for i in 1:3]; digits = 2)) +end +``` + +Order ``p+1``, as for the projection. The `Dirichlet` condition holds identically: + +```@example gal +bP, ûP, _ = poisson1d(3, 16) +uh = Spline(bP, ûP) +uh(0.0), uh(1.0), uh(0.5) +``` + +```@example gal +fig = Figure(size = (780, 340)) +xs = range(0, 1; length = 401) +ax1 = Axis(fig[1, 1]; xlabel = "x", title = "-u'' = π² sin(πx), p = 3, n = 16") +lines!(ax1, xs, exact.(xs); label = "exact") +lines!(ax1, xs, uh.(xs); linestyle = :dash, linewidth = 3, label = "Galerkin") +ax2 = Axis(fig[1, 2]; xlabel = "x", title = "error") +lines!(ax2, xs, uh.(xs) .- exact.(xs)) +Legend(fig[2, 1:2], ax1; orientation = :horizontal, framevisible = false) +fig +``` + +### A variable coefficient costs nothing extra + +```math +-\big(a(x) u'\big)' = f , \qquad a(x) = 1 + x^2 , +``` + +is the same problem with [`weighted_matrix`](@ref)`(q, a, 1, 1)` in place of the stiffness +matrix — one weighted contraction of the same tabulation. + +```@example gal +acoef(x) = 1 + x^2 +fvar(x) = π^2 * acoef(x) * sinpi(x) - 2x * π * cospi(x) # so that u = sin(πx) again + +function poisson1d_variable(p, n) + b = BSplineBasis(UniformMesh(n, 0 .. 1), p, Dirichlet()) + q = SplineQuadrature(b) + A = weighted_matrix(q, acoef, 1, 1) + load = basis_values(q, 0) * (quadrature_weights(q) .* fvar.(quadrature_nodes(q))) + û = Matrix(A) \ load + maximum(abs(evaluate(b, û, x) - exact(x)) for x in range(0, 1; length = 401)) +end + +for p in 2:4 + es = [poisson1d_variable(p, n) for n in (8, 16, 32, 64)] + println("p = ", p, " errors ", round.(es; sigdigits = 3), + " rates ", round.([log2(es[i] / es[i + 1]) for i in 1:3]; digits = 2)) +end +``` + +## 6. The Dirichlet eigenvalues + +```math +-u'' = \lambda u \ \text{ on } (0,1) , \qquad u(0) = u(1) = 0 , +\qquad \lambda_k = (k\pi)^2 , +``` + +is the generalised eigenproblem ``\mathbb{K} v = \lambda \mathbb{M} v``. The mass matrix is not +the identity — the basis is not orthonormal — so this is a genuine pencil and not a plain +eigenproblem. + +```@example gal +exact_λ = [(k * π)^2 for k in 1:6] +for p in 2:4 + b = BSplineBasis(UniformMesh(32, 0 .. 1), p, Dirichlet()) + q = SplineQuadrature(b) + λ = eigvals(Symmetric(Matrix(stiffness_matrix(q))), + Symmetric(Matrix(mass_matrix(q)))) + rel = abs.(λ[1:6] .- exact_λ) ./ exact_λ + println("p = ", p, " relative error on λ₁…λ₆: ", round.(rel; sigdigits = 2)) +end +``` + +Every eigenvalue is approximated from **above** — a Galerkin discretisation of a symmetric +positive problem restricts the Rayleigh quotient to a subspace — and the error grows with the +mode number, since the higher modes are the ones the mesh does not resolve. + +```@example gal +b6 = BSplineBasis(UniformMesh(32, 0 .. 1), 3, Dirichlet()) +q6 = SplineQuadrature(b6) +λ6, V6 = eigen(Symmetric(Matrix(stiffness_matrix(q6))), + Symmetric(Matrix(mass_matrix(q6)))) +all(λ6[1:6] .≥ exact_λ) +``` + +```@example gal +fig = Figure(size = (780, 300)) +ax1 = Axis(fig[1, 1]; xlabel = "x", title = "the first three eigenfunctions") +xs = range(0, 1; length = 401) +for k in 1:3 + v = V6[:, k] + v ./= maximum(abs, [evaluate(b6, v, x) for x in xs]) + lines!(ax1, xs, [evaluate(b6, v, x) for x in xs]; label = "k = $(k)") +end +ylims!(ax1, -1.1, 1.55) +axislegend(ax1; position = :rt, orientation = :horizontal, framevisible = false) +ax2 = Axis(fig[1, 2]; yscale = log10, xlabel = "k", ylabel = "relative error", + title = "eigenvalue error, n = 32") +for p in 2:4 + bb = BSplineBasis(UniformMesh(32, 0 .. 1), p, Dirichlet()) + qq = SplineQuadrature(bb) + ll = eigvals(Symmetric(Matrix(stiffness_matrix(qq))), + Symmetric(Matrix(mass_matrix(qq)))) + scatterlines!(ax2, 1:10, abs.(ll[1:10] .- [(k * π)^2 for k in 1:10]) ./ + [(k * π)^2 for k in 1:10]; label = "p = $(p)") +end +axislegend(ax2; position = :rb) +fig +``` + +## 7. Poisson on a box + +```math +-\Delta u = f \ \text{ on } (0,1)^2 , \qquad u = 0 \text{ on } \partial\Omega . +``` + +The Laplacian *is* separable, as a sum of ``D`` separable terms: with +``\mathbb{K}^{(k)}`` and ``\mathbb{M}^{(k)}`` the one-dimensional stiffness and mass matrices, + +```math +\mathbb{A} = \mathbb{M}^{(2)} \otimes \mathbb{K}^{(1)} + + \mathbb{K}^{(2)} \otimes \mathbb{M}^{(1)} , +``` + +in the first-axis-fastest ordering this package uses throughout. The right-hand side is the +[`contract`](@ref)ion of ``f`` sampled on the quadrature grid. + +```@example gal +exact2(x) = sinpi(x[1]) * sinpi(x[2]) +rhs2(x) = 2π^2 * sinpi(x[1]) * sinpi(x[2]) + +function poisson2d(p, n) + B = BSplineBasis(UniformMesh(n, 0 .. 1), p, Dirichlet()) ⊗ + BSplineBasis(UniformMesh(n, 0 .. 1), p, Dirichlet()) + Q = TensorProductQuadrature(B) + q1, q2 = quadratures(Q) + A = kron(mass_matrix(q2), stiffness_matrix(q1)) + + kron(stiffness_matrix(q2), mass_matrix(q1)) + L = contract(Q, quadrature_sample(Q, rhs2)) + û = reshape(Matrix(A) \ vec(L), size(B)) + g = range(0, 1; length = 41) + (B, û, maximum(abs(evaluate(B, û, (x, y)) - exact2((x, y))) for x in g, y in g)) +end + +for n in (4, 8, 16) + println("n = ", rpad(n, 3), " error = ", poisson2d(3, n)[3]) +end +``` + +```@example gal +B7, û7, err7 = poisson2d(3, 16) +gg = range(0, 1; length = 101) +Z7 = [evaluate(B7, û7, (x, y)) for x in gg, y in gg] + +fig = Figure(size = (780, 320)) +ax1 = Axis3(fig[1, 1]; xlabel = "x₁", ylabel = "x₂", zlabel = "u", + title = "-Δu = 2π² sin πx₁ sin πx₂") +surface!(ax1, gg, gg, Z7) +ax2 = Axis(fig[1, 2]; xlabel = "x₁", ylabel = "x₂", aspect = 1, title = "error") +he = heatmap!(ax2, gg, gg, Z7 .- [exact2((x, y)) for x in gg, y in gg]; + colormap = :balance) +Colorbar(fig[1, 3], he) +fig +``` + +The Dirichlet condition holds on all four sides, again by construction: + +```@example gal +maximum(abs, [evaluate(B7, û7, pt) + for pt in ((0.0, 0.37), (1.0, 0.37), (0.42, 0.0), (0.42, 1.0))]) +``` + +Note what was *not* done here: the ``\mathbb{A}`` above was assembled densely, because the +example solves it with a dense factorisation. That is fine at ``n = 16`` and is not how a real +solver would work — the Laplacian is a sum of Kronecker products, so it is amenable to a +fast diagonalisation or a matrix-free Krylov method. The **mass** matrix, by contrast, is never +assembled at all: + +```@example gal +op7 = mass_operator(TensorProductQuadrature(B7)) +size(op7), map(f -> nameof(typeof(f)), mass_factors(op7)) +``` + +## 8. Depositing particles onto a basis + +The reverse of evaluation: given weights at scattered positions, accumulate them into +coefficients. This is the loop [`evaluate_all!`](@ref) and [`findcell`](@ref) exist for, and it +allocates nothing per particle. + +```@example gal +using Random +Random.seed!(42) + +bdep = BSplineBasis(UniformMesh(32, 0 .. 1), 3, Periodic()) +qdep = SplineQuadrature(bdep) + +positions = rand(20_000) +weights = fill(1 / 20_000, 20_000) + +function deposit!(coeffs, b, positions, weights) + fill!(coeffs, 0) + scratch = zeros(local_width(b)) + for (x, w) in zip(positions, weights) + j₀ = evaluate_all!(scratch, b, x) + for t in eachindex(scratch) + i = basis_index(b, j₀ + t - 1) + 1 ≤ i ≤ nbasis(b) || continue + coeffs[i] += w * scratch[t] + end + end + return coeffs +end + +load = deposit!(zeros(nbasis(bdep)), bdep, positions, weights) +sum(load), (@allocated deposit!(load, bdep, positions, weights)) +``` + +The sum is exactly the total weight, by the partition of unity, and the whole sweep over +20 000 particles allocates only the one scratch buffer. What comes out is the *load vector* +``\int f \phi_i``, so a density needs one mass solve: + +```@example gal +density = mass_operator(qdep) \ load +fdep = Spline(bdep, density) +mean_density = 1.0 # uniform on [0,1], total weight 1 +maximum(abs(fdep(x) - mean_density) for x in range(0, 1; length = 401)) +``` + +```@example gal +fig = Figure(size = (780, 340)) +xs = range(0, 1; length = 401) +ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "f(x)", + title = "20 000 uniform particles deposited onto a cubic periodic basis, n = 32") +lines!(ax, xs, fdep.(xs); label = "reconstructed density") +hlines!(ax, [1.0]; color = :black, linestyle = :dash, label = "exact") +Legend(fig[2, 1], ax; orientation = :horizontal, framevisible = false) +fig +``` + +The wobble is Monte-Carlo noise, not a discretisation error: the deposition itself is exact, +and the sum above is the total weight to the last bit. Its root-mean-square size is +``\sqrt{n/N_p}`` — about four per cent for ``20\,000`` particles in ``32`` cells, with the +maximum over the domain a few times that. Refining the basis makes it worse rather than +better, since it puts fewer particles in each cell. diff --git a/docs/src/index.md b/docs/src/index.md index 4186e83..1b717a5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,13 +4,27 @@ CurrentModule = SimpleSplines # SimpleSplines.jl -Periodic 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. It provides one basis — the periodic B-spline basis of -arbitrary degree on a uniform or non-uniform mesh — and one assembly table built on it, from -which the mass, stiffness, derivative and variable-coefficient matrices all follow as single -weighted contractions. +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`](@ref)), periodic ([`PeriodicBSplineBasis`](@ref)), and recombined ([`RecombinedBSplineBasis`](@ref)) — of any degree ``p \ge 0`` | +| **Meshes** | [`UniformMesh`](@ref), [`GradedMesh`](@ref), [`RandomMesh`](@ref), [`GeneralMesh`](@ref) | +| **Boundary conditions** | [`Free`](@ref), [`Dirichlet`](@ref), [`Neumann`](@ref), [`Natural`](@ref), [`Robin`](@ref), and the general [`Constraint`](@ref) — per end, plus [`Periodic`](@ref) for the whole basis | +| **Assembly** | [`SplineQuadrature`](@ref) and, from it, [`mass_matrix`](@ref), [`stiffness_matrix`](@ref), [`derivative_matrix`](@ref), [`mixed_matrix`](@ref), [`weighted_matrix`](@ref), [`basis_integrals`](@ref) | +| **Mass solves** | [`CirculantMass`](@ref) (FFT), [`BandedMass`](@ref) (banded Cholesky), [`FactorizedMass`](@ref) (sparse Cholesky), [`KroneckerMass`](@ref) (factored, ``D``-dimensional) | +| **Tensor products** | [`TensorProductBasis`](@ref) and [`TensorProductQuadrature`](@ref) in any number of dimensions, with degree, mesh, domain and boundary condition **per axis** | +| **Functions** | [`Spline`](@ref), [`SplineDerivative`](@ref), [`l2_projection`](@ref) | +| **Particle path** | [`findcell`](@ref) and [`evaluate_all!`](@ref), allocation-free, for depositing onto the basis | ## Installation @@ -19,113 +33,84 @@ using Pkg Pkg.add(url = "https://github.com/JuliaDEC/SimpleSplines.jl") ``` -## Basic usage +## In thirty seconds -A basis is a [`Mesh`](@ref) and a degree: +A basis is a [`Mesh`](@ref), a degree, and a [`BoundaryCondition`](@ref): ```@example intro using SimpleSplines -b = PeriodicBSplineBasis(UniformMesh(16, 2π), 3) -nbasis(b), degree(b), order(b) -``` - -Note that the number of degrees of freedom is the number of *cells*, `N = n`, not `n + p`: -on a torus the `p` extra functions of the bounded case are the ones the clamping introduces -at the two ends, and the periodic wrap identifies them in pairs. - -Basis functions are indexed as `b[x, j]`, and derivatives of arbitrary order come from -[`evaluate`](@ref): - -```@example intro -b[0.7, 3], evaluate(b, 3, 0.7, 1), evaluate(b, 3, 0.7, 2) +b = BSplineBasis(UniformMesh(16, 0 .. 1), 3, Dirichlet()) +nbasis(b), degree(b), order(b), polynomial_reproduction(b) ``` Assembly goes through a [`SplineQuadrature`](@ref), which tabulates the basis and its -derivatives at the global Gauß-Legendre points: +derivatives at the global Gauß-Legendre points and hands back every matrix built on them: ```@example intro q = SplineQuadrature(b) M = mass_matrix(q) -S = derivative_matrix(q) -maximum(abs, S + S') # antisymmetric on a periodic mesh +K = stiffness_matrix(q) +size(M), size(K) ``` -Projecting a function onto the spline space: +Fitting a function to the space is an ``L^2`` projection, and the result is callable: ```@example intro -û = l2_projection(q, sin) -abs(evaluate(b, û, 1.0) - sin(1.0)) +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 ``` -## Meshes - -Three families are provided, and which one to use is a question about what is being tested -rather than about the discretisation: - -| mesh | refining `n` gives | use it for | -|:--|:--|:--| -| [`UniformMesh`](@ref) | a uniform refinement | ordinary computation; but its assemblies are *circulant*, so it confirms some identities for the wrong reason | -| [`GradedMesh`](@ref) | a genuine mesh family | convergence rates | -| [`RandomMesh`](@ref) | an unrelated mesh | properties that must hold on *any* mesh | - -## The mass matrix - -Solves against the mass matrix go through a [`MassOperator`](@ref), and which representation -is built is decided by the **basis** — a uniform mesh is necessary for circulance but not -sufficient, since a clamped basis on one has ``p`` boundary functions at each end that are -not translates of anything: - - - for a [`PeriodicBSplineBasis`](@ref) on a [`UniformMesh`](@ref) the basis functions are - translates of a single cardinal spline, so ``\mathbb{M}`` is **circulant** and - diagonalised by the discrete Fourier transform. A solve is two planned transforms and a - pointwise division, and allocates nothing beyond its result. This is a - [`CirculantMass`](@ref). - - for a bounded basis — [`BSplineBasis`](@ref) or [`RecombinedBSplineBasis`](@ref) — there - is no seam, so the overlaps are contiguous and ``\mathbb{M}`` is **banded** outright. A - banded Cholesky solves it in ``O(Np)`` and allocates nothing at all. This is a - [`BandedMass`](@ref). - - for a periodic basis on a [`GradedMesh`](@ref) or [`RandomMesh`](@ref) it is banded - modulo ``N`` but *not* circulant — the basis functions are no longer translates of one - another — so there is nothing for a transform to diagonalise, and the wrap-around entries - put it outside the banded representation too. A sparse Cholesky factorisation is what is - left, a [`FactorizedMass`](@ref). +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: ```@example intro -mass_operator(q) +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)) ``` -The construction *verifies* circulance rather than assuming it: 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. - -## Storage - -The basis tabulation is **sparse**. A basis function is supported on `p+1` cells, so only -`(p+1)·nq` of the `n·nq` entries in its row are structurally nonzero, and storing ``\Phi`` -densely would make every contraction ``\Phi \, \mathrm{diag}(fw) \, \Phi^T`` cost -``O(N^2 n n_q)`` instead of ``O(N p^2 n_q)``. +## How the pieces fit -The constant assemblies — [`mass_matrix`](@ref), [`stiffness_matrix`](@ref), -[`derivative_matrix`](@ref) and any [`mixed_matrix`](@ref) — are memoised on first use, since -they do not depend on the field but are asked for inside every Newton iteration of a -downstream time integrator. - -## Quadrature - -The default rule, [`quadrature_order`](@ref), integrates degree ``3p-1`` exactly — what the -consistency of a Galerkin discretisation of a quadratic nonlinearity needs. Two lower -thresholds are worth knowing: - - - the mass matrix needs degree ``2p``, i.e. `nq ≥ p+1`; - - the antisymmetry of [`derivative_matrix`](@ref) needs degree ``2p-1``, i.e. `nq ≥ p`, - and fails outright below it on a non-uniform mesh. - -## Library - -```@index ``` - -```@autodocs -Modules = [SimpleSplines] +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`](@ref)`(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](@ref)** — one worked example from a bare `using` to a solved problem. + - **Theory** — [B-Splines](@ref) for the recursion, the derivative formula and the dimension + counts; [Boundary Conditions](@ref theory-boundary) for recombination; [Tensor + Products](@ref theory-tensorproduct) for the Kronecker structure. + - **Usage** — [Meshes](@ref usage-meshes), [Bases](@ref usage-bases), [Boundary + Conditions](@ref usage-boundary), [Assembly](@ref usage-assembly), [Tensor Products](@ref + usage-tensorproduct) and [Splines](@ref usage-splines) give the constructors, the + accessors and the traps. + - **[Gallery](@ref)** — eight solved problems with their measured errors. + - **[Library](@ref)** — every exported name. + +## Related packages + +[CompactBasisFunctions.jl](https://github.com/JuliaGNI/CompactBasisFunctions.jl) provides the +`Basis` hierarchy these bases join, along with the Lagrange, Chebyshev, Legendre and Bernstein +families; [QuadratureRules.jl](https://github.com/JuliaGNI/QuadratureRules.jl) provides the +Gauß-Legendre nodes and weights; and `GeometricBase` owns the shared accessors +[`basis`](@ref), [`degree`](@ref), [`nodes`](@ref), `nnodes` and [`order`](@ref), so that one +generic function per accessor is extended across the ecosystem rather than one defined per +package. diff --git a/docs/src/library.md b/docs/src/library.md new file mode 100644 index 0000000..eaa5364 --- /dev/null +++ b/docs/src/library.md @@ -0,0 +1,148 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# Library + +The complete API of *SimpleSplines.jl*. See the [Tutorial](@ref) for how the pieces fit +together, the Theory pages for the mathematics, and the Usage pages for the details of each +object. + +Seven accessors are **shared with the rest of the ecosystem** rather than defined here. +[`basis`](@ref), [`degree`](@ref), [`nodes`](@ref), [`nnodes`](@ref) and [`order`](@ref) +belong to `GeometricBase`, [`grid`](@ref) to `ContinuumArrays`, and [`nbasis`](@ref) to +[CompactBasisFunctions.jl](https://github.com/JuliaGNI/CompactBasisFunctions.jl), so that one +generic function per accessor is extended across the packages instead of one being defined in +each. Four further names are **re-exported unchanged** so that `using SimpleSplines` is enough +to write a domain: `..`, `leftendpoint` and `rightendpoint` come from `IntervalSets` and +`DomainSets`, and `Basis` from `CompactBasisFunctions`; their documentation lives in those +packages. + +```@index +``` + +## Meshes + +```@docs +Mesh +UniformMesh +GradedMesh +RandomMesh +GeneralMesh +breakpoints +ncells +domain +domainlength +meshwidth +findcell +``` + +## Boundary conditions + +```@docs +BoundaryCondition +Free +Periodic +Dirichlet +Neumann +Natural +Robin +Constraint +boundary_conditions +constraint_coefficients +constraint_order +nconstraints +``` + +## Bases + +```@docs +AbstractBSplineBasis +BSplineBasis +PeriodicBSplineBasis +RecombinedBSplineBasis +recombination_matrix +``` + +### Accessors + +```@docs +nbasis +degree +order +mesh +knotvector +nodes +nnodes +grid +boundary +basis_index +local_width +local_indices +polynomial_reproduction +``` + +### Evaluation + +```@docs +evaluate +evaluate_all +evaluate_all! +BSplineDerivative +PeriodicBSplineDerivative +``` + +## Assembly + +```@docs +SplineQuadrature +quadrature_order +quadrature_nodes +quadrature_weights +basis_values +basis +mass_matrix +mass_factorization +stiffness_matrix +derivative_matrix +mixed_matrix +weighted_matrix +basis_integrals +l2_projection +l2_projection! +``` + +## Mass operators + +```@docs +MassOperator +CirculantMass +BandedMass +FactorizedMass +KroneckerMass +mass_operator +mass_factors +mass_solve! +``` + +## Tensor products + +```@docs +TensorProductBasis +⊗ +bases +TensorProductQuadrature +quadratures +quadrature_grid_size +quadrature_sample +contract +``` + +## Splines + +```@docs +Spline +SplineDerivative +coefficients +derivative +``` diff --git a/docs/src/references.bib b/docs/src/references.bib new file mode 100644 index 0000000..deeb799 --- /dev/null +++ b/docs/src/references.bib @@ -0,0 +1,89 @@ +@book{deBoor:2001:PracticalGuideToSplines, + title = {A Practical Guide to Splines}, + author = {de Boor, Carl}, + year = {2001}, + edition = {Revised}, + series = {Applied Mathematical Sciences}, + volume = {27}, + publisher = {Springer}, + address = {New York}, + isbn = {978-0-387-95366-3} +} + +@book{Schumaker:2007:SplineFunctions, + title = {Spline Functions: Basic Theory}, + author = {Schumaker, Larry L.}, + year = {2007}, + edition = {Third}, + series = {Cambridge Mathematical Library}, + publisher = {Cambridge University Press}, + address = {Cambridge}, + doi = {10.1017/CBO9780511618994}, + isbn = {978-0-521-70512-7} +} + +@book{PieglTiller:1997:NURBSBook, + title = {The {NURBS} Book}, + author = {Piegl, Les and Tiller, Wayne}, + year = {1997}, + edition = {Second}, + series = {Monographs in Visual Communication}, + publisher = {Springer}, + address = {Berlin}, + doi = {10.1007/978-3-642-59223-2}, + isbn = {978-3-540-61545-3} +} + +@book{Hoellig:2003:FiniteElementMethodsWithBSplines, + title = {Finite Element Methods with {B}-Splines}, + author = {H{\"o}llig, Klaus}, + year = {2003}, + series = {Frontiers in Applied Mathematics}, + volume = {26}, + publisher = {SIAM}, + address = {Philadelphia}, + doi = {10.1137/1.9780898717532}, + isbn = {978-0-89871-699-3} +} + +@book{Cottrell:2009:IsogeometricAnalysis, + title = {Isogeometric Analysis: Toward Integration of {CAD} and {FEA}}, + author = {Cottrell, J. Austin and Hughes, Thomas J. R. and Bazilevs, Yuri}, + year = {2009}, + publisher = {Wiley}, + address = {Chichester}, + doi = {10.1002/9780470749081}, + isbn = {978-0-470-74873-2} +} + +@book{Boyd:2001:ChebyshevFourier, + title = {Chebyshev and Fourier Spectral Methods}, + author = {Boyd, John P.}, + year = {2001}, + edition = {Second}, + publisher = {Dover}, + address = {Mineola, New York}, + isbn = {978-0-486-41183-1} +} + +@article{Cox:1972:NumericalEvaluationOfBSplines, + title = {The numerical evaluation of {B}-splines}, + author = {Cox, Maurice G.}, + year = {1972}, + journal = {IMA Journal of Applied Mathematics}, + volume = {10}, + number = {2}, + pages = {134--149}, + doi = {10.1093/imamat/10.2.134} +} + +@article{deBoor:1972:OnCalculatingWithBSplines, + title = {On calculating with {B}-splines}, + author = {de Boor, Carl}, + year = {1972}, + journal = {Journal of Approximation Theory}, + volume = {6}, + number = {1}, + pages = {50--62}, + doi = {10.1016/0021-9045(72)90080-9} +} diff --git a/docs/src/references.md b/docs/src/references.md new file mode 100644 index 0000000..edb70bc --- /dev/null +++ b/docs/src/references.md @@ -0,0 +1,24 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# References + +The standard references for spline theory are de Boor +[deBoor:2001:PracticalGuideToSplines](@cite) and Schumaker +[Schumaker:2007:SplineFunctions](@cite). The recursion this package evaluates is Cox's and de +Boor's [Cox:1972:NumericalEvaluationOfBSplines](@cite), +[deBoor:1972:OnCalculatingWithBSplines](@cite). Piegl and Tiller +[PieglTiller:1997:NURBSBook](@cite) treat the geometric-modelling side, which this package +deliberately does not; Höllig [Hoellig:2003:FiniteElementMethodsWithBSplines](@cite) and +Cottrell, Hughes and Bazilevs [Cottrell:2009:IsogeometricAnalysis](@cite) treat the +finite-element and isogeometric settings it is written for. Basis recombination as a device for +imposing boundary conditions is discussed by Boyd [Boyd:2001:ChebyshevFourier](@cite). + +Every identifier below was checked against Crossref. Two entries carry no DOI because none is +registered for the edition cited: de Boor's revised edition of 2001 — the only Crossref record +for that title is the 1978 first edition's eBook, whose landing page no longer resolves — and +Boyd's Dover second edition, for which only the 1989 Springer first edition is registered. + +```@bibliography +``` diff --git a/docs/src/theory/boundary.md b/docs/src/theory/boundary.md new file mode 100644 index 0000000..c5cc527 --- /dev/null +++ b/docs/src/theory/boundary.md @@ -0,0 +1,341 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Boundary Conditions](@id theory-boundary) + +Two quite different things go by the name *boundary condition*, and separating them is what +makes the general case tractable. + +```@example bnd +using SimpleSplines +using CairoMakie +using LinearAlgebra +using Random +Random.seed!(1234) +CairoMakie.activate!(type = "png") # hide +nothing # hide +``` + +## Two kinds of closure + +**How the knot vector is closed** fixes the spline *space* and its dimension. There are two +choices, and they are different constructions rather than variants of one: + + - [`Free`](@ref) repeats each end knot ``p+1`` times — the clamped knot vector — giving + ``N = n + p`` and a basis that is interpolatory at the ends. + - [`Periodic`](@ref) continues the breakpoints periodically and wraps the basis onto the + torus, giving ``N = n``. + +Periodicity is not a per-end setting. It *identifies* the two ends, so it applies to a whole +axis or not at all; there is no such thing as being periodic at the left end only, and +`boundary_conditions` rejects the attempt. + +**A homogeneous linear constraint at one end** is something else: a condition *on* the clamped +space. [`Dirichlet`](@ref) is ``u = 0``, [`Neumann`](@ref) is ``u' = 0``, +[`Natural`](@ref) is ``u'' = 0``, [`Robin`](@ref) is ``\alpha u + \beta u' = 0``, and +[`Constraint`](@ref) is the general + +```math +L u \big\vert_{\text{end}} = 0 , +\qquad L = \sum_{k=0}^{m} c_k \, D^k , +\qquad c_m \ne 0 , +``` + +of order ``m``. These are imposed by *recombination*, and each one costs exactly one degree of +freedom at the end it applies to. + +All seven are subtypes of [`BoundaryCondition`](@ref), and all of them are given their meaning +in one place — [`constraint_coefficients`](@ref), which returns the ``(c_0, c_1, \dots)`` above +— so that the recombination has a single implementation rather than one method per condition. + +```@example bnd +[(bc, constraint_coefficients(bc), constraint_order(bc), nconstraints(bc)) + for bc in (Free(), Periodic(), Dirichlet(), Neumann(), Natural(), + Robin(2.0, 3.0), Constraint(1, 0, -2))] +``` + +## Recombination + +Take a clamped basis ``\varphi_1, \dots, \varphi_{N_p}`` of degree ``p`` and a condition of +order ``m`` at the left end ``a``. Put + +```math +a_i = (L \varphi_i)(a) . +``` + +A clamped knot vector has ``D^k \varphi_i(a) = 0`` for ``i > k+1``, so only the first ``m+1`` +functions can violate the condition, and of those only ``\varphi_{m+1}`` contributes to + +```math +a_{m+1} = c_m \, D^m \varphi_{m+1}(a) , +``` + +which is nonzero whenever the leading coefficient ``c_m`` is. The ``m`` combinations + +```math +\psi_t = \varphi_t - \frac{a_t}{a_{m+1}} \, \varphi_{m+1} , +\qquad t = 1, \dots, m , +``` + +are annihilated by ``L`` by construction, and ``\varphi_{m+2}, \varphi_{m+3}, \dots`` pass +through unchanged because they already satisfy the condition. The right end is the mirror +image, anchored on ``\varphi_{N_p - m}``. + +**The accounting.** ``m+1`` functions go in and ``m`` come out, so a constrained end costs one +degree of freedom *whatever the order of its condition*. That is worth stating separately from +the disjointness requirement, because the two are easy to conflate: `Natural()` reaches three +functions and removes one. + + - ``m = 0`` ([`Dirichlet`](@ref)): there are no combinations at all and ``\varphi_1`` is + simply dropped — the textbook elimination. + - ``m = 1`` ([`Neumann`](@ref)): one genuine combination of ``\varphi_1`` and + ``\varphi_2``. + - ``m = 2`` ([`Natural`](@ref)): ``\varphi_1`` and ``\varphi_2`` already have vanishing + second derivative at ``a`` and pass through; ``\varphi_3`` is the one that goes. + +The transformation is stored as a sparse matrix ``R`` with ``\psi_j = \sum_i R_{ij} +\varphi_i``, of size ``N_p \times N``, available as [`recombination_matrix`](@ref): + +```@example bnd +bclamped = BSplineBasis(UniformMesh(6, 0 .. 1), 3) +b = BSplineBasis(UniformMesh(6, 0 .. 1), 3, Neumann()) +R = recombination_matrix(b) +size(R), nbasis(bclamped), nbasis(b) +``` + +```@example bnd +Matrix(R)[1:4, 1:3] +``` + +The first column is ``\psi_1 = \varphi_1 - (a_1/a_2)\varphi_2``, and here ``a_1/a_2 = -1``: +the two end derivatives of a clamped basis are ``\mp p/h``, so the Neumann combination is +``\varphi_1 + \varphi_2``. That the anchor coefficient is itself often exactly ``1`` is why +the construction records which parent row each column carries with unit coefficient rather +than searching a column for the value ``1``. + +What the recombined functions look like, against the parent functions they are built from: + +```@example bnd +xs = range(0, 1; length = 801) +fig = Figure(size = (780, 560)) +for (row, bc) in enumerate((Dirichlet(), Neumann(), Natural())) + br = BSplineBasis(UniformMesh(6, 0 .. 1), 3, bc) + ax1 = Axis(fig[row, 1]; ylabel = string(nameof(typeof(bc))), + title = row == 1 ? "parent (clamped)" : "") + for j in 1:4 + lines!(ax1, xs, [evaluate(bclamped, j, x) for x in xs]) + end + ax2 = Axis(fig[row, 2]; title = row == 1 ? "recombined" : "") + for j in 1:4 + lines!(ax2, xs, [evaluate(br, j, x) for x in xs]) + end + for ax in (ax1, ax2) + xlims!(ax, 0, 0.7) + vlines!(ax, breakpoints(bclamped); color = (:black, 0.2), linestyle = :dash) + end +end +fig +``` + +The Dirichlet basis is the parent's with ``\varphi_1`` removed; the Neumann one has a first +function with a horizontal tangent at ``x = 0``; the Natural one has a first function whose +curvature vanishes there. + +```@example bnd +[(nameof(typeof(bc)), + evaluate(BSplineBasis(UniformMesh(6, 0 .. 1), 3, bc), 1, 0.0, d)) + for (bc, d) in ((Dirichlet(), 0), (Neumann(), 1), (Natural(), 2))] +``` + +Each is zero to round-off, and — the point of doing it this way — it is zero for *every* +coefficient vector, because the condition holds function by function: + +```@example bnd +bd = BSplineBasis(UniformMesh(6, 0 .. 1), 3, Dirichlet()) +û = randn(nbasis(bd)) +evaluate(bd, û, 0.0), evaluate(bd, û, 1.0) +``` + +## Every assembly is the parent's, conjugated + +Because ``\psi_j = \sum_i R_{ij} \varphi_i`` is linear, every bilinear form on the recombined +basis is the parent's form conjugated by ``R``: + +```math +\tilde{\mathbb{M}} = R^{\mathsf T} \mathbb{M} R , +\qquad +\tilde{\Phi} = R^{\mathsf T} \Phi , +\qquad +\tilde{A}_{a,b} = R^{\mathsf T} A_{a,b} R . +``` + +This is why [`SplineQuadrature`](@ref) needs no separate implementation for the recombined +case, and it is checked rather than assumed: + +```@example bnd +q = SplineQuadrature(bd) +qp = SplineQuadrature(parent(bd)) +Rd = recombination_matrix(bd) +maximum(abs, Matrix(mass_matrix(q)) - Matrix(Rd' * mass_matrix(qp) * Rd)) +``` + +``R`` has at most two nonzeros per column, so ``R^{\mathsf T} \mathbb{M} R`` is still banded +with half-bandwidth ``p`` up to the end blocks, and the mass matrix stays symmetric positive +definite. Every solve is therefore unaffected, which is the practical reason to prefer +recombination to the alternatives: + +| approach | banding | dimension | +|:--|:--|:--| +| **recombination** | preserved | ``N`` = the constrained dimension | +| row/column elimination | preserved, but only for ``u = 0`` | same | +| Lagrange multipliers | destroyed — a dense saddle-point block | ``N_p + \#`` constraints | +| penalty | preserved, but the condition holds only approximately and conditioning degrades | ``N_p`` | + +Basis recombination as a general device for imposing boundary conditions on a spectral or +spline basis is discussed by Boyd [Boyd:2001:ChebyshevFourier](@cite); the spline case is +treated in [Hoellig:2003:FiniteElementMethodsWithBSplines](@cite). + +## What is lost + +Two properties the clamped basis has are not guaranteed to survive recombination, and which +of them actually goes depends on the condition. Measured, on cubics over six cells: + +```@example bnd +xs2 = range(0, 1; length = 401) +for bc in (Free(), Dirichlet(), Neumann(), Natural(), Robin(1.5, -0.7)) + br = BSplineBasis(UniformMesh(6, 0 .. 1), 3, bc) + least = minimum(evaluate(br, j, x) for j in eachindex(br), x in xs2) + defect = maximum(abs, [sum(evaluate(br, j, x) for j in eachindex(br)) - 1 + for x in xs2]) + println(rpad(nameof(typeof(bc)), 11), " min ψ = ", rpad(round(least; digits = 4), 9), + " |Σψ - 1| = ", defect) +end +``` + +**The partition of unity survives exactly when the condition annihilates the constants**, that +is when ``c_0 = 0``. The reason is a row sum: ``\sum_j \psi_j = \sum_i (\sum_j R_{ij}) +\varphi_i``, and the anchor row contributes ``-\sum_{t \le m} a_t / a_{m+1}``, which equals one +precisely when ``\sum_i a_i = (L\mathbb{1})(a) = c_0`` vanishes. So [`Neumann`](@ref) and +[`Natural`](@ref) keep it, while [`Dirichlet`](@ref) loses it outright — the defect is ``1``, +since every function vanishes at the ends and the constant is simply not in the space — and a +[`Robin`](@ref) condition with ``\alpha \ne 0`` loses it partially. + +**Non-negativity** can go whichever way. ``\psi_t = \varphi_t - (a_t/a_{m+1})\varphi_{m+1}`` +is a difference, so nothing guarantees it, and above it is `Natural` that loses it while +`Dirichlet`, `Neumann` and that particular `Robin` happen to keep it. A scheme that needs a +projected density to stay non-negative should check rather than assume. + +The consequence is not cosmetic. A conservation law survives a discretisation only if the +conserved density lies in the *span* of the basis, and [`polynomial_reproduction`](@ref) is +the number that says which do: + +```math +\text{polynomial\_reproduction} = +\min_{\text{ends}} \big( \text{order of the lowest derivative the condition involves} \big) - 1 , +``` + +with an unconstrained end contributing ``p``. The derivation is short: the recombined space is +the parent space intersected with the two conditions, and in the local coordinate ``x - a`` one +has ``D^k u(a) = k!\, b_k``, so ``L = \sum_k c_k D^k`` reads off the coefficients +``b_0, \dots, b_m`` one at a time and annihilates the whole of ``\mathbb{P}_m`` exactly when +``c_0 = \dots = c_m = 0``. + +| condition | lowest derivative in ``L`` | `polynomial_reproduction` | +|:--|:--|--:| +| [`Free`](@ref) | — | ``p`` | +| [`Periodic`](@ref) | — | ``0`` — the constants only, since ``x`` is not periodic | +| [`Dirichlet`](@ref) | ``u`` | ``-1`` | +| [`Neumann`](@ref), `Robin(0, β)` | ``u'`` | ``0`` | +| [`Natural`](@ref) | ``u''`` | ``1`` | +| `Robin(α, β)` with ``\alpha \ne 0`` | ``u`` | ``-1`` | +| `Constraint(0, 0, 0, 1)` | ``u'''`` | ``2`` | +| a pair `(left, right)` | — | the minimum of the two ends | + +```@example bnd +m = UniformMesh(8, 0 .. 1) +[(bc, polynomial_reproduction(BSplineBasis(m, 4, bc))) + for bc in (Free(), Periodic(), Dirichlet(), Neumann(), Natural(), + Robin(1.0, 1.0), Robin(0.0, 1.0), Constraint(0, 0, 0, 1), + (Dirichlet(), Neumann()))] +``` + +Two entries there are worth a second look. `Robin(1.0, 1.0)` reproduces nothing — ``-1``, the +same as `Dirichlet` — because what matters is the *lowest* derivative the condition involves, +and a Robin condition with ``\alpha \ne 0`` involves ``u`` itself: no non-zero constant +satisfies ``u + u' = 0`` at an endpoint. And `Robin(0.0, 1.0)` is `Neumann` in disguise, so it +keeps the constants; the named type is preferable where it applies, since it says what is +meant and `Dirichlet` in addition takes the cheaper elimination path. + +!!! warning "The boundary condition is not a free choice" + Choosing [`Dirichlet`](@ref) because the solution happens to decay at the edge of the + domain removes the constants from the span, and with them the conservation of mass. The + decay is a property of the solution; imposing it on the *space* is a different and much + stronger statement. A scheme whose conservation proof needs ``1``, ``v`` and ``v^2`` in the + span needs `polynomial_reproduction(b) ≥ 2`, and can assert it. + +## The periodic alternative + +Where the geometry allows it, the periodic closure gives up nothing at the boundary because +there is no boundary. Every basis function spans ``p+1`` cells, the basis is +``\mathcal{C}^{p-1}`` across the seam ``a \equiv b`` as well as inside, and integration by +parts over ``\Omega`` leaves no boundary terms. That last point is what makes +[`derivative_matrix`](@ref) exactly antisymmetric on a periodic mesh: + +```math +\int_\Omega \partial_x (\phi_k \phi_l) \, \mathrm{d}x = 0 +\quad \Longrightarrow \quad +S_{kl} + S_{lk} = 0 . +``` + +```@example bnd +qper = SplineQuadrature(BSplineBasis(UniformMesh(16, 0 .. 2π), 3, Periodic())) +S = derivative_matrix(qper) +maximum(abs, S + S') +``` + +On a clamped basis the same integral is ``[\phi_k \phi_l]_a^b``, which does not vanish, and +``S`` is not antisymmetric: + +```@example bnd +qcl = SplineQuadrature(BSplineBasis(UniformMesh(16, 0 .. 2π), 3)) +maximum(abs, derivative_matrix(qcl) + derivative_matrix(qcl)') +``` + +What the periodic closure costs is polynomial reproduction beyond the constants. + +## Two things a condition cannot be + +**Nonlocal.** Every condition here involves the solution at one endpoint only, which is what +keeps ``R`` sparse and the mass matrix banded. A multi-point or integral condition — +``u(a) = u(b)`` other than through [`Periodic`](@ref), or ``\int_\Omega u = 0`` — would need +the nullspace of a dense constraint matrix, and the banding that every assembly here relies on +would be gone. This is deliberately out of scope. + +**Of order higher than the degree.** A condition involving ``D^k`` with ``k > p`` constrains +nothing, since the ``k``-th derivative of a degree-``p`` spline vanishes identically. That is +reported rather than silently accepted: + +```@example bnd +try + BSplineBasis(UniformMesh(8, 0 .. 1), 1, Natural()) +catch err + err +end +``` + +## A caveat on ill-scaled Robin conditions + +The recombination anchors each end block on ``\varphi_{m+1}``, the only one of the ``m+1`` +candidates that contributes to ``a_{m+1} = c_m D^m \varphi_{m+1}(a)``. That anchor is nonzero +whenever the leading coefficient ``c_m`` is — but *how* nonzero is the caller's business, not +the construction's. + +`Robin(1.0, 1e-20)` puts a vanishing coefficient on the highest derivative. The anchor is then +numerically zero and the recombination coefficients ``-a_t / a_{m+1}`` become huge, and the +failure is silent in an unpleasant way: the mass matrix comes out finite with a condition number +already `Inf`, `cholesky(…; check = false)` reports success on it anyway, and the projection +that follows looks plausible. Prefer the named type where one applies — +`Robin(1.0, 0.0)` is [`Dirichlet`](@ref) and `Robin(0.0, 1.0)` is [`Neumann`](@ref) — and keep +the two coefficients of a genuine Robin condition within a few orders of magnitude of each +other. diff --git a/docs/src/theory/bsplines.md b/docs/src/theory/bsplines.md new file mode 100644 index 0000000..1942692 --- /dev/null +++ b/docs/src/theory/bsplines.md @@ -0,0 +1,347 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# B-Splines + +What a B-spline basis is, how this package computes it, and which of its properties the rest +of the manual leans on. The standard references are de Boor +[deBoor:2001:PracticalGuideToSplines](@cite) and Schumaker +[Schumaker:2007:SplineFunctions](@cite); Piegl and Tiller +[PieglTiller:1997:NURBSBook](@cite) is the geometric-modelling counterpart, and Höllig +[Hoellig:2003:FiniteElementMethodsWithBSplines](@cite) and Cottrell, Hughes and Bazilevs +[Cottrell:2009:IsogeometricAnalysis](@cite) treat the finite-element setting this package is +written for. + +```@example theory +using SimpleSplines +using CairoMakie +using LinearAlgebra +CairoMakie.activate!(type = "png") # hide +nothing # hide +``` + +## The spline space + +Fix a closed interval ``\Omega = [a,b]`` and a subdivision of it into ``n`` cells by +breakpoints + +```math +a = y_1 < y_2 < \dots < y_{n+1} = b . +``` + +The spline space of degree ``p`` on that subdivision is + +```math +\mathcal{S}^p + = \left\{ u \in \mathcal{C}^{p-1}(\Omega) : + u \big\vert_{[y_k, y_{k+1}]} \in \mathbb{P}_p \ \text{ for } k = 1, \dots, n \right\} , +``` + +the piecewise polynomials of degree ``p`` that are ``p-1`` times continuously differentiable +across every interior breakpoint. Counting: ``n(p+1)`` coefficients, less ``p`` continuity +conditions at each of the ``n-1`` interior breakpoints, leaves + +```math +\dim \mathcal{S}^p = n(p+1) - (n-1)p = n + p . +``` + +That is the dimension the clamped basis realises. The other two closures in this package +change the space rather than the basis, and change the count with it — see [Dimension +counting](@ref) below. + +## Knot vectors + +A B-spline basis is not built from the breakpoints directly but from a *knot vector*, a +non-decreasing sequence in which a value may repeat. Multiplicity is what encodes reduced +continuity: a knot of multiplicity ``\mu`` drops the smoothness there to +``\mathcal{C}^{p-\mu}``. + +The **clamped** (or open) knot vector repeats each end knot ``p+1`` times, + +```math +\big[\, \underbrace{a, \dots, a}_{p+1},\ y_2, \dots, y_n,\ + \underbrace{b, \dots, b}_{p+1} \,\big] , +``` + +of length ``n + 2p + 1``. The continuity at the two ends is thereby ``\mathcal{C}^{-1}``, +which is what lets a basis function be nonzero *at* the endpoint at all; inside, every knot is +simple and the space is ``\mathcal{C}^{p-1}`` as it should be. + +```@example theory +b = BSplineBasis(UniformMesh(5, 0 .. 1), 2) +knotvector(b) +``` + +The **periodic** closure does something different. The right endpoint is dropped — +``y_{n+1} \equiv y_1`` on a torus — and the remaining breakpoints are continued periodically, +``y_{i+n} = y_i + L`` with ``L = b - a``, so that the recursion runs on a bi-infinite sequence +with no repeated knots anywhere. The resulting splines satisfy +``\phi_{j+n}(x) = \phi_j(x - L)``, so only ``n`` of them are distinct on ``\Omega``, and the +periodic basis is their wrapping onto the torus, + +```math +\phi_j \big\vert_\Omega (x) = \sum_{r \in \mathbb{Z}} \phi_j(x + rL) , +\qquad 1 \le j \le n , +``` + +a finite sum since each ``\phi_j`` is supported on ``p+1`` cells. In practice five periodic +images of the breakpoints are stored, which is what keeps every index the recursion reaches +inside the array: + +```@example theory +c = BSplineBasis(UniformMesh(5, 0 .. 1), 2, Periodic()) +knotvector(c) +``` + +## The Cox-de Boor recursion + +Given a knot vector ``(x_j)``, the B-splines of degree ``0`` are the indicators of the knot +spans and those of degree ``p`` are built from those of degree ``p-1`` +[Cox:1972:NumericalEvaluationOfBSplines](@cite), [deBoor:1972:OnCalculatingWithBSplines](@cite): + +```math +\phi_j^0 = \mathbb{1}_{[x_j, x_{j+1})} , +\qquad +\phi_j^p = w_j^p \, \phi_j^{p-1} + \left(1 - w_{j+1}^p\right) \phi_{j+1}^{p-1} , +\qquad +w_j^p = \frac{x - x_j}{x_{j+p} - x_j} . +``` + +A knot span of zero length contributes nothing, which is what makes the recursion well +defined at a repeated knot: the guard is on the span being positive, since it is the division +and not the numerator that fails. + +Raising the degree by one widens the support by one span and adds one order of smoothness. +The four degrees on one mesh: + +```@example theory +fig = Figure(size = (780, 480)) +xs = range(0, 1; length = 801) +for (i, p) in enumerate(0:3) + ax = Axis(fig[(i - 1) ÷ 2 + 1, (i - 1) % 2 + 1]; + title = "p = $(p), N = $(nbasis(BSplineBasis(UniformMesh(5, 0 .. 1), p)))", + xlabel = i > 2 ? "x" : "") + bp = BSplineBasis(UniformMesh(5, 0 .. 1), p) + for j in eachindex(bp) + lines!(ax, xs, bp[xs, j]) + end + vlines!(ax, breakpoints(bp); color = (:black, 0.2), linestyle = :dash) +end +fig +``` + +`SimpleSplines` computes single values from the recursion exactly as written — see the +internal `_bspline`, which the test suite uses as the reference — and the whole local block by +de Boor's triangular scheme, which is ``O(p^2)`` for all ``p+1`` nonzero functions together +rather than per function. [`evaluate`](@ref) is the first path, [`evaluate_all!`](@ref) the +second. + +## Properties + +Four properties, each of which something in this package relies on. + +**Local support.** ``\phi_j^p`` is supported on ``[x_j, x_{j+p+1}]``, which is ``p+1`` +consecutive cells. So at most ``p+1`` basis functions are nonzero at any point, at most +``2p+1`` functions overlap any given one, and every assembled matrix is banded with +half-bandwidth ``p``. This is the property that makes the whole package cheap: the tabulation +is sparse, the mass matrix is banded, and a particle deposition costs ``O(p^2)`` rather than +``O(N)``. + +```@example theory +b3 = BSplineBasis(UniformMesh(8, 0 .. 1), 3) +count(j -> abs(b3[0.3, j]) > 0, eachindex(b3)), local_width(b3) +``` + +**Non-negativity.** ``\phi_j^p \ge 0`` everywhere. Together with the next property this makes +a spline with non-negative coefficients non-negative, which is why a projected distribution +function stays a distribution function on a clamped or periodic basis. + +**Partition of unity.** ``\sum_j \phi_j^p \equiv 1`` on ``\Omega``. This is what makes the +constants lie in the span, hence ``\int_\Omega u_h = \mathbf{1}^{\mathsf T} \mathbb{M} \hat{u}`` +and [`basis_integrals`](@ref)` == mass_matrix(q) * ones(N)`. + +```@example theory +maximum(abs, [sum(b3[x, :]) - 1 for x in range(0, 1; length = 401)]) +``` + +**Smoothness.** ``\phi_j^p \in \mathcal{C}^{p-1}`` across an interior breakpoint, and no +better: the ``p``-th derivative jumps. Both halves matter — the first is why a +``\mathcal{C}^2`` cubic basis can carry a third-order operator after one integration by parts, +the second is why it cannot carry one without. + +```@example theory +b4 = BSplineBasis(UniformMesh(6, 0 .. 1), 3) +y = breakpoints(b4)[3] +ε = 1e-7 +[(d, abs(evaluate(b4, 3, y - ε, d) - evaluate(b4, 3, y + ε, d))) for d in 0:3] +``` + +For ``d \le 2`` the difference across the breakpoint is proportional to ``\varepsilon``, which +is what a continuous derivative probed at ``y \pm \varepsilon`` gives; at ``d = 3`` it is +``O(1)``, and that is the jump. + +```@example theory +fig = Figure(size = (780, 300)) +xs2 = range(0, 1; length = 1201) +for d in 0:2 + ax = Axis(fig[1, d + 1]; xlabel = "x", title = "D^$(d) φ₃, p = 3") + lines!(ax, xs2, [evaluate(b4, 3, x, d) for x in xs2]) + vlines!(ax, breakpoints(b4); color = (:black, 0.2), linestyle = :dash) +end +fig +``` + +## Dimension counting + +The number of degrees of freedom is a property of the closure, not of the degree alone. +[`nbasis`](@ref) reports it. + +| basis | ``N`` | why | +|:--|:--|:--| +| [`BSplineBasis`](@ref) | ``n + p`` | the dimension of ``\mathcal{S}^p`` computed above | +| [`PeriodicBSplineBasis`](@ref) | ``n`` | on a torus there are only ``n`` breakpoints and no clamping, so the ``p`` extra functions of the bounded case are the ones the clamping introduced; the wrap identifies them in pairs | +| [`RecombinedBSplineBasis`](@ref) | ``n + p - \#\{\text{constrained ends}\}`` | each homogeneous condition is one linear functional on ``\mathcal{S}^p``, hence one dimension, whatever the order of the derivative it involves | + +```@example theory +m = UniformMesh(8, 0 .. 1) +[(bc, nbasis(BSplineBasis(m, 3, bc))) + for bc in (Free(), Periodic(), Dirichlet(), Neumann(), Natural(), + (Dirichlet(), Free()))] +``` + +A periodic basis also needs ``n > p``: a basis function spans ``p+1`` cells, so with +``n \le p`` it would wrap onto itself and the sum over periodic images would not terminate. +The clamped basis has no such requirement — a single cell carries the whole of +``\mathbb{P}_p``. + +## Derivatives + +Differentiating the recursion once gives a recursion of its own, + +```math +\frac{\mathrm{d}}{\mathrm{d}x} \phi_j^p + = p \left( \frac{\phi_j^{p-1}}{x_{j+p} - x_j} + - \frac{\phi_{j+1}^{p-1}}{x_{j+p+1} - x_{j+1}} \right) , +``` + +so the derivative of a degree-``p`` B-spline is a difference of two degree-``(p-1)`` +B-splines. Iterating, ``D^d \phi_j^p`` is a combination of ``d+1`` B-splines of degree +``p-d``, and for ``d > p`` it vanishes identically. + +Two consequences worth stating plainly: + + - **Any derivative order is available**, as the fourth positional argument of + [`evaluate`](@ref). This is not the same as stacking lazy `Derivative` products: `b'` gives + a [`BSplineDerivative`](@ref), which is first order only and does not compose, whereas + `evaluate(b, j, x, 3)` is the third derivative directly. A third-order operator needs the + latter. + - **``d > p`` returns zero rather than an error.** The ``p``-th derivative of a degree-``p`` + spline is piecewise constant and everything above it is zero, so a loop over derivative + orders needs no special case. + +```@example theory +b1 = BSplineBasis(UniformMesh(4, 0 .. 1), 1) +[evaluate(b1, 2, 0.3, d) for d in 0:3] +``` + +Internally, [`evaluate_all!`](@ref) does not apply the derivative recursion to a finished +block. It runs de Boor's scheme at degree ``q = p - d`` — cheaper, since the block is +smaller — and then *lifts* the result to degree ``p`` by applying the recursion above ``d`` +times, the block growing by one function at each step. The two routes agree to round-off, and +the test suite checks that at every degree and order. + +## Greville abscissae and polynomial reproduction + +The natural interpolation points of a B-spline basis are the *Greville abscissae*, the +averages of the ``p`` interior knots of each function, + +```math +\xi_j = \frac{1}{p} \sum_{i=1}^{p} x_{j+i} . +``` + +Each ``\xi_j`` lies in the support of ``\phi_j``, and for the clamped basis the collocation +matrix ``\phi_j(\xi_i)`` is invertible — the Schoenberg-Whitney condition. [`nodes`](@ref) +returns them, and `ContinuumArrays`' `grid` is the same thing. + +```@example theory +nodes(BSplineBasis(UniformMesh(6, 0 .. 1), 3)) +``` + +```@example theory +bd = BSplineBasis(UniformMesh(6, 0 .. 1), 3) +C = [evaluate(bd, j, ξ) for ξ in nodes(bd), j in eachindex(bd)] +cond(C) +``` + +A separate question is which polynomials lie in the *span*. +[`polynomial_reproduction`](@ref) answers it with a number, and the answer is not always +``p``: the periodic closure keeps only the constants, because ``x`` is not periodic, and a +Dirichlet-recombined basis keeps nothing at all, because every one of its functions vanishes +at the ends. This decides which moments an ``L^2`` projection preserves, so it is not a matter +of taste — see [Boundary Conditions](@ref theory-boundary). + +## Two conventions to know + +**Cells are half-open, except the last.** ``\phi^0_j`` is the indicator of +``[x_j, x_{j+1})``, which makes every B-spline right-continuous and makes the partition of +unity hold on ``[a,b)`` but fail at ``b``, where every indicator is false. That would be the +wrong answer for a clamped basis, which is meant to be interpolatory at both ends. The +topmost span is therefore closed on the right, and only there: the value at ``b`` is exactly +the limit from the left, and ``\phi_N(b) = 1``. + +```@example theory +bc3 = BSplineBasis(UniformMesh(4, 0 .. 1), 3) +sum(bc3[1.0, :]), bc3[1.0, nbasis(bc3)] +``` + +**Outside the domain, a bounded basis evaluates to zero.** Not an error, and not an +extrapolated polynomial — zero, which is the mathematically correct value of a compactly +supported function. This is deliberate, and it is the behaviour a particle method wants: a +particle that leaves the domain stops contributing. It is also a trap worth knowing, since +nothing is raised. + +```@example theory +evaluate(bc3, degree(bc3) + 1, -0.5), sum(bc3[1.5, :]) +``` + +A [`PeriodicBSplineBasis`](@ref) is different: any real argument is accepted and reduced onto +the domain first, so the result is the periodic extension and no point is ever outside. + +```@example theory +cp = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Periodic()) +evaluate(cp, 3, 0.4), evaluate(cp, 3, 0.4 + 5.0), sum(cp[7.3, :]) +``` + +## Cost + +| operation | cost | +|:--|:--| +| one basis function at a point, [`evaluate`](@ref) | ``O(2^p)`` — the recursion as written, unmemoised | +| all nonzero functions at a point, [`evaluate_all!`](@ref) | ``O(p^2)`` for the whole block, no allocation | +| a spline at a point, `evaluate(b, û, x)` | ``O(p^2)``, not ``O(N)`` — it goes through the block | +| assembling ``\int f D^a\phi_k D^b\phi_l``, [`weighted_matrix`](@ref) | ``O(N p^2 n_q)`` on the sparse tabulation, against ``O(N^2 n n_q)`` on a dense one | +| a mass solve | ``O(N \log N)`` circulant, ``O(Np)`` banded | + +The first row is worth dwelling on, because it is a deliberate choice and not an oversight. +[`evaluate`](@ref) runs the Cox-de Boor recursion literally, splitting into two subproblems at +every level with no memoisation, so that the assembly built on it doubles as a check of the +formula. That costs ``O(2^p)`` per value. [`evaluate_all!`](@ref) uses de Boor's triangular +scheme, which is ``O(p^2)`` for *all* ``p+1`` functions together — so the block form is not +merely `p+1` times cheaper, it has a different exponent: + +| ``p`` | `evaluate`, one function | `evaluate_all`, the whole block | ratio | +|--:|--:|--:|--:| +| 3 | 26 ns | 21 ns | 1.3 | +| 6 | 241 ns | 39 ns | 6.1 | +| 8 | 963 ns | 58 ns | 16.7 | +| 12 | 15.5 µs | 104 ns | 149 | + +Absolute timings are a property of the machine; the exponents are not. Doubling the degree +from ``6`` to ``12`` costs the single value 64 times more and the block 2.6 times more, which +is ``2^6`` against ``(12/6)^2``. The measurement is archived as +`scripts/evaluate_cost_scaling.jl`, which fails if the two stop being distinguishable. + +At the cubics almost everything here uses the difference hardly matters; at ``p = 12`` it is a +factor of 150. Either way, a loop that wants every nonzero function at a point should ask for +the block. diff --git a/docs/src/theory/tensorproduct.md b/docs/src/theory/tensorproduct.md new file mode 100644 index 0000000..942b1e7 --- /dev/null +++ b/docs/src/theory/tensorproduct.md @@ -0,0 +1,245 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Tensor Products](@id theory-tensorproduct) + +On a box ``\Omega = \Omega_1 \times \dots \times \Omega_D`` the natural spline space is the +tensor product of ``D`` one-dimensional spaces. Everything about it follows from the +one-dimensional case, but the *bookkeeping* does not — and the bookkeeping is where a +hand-rolled tensor-product spline reliably goes wrong. + +```@example tp +using SimpleSplines +using CairoMakie +using LinearAlgebra +using Random +Random.seed!(1234) +CairoMakie.activate!(type = "png") # hide +nothing # hide +``` + +## The space + +Given one basis per axis, the product basis functions are the products + +```math +\Phi_{i_1 \dots i_D}(x) = \prod_{d=1}^{D} \phi^{(d)}_{i_d}(x_d) , +\qquad 1 \le i_d \le N_d , +``` + +and they span + +```math +\mathcal{S}^{p_1} \otimes \dots \otimes \mathcal{S}^{p_D} , +\qquad +N = \dim = \prod_{d=1}^{D} N_d . +``` + +Each factor is an independent basis, so **degree, mesh, domain and boundary condition are all +per axis**, and nothing is shared between them: + +```@example tp +B = BSplineBasis(UniformMesh(12, 0 .. 2π), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(9, -10 .. 10), 4, Dirichlet()) +ndims(B), size(B), nbasis(B), degree(B), boundary(B) +``` + +The dimension count is per axis too: ``12`` for the periodic cubic axis (``n``), and +``9 + 4 - 2 = 11`` for the Dirichlet quartic one (``n + p`` less one per constrained end). + +```@example tp +xs = range(0, 1; length = 61) +Bp = BSplineBasis(UniformMesh(5, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(5, 0 .. 1), 2) +Z = [evaluate(Bp, (4, 3), (x, y)) for x in xs, y in xs] +fig = Figure(size = (760, 320)) +ax1 = Axis3(fig[1, 1]; xlabel = "x₁", ylabel = "x₂", zlabel = "Φ", + title = "Φ₄,₃ = φ⁽¹⁾₄ ⋅ φ⁽²⁾₃") +surface!(ax1, xs, xs, Z) +ax2 = Axis(fig[1, 2]; xlabel = "x₁", ylabel = "x₂", aspect = 1, + title = "its support: 4 × 3 cells") +heatmap!(ax2, xs, xs, Z) +vlines!(ax2, breakpoints(bases(Bp)[1]); color = (:white, 0.5)) +hlines!(ax2, breakpoints(bases(Bp)[2]); color = (:white, 0.5)) +fig +``` + +A product basis function is supported on a box of ``\prod_d (p_d + 1)`` cells, so that is how +many are nonzero at any point — ``64`` for cubics in three dimensions, against +``N = \prod_d N_d`` in total. + +## Coefficients are an array, and the flattening matters + +A spline in this basis is ``u_h = \sum_I \hat{u}_I \Phi_I`` with ``\hat{u}`` a +``D``-dimensional array of size `size(B)`. That is the shape in which the Kronecker structure +of every operator is visible, and it is what [`l2_projection`](@ref) returns and +[`Spline`](@ref) holds. + +Sooner or later the array has to be flattened — to hand it to a linear solver, or to a time +integrator written against vectors. The convention here is that the **first axis varies +fastest**, which is `vec`, `LinearIndices(B)`, and Julia's column-major order: + +```@example tp +LinearIndices(Bp)[1:3, 1:3] +``` + +Three things have to agree on that convention or the results are silently wrong: the +evaluation, the mass matrix and the projection. They do, and the visible consequence is the +*reversed* factor order in the assembled Kronecker product: + +```math +\mathbb{M} = \mathbb{M}^{(D)} \otimes \dots \otimes \mathbb{M}^{(1)} , +``` + +with axis ``1`` innermost. This is not a choice that can be made differently in one place: + +```@example tp +Q = TensorProductQuadrature(Bp) +op = mass_operator(Q) +M1, M2 = map(mass_matrix, mass_factors(op)) +maximum(abs, Matrix(op) - kron(Matrix(M2), Matrix(M1))) +``` + +```@example tp +û = randn(size(Bp)...) +maximum(abs, Matrix(op) * vec(û) - vec(op * û)) +``` + +## The Kronecker structure is used, not just noted + +The mass matrix is never assembled. Its inverse factorises exactly — + +```math +\mathbb{M}^{-1} = \left( \mathbb{M}^{(D)} \right)^{-1} \otimes \dots \otimes + \left( \mathbb{M}^{(1)} \right)^{-1} +``` + +— which is an identity, not an approximate splitting, so a solve is ``D`` sweeps of +one-dimensional solves, each factor keeping whatever representation it earned on its own axis: + +```@example tp +map(f -> nameof(typeof(f)), mass_factors(mass_operator( + TensorProductQuadrature(B)))) +``` + +```@example tp +maximum(abs, op \ (op * û) - û) +``` + +What that buys, for a cubic basis with ``N_d = n`` per axis: + +| | dense Cholesky of ``\mathbb{M}`` | [`KroneckerMass`](@ref) | +|:--|:--|:--| +| storage | ``O(N^2) = O(n^{2D})`` | ``O\!\left(\sum_d N_d p_d\right) = O(D n p)`` | +| factorisation | ``O(N^3)`` | ``D`` banded factorisations, ``O(D n p^2)`` | +| one solve | ``O(N^2)`` | ``O(N \sum_d p_d)`` | + +At ``D = 2`` and ``n = 41`` that is a ``1681 \times 1681`` dense factorisation avoided; at +``D = 3`` it is a ``68921 \times 68921`` one, which does not fit. The factored form is the +only reason a three-dimensional velocity space is tractable at all. + +```@example tp +[(D, prod(ntuple(_ -> 41, D)), prod(ntuple(_ -> 41, D))^2 * 8 / 2^30) for D in 1:3] +``` + +The third column is the memory, in gibibytes, that the dense mass matrix of a +``41``-per-axis basis would need. + +## The load vector factorises too + +An ``L^2`` projection needs ``L_I = \int_\Omega f \, \Phi_I``, and since ``\Phi_I`` is a +product the quadrature over the ``D``-dimensional grid contracts one axis at a time. With +``\Phi_k`` the sparse ``N_k \times Q_k`` tabulation of axis ``k``, + +```math +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} , +``` + +which [`contract`](@ref) evaluates as ``D`` sparse matrix products, reshaping so that the axis +being contracted is first and cycling it to the back afterwards. After ``D`` steps the axes are +back in order. There is no ``D``-dimensional tabulation anywhere. + +**The integrand need not be separable — only the basis is, and that is enough.** This is the +point most easily missed. A non-separable ``f`` is sampled on the product grid and contracted; +nothing about the method assumes it factorises: + +```@example tp +g(x) = sin(x[1]) * exp(-x[2]^2 / 8) + cos(2x[1]) * x[2] / 20 # not separable +errs = Float64[] +pts = [(2π * rand(), 20rand() - 10) for _ in 1:400] +for n in (8, 16, 32) + Bn = BSplineBasis(UniformMesh(n, 0 .. 2π), 3, Periodic()) ⊗ + BSplineBasis(UniformMesh(n, -10 .. 10), 3) + Qn = TensorProductQuadrature(Bn) + ĝ = l2_projection(Qn, g) + push!(errs, maximum(abs(evaluate(Bn, ĝ, x) - g(x)) for x in pts)) +end +errs, errs[1] / errs[2], errs[2] / errs[3] +``` + +Cubics give fourth-order convergence, so the error should fall by roughly ``2^4 = 16`` per +refinement. + +## Derivatives are a multi-index + +There is no scalar notion of "the derivative order" on a product basis: +``\partial_1 \Phi`` and ``\partial_2 \Phi`` are different objects. The order is therefore a +per-axis tuple, and + +```math +D^{d} \Phi_I = \prod_{k=1}^{D} D^{d_k} \phi^{(k)}_{i_k} , +``` + +so the gradient component ``\partial_k`` is `d = ntuple(i -> i == k ? 1 : 0, D)`. For a +[`Spline`](@ref), `derivative(s, k::Integer)` builds that tuple. + +```@example tp +b1, b2 = bases(Bp) +evaluate(Bp, (4, 3), (0.3, 0.4), (1, 2)) ≈ + evaluate(b1, 4, 0.3, 1) * evaluate(b2, 3, 0.4, 2) +``` + +## What separates and what does not + +| quantity | separable? | +|:--|:--| +| the mass matrix ``\int \Phi_I \Phi_J`` | yes — ``\mathbb{M}^{(D)} \otimes \dots \otimes \mathbb{M}^{(1)}`` | +| a pure derivative matrix ``\int D^a \Phi_I \, D^b \Phi_J`` | yes — one factor is a one-dimensional mixed matrix, the rest are mass matrices | +| the Laplacian ``\int \nabla\Phi_I \cdot \nabla\Phi_J`` | yes as a **sum** of ``D`` separable terms, ``\sum_k \mathbb{M} \otimes \dots \otimes \mathbb{K}^{(k)} \otimes \dots \otimes \mathbb{M}`` | +| the load vector ``\int f \, \Phi_I`` for arbitrary ``f`` | yes as a *contraction*, though not as a product | +| ``\int f \, D^a\Phi_I \, D^b\Phi_J`` with a non-separable ``f`` | **no** | + +The last row is where the tensor-product structure stops paying. Such a matrix is still +``\Phi_a \, \mathrm{diag}(f \odot w) \, \Phi_b^{\mathsf T}`` with a ``D``-dimensional +tabulation ``\Phi``, but no sequence of one-dimensional contractions gives it, so it needs the +Kronecker product of the one-dimensional tables formed explicitly. `SimpleSplines` does not +provide that; a package that needs it builds it from +[`basis_values`](@ref)`(quadratures(q)[k], d)` and pays the storage. + +## One boundary condition per axis + +The condition is a property of an axis. Passing a ``D``-tuple gives one condition per axis, +**not** the two ends of one axis: + +```@example tp +B2 = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet()) ⊗ + BSplineBasis(UniformMesh(8, 0 .. 1), 3, Neumann()) +boundary(B2), size(B2) +``` + +To impose different conditions at the two ends of a single axis, build that axis's basis with +a pair and then take the product: + +```@example tp +axis = BSplineBasis(UniformMesh(8, 0 .. 1), 3, (Dirichlet(), Neumann())) +boundary(axis), nbasis(axis) +``` + +[`polynomial_reproduction`](@ref) of the product is the minimum over the axes, because a +polynomial lies in the span only if its restriction to each axis lies in that axis's span. One +periodic or Dirichlet axis is enough to lose it: + +```@example tp +polynomial_reproduction(B2), polynomial_reproduction(Bp), polynomial_reproduction(B) +``` diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md new file mode 100644 index 0000000..8ce181b --- /dev/null +++ b/docs/src/tutorial.md @@ -0,0 +1,268 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# Tutorial + +One example, worked from a bare `using` to a solved boundary-value problem. Everything on this +page runs; the numbers and figures are produced when the manual is built. + +```@example tutorial +using SimpleSplines +using CairoMakie +CairoMakie.activate!(type = "png") # hide +nothing # hide +``` + +## A mesh + +A [`Mesh`](@ref) is a subdivision of a closed interval into cells, and nothing more — it +carries no degree and no boundary condition: + +```@example tutorial +m = UniformMesh(8, 0 .. 1) +ncells(m), domain(m), meshwidth(m) +``` + +```@example tutorial +breakpoints(m) +``` + +There are `n+1` breakpoints for `n` cells, both endpoints included. The domain may be written +as an interval `0 .. 1`, as a tuple `(0, 1)`, or — the common periodic case — as a single +number standing for ``[0, L]``, so `UniformMesh(8, 2π)` is ``[0, 2\pi]``. + +## A basis + +A basis is a mesh plus a degree. With no boundary condition given it is the *clamped* basis, +of dimension ``N = n + p``: + +```@example tutorial +b = BSplineBasis(m, 3) +nbasis(b), degree(b), order(b) +``` + +Basis functions are indexed `b[x, j]` — the point first, as for a matrix of samples — and +there are `p+1` of them nonzero at any point: + +```@example tutorial +fig = Figure(size = (760, 300)) +ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "φⱼ(x)", + title = "the 11 cubic B-splines of a clamped basis on 8 cells") +xs = range(0, 1; length = 601) +for j in eachindex(b) + lines!(ax, xs, b[xs, j]) +end +vlines!(ax, breakpoints(b); color = (:black, 0.25), linestyle = :dash) +fig +``` + +They are non-negative and sum to one at every point — a *partition of unity* — and the +clamping makes the first and last interpolatory at the ends: + +```@example tutorial +maximum(abs, [sum(b[x, :]) - 1 for x in xs]), b[0.0, 1], b[1.0, nbasis(b)] +``` + +## A spline + +A [`Spline`](@ref) is a basis together with the coefficients of one element of its span, and +it is callable: + +```@example tutorial +û = zeros(nbasis(b)) +û[4] = 1.0 +û[7] = -0.6 +s = Spline(b, û) +s(0.35) +``` + +A second argument is the derivative order, and [`derivative`](@ref) turns the derivative into +a callable of its own so that it broadcasts: + +```@example tutorial +fig = Figure(size = (760, 300)) +ax = Axis(fig[1, 1]; xlabel = "x", title = "a spline and its first two derivatives") +lines!(ax, xs, s.(xs); label = "s") +lines!(ax, xs, derivative(s, 1).(xs) ./ 10; label = "s′ / 10") +lines!(ax, xs, derivative(s, 2).(xs) ./ 100; label = "s″ / 100") +axislegend(ax; position = :rb) +fig +``` + +The second derivative is piecewise linear and the third would be piecewise constant: a +degree-``p`` spline is ``\mathcal{C}^{p-1}`` across an interior breakpoint, so exactly `p-1` +derivatives are continuous. + +## Fitting a function + +To put a given function into the space, project it. That needs a quadrature — a +[`SplineQuadrature`](@ref), which tabulates the basis and its derivatives at the +Gauß-Legendre points of every cell: + +```@example tutorial +q = SplineQuadrature(b) +length(quadrature_nodes(q)), sum(quadrature_weights(q)) +``` + +The weights sum to the length of the domain. [`l2_projection`](@ref) returns the coefficient +vector: + +```@example tutorial +g(x) = exp(-8 * (x - 0.35)^2) * sinpi(3x) +ĝ = l2_projection(q, g) +sg = Spline(b, ĝ) +maximum(abs(sg(x) - g(x)) for x in xs) +``` + +```@example tutorial +fig = Figure(size = (760, 320)) +ax1 = Axis(fig[1, 1]; xlabel = "x", title = "L² projection, p = 3, n = 8") +lines!(ax1, xs, g.(xs); label = "g") +lines!(ax1, xs, sg.(xs); linestyle = :dash, linewidth = 3, label = "projection") +axislegend(ax1; position = :rb) +ax2 = Axis(fig[1, 2]; xlabel = "x", title = "error") +lines!(ax2, xs, sg.(xs) .- g.(xs)) +vlines!(ax2, breakpoints(b); color = (:black, 0.2), linestyle = :dash) +fig +``` + +Refining the mesh reduces that error at order ``p+1``: + +```@example tutorial +for n in (8, 16, 32, 64) + bn = BSplineBasis(UniformMesh(n, 0 .. 1), 3) + qn = SplineQuadrature(bn) + ĝn = l2_projection(qn, g) + err = maximum(abs(evaluate(bn, ĝn, x) - g(x)) for x in xs) + println("n = ", lpad(n, 2), " error = ", err) +end +``` + +## Imposing a boundary condition + +The same construction takes a [`BoundaryCondition`](@ref) as a third argument, and returns +whichever of the three bases that condition calls for: + +```@example tutorial +bfree = BSplineBasis(m, 3, Free()) +bdir = BSplineBasis(m, 3, Dirichlet()) +bper = BSplineBasis(m, 3, Periodic()) +(typeof(bfree).name.name, nbasis(bfree)), +(typeof(bdir).name.name, nbasis(bdir)), +(typeof(bper).name.name, nbasis(bper)) +``` + +`Dirichlet()` costs one degree of freedom per end, and the condition ``u(0) = u(1) = 0`` then +holds for *every* coefficient vector — it is a property of the space, not something enforced +afterwards: + +```@example tutorial +using Random +Random.seed!(1234) +v̂ = randn(nbasis(bdir)) +evaluate(bdir, v̂, 0.0), evaluate(bdir, v̂, 1.0) +``` + +```@example tutorial +fig = Figure(size = (760, 300)) +ax = Axis(fig[1, 1]; xlabel = "x", ylabel = "ψⱼ(x)", + title = "the 9 recombined functions of the Dirichlet basis") +for j in eachindex(bdir) + lines!(ax, xs, [evaluate(bdir, j, x) for x in xs]) +end +vlines!(ax, breakpoints(bdir); color = (:black, 0.25), linestyle = :dash) +fig +``` + +Every one of them vanishes at both ends. Note that they are no longer a partition of unity: for +`Dirichlet` recombination just drops the one function at each end that fails to vanish there, so +the constants have left the space. [`polynomial_reproduction`](@ref) reports that as a number. + +```@example tutorial +polynomial_reproduction.((bfree, bdir, bper)) +``` + +## Solving a boundary-value problem + +With the boundary condition in the space, the Galerkin discretisation of +``-u'' = f`` on ``(0,1)`` with ``u(0) = u(1) = 0`` is the linear system +``\mathbb{K} \hat{u} = \mathbf{b}``, with ``\mathbb{K}`` the stiffness matrix and +``b_i = \int f \psi_i``: + +```@example tutorial +qd = SplineQuadrature(bdir) +K = stiffness_matrix(qd) + +f(x) = π^2 * sinpi(x) # so that u(x) = sin(πx) +rhs = basis_values(qd, 0) * (quadrature_weights(qd) .* f.(quadrature_nodes(qd))) +û = Matrix(K) \ rhs + +uh = Spline(bdir, û) +maximum(abs(uh(x) - sinpi(x)) for x in xs) +``` + +```@example tutorial +fig = Figure(size = (760, 300)) +ax = Axis(fig[1, 1]; xlabel = "x", title = "-u'' = π² sin(πx), u(0) = u(1) = 0") +lines!(ax, xs, sinpi.(xs); label = "exact") +scatter!(ax, nodes(bdir), uh.(nodes(bdir)); label = "solution at the nodes") +axislegend(ax; position = :cb) +fig +``` + +No rows were deleted and no values were substituted into the system: the basis has one +function fewer per constrained end, so `K` is already the matrix of the constrained problem. + +## Two dimensions + +A tensor product is the `⊗` of one-dimensional bases, each with its own degree, mesh and +boundary condition: + +```@example tutorial +B = BSplineBasis(UniformMesh(12, 0 .. 1), 3, Dirichlet()) ⊗ + BSplineBasis(UniformMesh(12, 0 .. 1), 3, Dirichlet()) +ndims(B), size(B), degree(B) +``` + +Coefficients are a `D`-dimensional array of size `size(B)`, and everything else follows the +one-dimensional pattern: + +```@example tutorial +Q = TensorProductQuadrature(B) +h(x) = sinpi(x[1]) * sinpi(2 * x[2]) * (1 - x[1]) +Ĥ = l2_projection(Q, h) +size(Ĥ), maximum(abs(evaluate(B, Ĥ, (x, y)) - h((x, y))) + for x in range(0, 1; length = 41), y in range(0, 1; length = 41)) +``` + +```@example tutorial +gr = range(0, 1; length = 121) +Z = [evaluate(B, Ĥ, (x, y)) for x in gr, y in gr] +fig = Figure(size = (760, 320)) +ax1 = Axis(fig[1, 1]; xlabel = "x₁", ylabel = "x₂", title = "projection", aspect = 1) +hm = heatmap!(ax1, gr, gr, Z) +Colorbar(fig[1, 2], hm) +ax2 = Axis(fig[1, 3]; xlabel = "x₁", ylabel = "x₂", title = "error", aspect = 1) +he = heatmap!(ax2, gr, gr, Z .- [h((x, y)) for x in gr, y in gr]; colormap = :balance) +Colorbar(fig[1, 4], he) +fig +``` + +The mass matrix of that basis is never assembled: it is the Kronecker product of the two +one-dimensional ones, and a solve is two sequences of one-dimensional solves. + +```@example tutorial +op = mass_operator(Q) +nameof(typeof(op)), map(f -> nameof(typeof(f)), mass_factors(op)), size(op) +``` + +## Where to go next + + - The mathematics: [B-Splines](@ref), [Boundary Conditions](@ref theory-boundary), + [Tensor Products](@ref theory-tensorproduct). + - The details of each object, its keyword arguments and its traps: [Meshes](@ref + usage-meshes), [Bases](@ref usage-bases), [Boundary Conditions](@ref usage-boundary), + [Assembly](@ref usage-assembly), [Tensor Products](@ref usage-tensorproduct), + [Splines](@ref usage-splines). + - More solved problems, with convergence rates: [Gallery](@ref). diff --git a/docs/src/usage/bases.md b/docs/src/usage/bases.md new file mode 100644 index 0000000..389ec5f --- /dev/null +++ b/docs/src/usage/bases.md @@ -0,0 +1,251 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Bases](@id usage-bases) + +Three concrete bases implement [`AbstractBSplineBasis`](@ref), and which one you get is +decided by the boundary condition rather than chosen by name. + +```@example bases +using SimpleSplines +``` + +## Construction + +One constructor covers all three: + +```julia +BSplineBasis(mesh, p) # clamped, the default +BSplineBasis(mesh, p, bc) # dispatches on bc +``` + +```@example bases +m = UniformMesh(8, 0 .. 1) +[(bc, nameof(typeof(BSplineBasis(m, 3, bc))), nbasis(BSplineBasis(m, 3, bc))) + for bc in (Free(), (Free(), Free()), Periodic(), :periodic, + Dirichlet(), (Dirichlet(), Free()), Neumann(), Natural())] +``` + +`Periodic()` gives a [`PeriodicBSplineBasis`](@ref), two `Free` ends give the plain clamped +[`BSplineBasis`](@ref) unwrapped, and anything else gives a +[`RecombinedBSplineBasis`](@ref). This is deliberate: it lets one call site select any of the +three, which is exactly what a tensor product with a different condition per axis needs. + +The named constructors are also available directly, and `PeriodicBSplineBasis` has a +convenience form that builds its own uniform mesh: + +```julia +BSplineBasis(mesh, p) BSplineBasis{T}(mesh, p) +PeriodicBSplineBasis(mesh, p) PeriodicBSplineBasis{T}(mesh, p) +PeriodicBSplineBasis(n, p; L = 2π) # = PeriodicBSplineBasis(UniformMesh(n, L), p) +RecombinedBSplineBasis(parent, left, right) +``` + +```@example bases +PeriodicBSplineBasis(16, 3) == BSplineBasis(UniformMesh(16, 2π), 3, Periodic()) +``` + +### What is required, and what is rejected + +| basis | requires | +|:--|:--| +| [`BSplineBasis`](@ref) | ``p \ge 0``, ``n \ge 1`` — a single cell carries the whole of ``\mathbb{P}_p`` | +| [`PeriodicBSplineBasis`](@ref) | ``p \ge 0`` and ``n > p``, so that a basis function does not wrap onto itself | +| [`RecombinedBSplineBasis`](@ref) | a **clamped** parent; each condition of order ``\le p``; the two end blocks disjoint; at least one function left | + +```@example bases +for bad in (() -> BSplineBasis(UniformMesh(8, 0 .. 1), -1), + () -> PeriodicBSplineBasis(UniformMesh(3, 0 .. 1), 3), + () -> BSplineBasis(UniformMesh(8, 0 .. 1), 1, Natural()), + () -> BSplineBasis(UniformMesh(1, 0 .. 1), 0, Dirichlet()), + () -> RecombinedBSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3), + Periodic(), Dirichlet())) + try + bad() + catch err + println(err.msg, "\n") + end +end +``` + +## Accessors + +```@example bases +b = BSplineBasis(UniformMesh(8, 0 .. 1), 3) +(nbasis(b), degree(b), order(b), ncells(b), local_width(b), + polynomial_reproduction(b), boundary(b)) +``` + +| call | returns | +|:--|:--| +| [`nbasis`](@ref)`(b)` | the dimension ``N``, and the length a coefficient vector must have | +| [`degree`](@ref)`(b)` | the polynomial degree ``p`` | +| [`order`](@ref)`(b)` | ``p + 1`` | +| [`ncells`](@ref)`(b)`, [`mesh`](@ref)`(b)`, [`domain`](@ref)`(b)`, [`domainlength`](@ref)`(b)`, [`meshwidth`](@ref)`(b)` | forwarded to the mesh | +| [`breakpoints`](@ref)`(b)` | the `n+1` cell boundaries, **cached in the basis** | +| [`knotvector`](@ref)`(b)` | the knot sequence the recursion runs on | +| [`boundary`](@ref)`(b)` | `Periodic()`, or the pair `(left, right)` | +| [`nodes`](@ref)`(b)` | the Greville abscissae; `ContinuumArrays`' `grid` is the same | +| `nnodes(b)` | `nbasis(b)` | +| [`local_width`](@ref)`(b)` | how many entries [`evaluate_all!`](@ref) writes: ``p+1``, or more near a recombined end | +| [`polynomial_reproduction`](@ref)`(b)` | the largest degree the span reproduces exactly, or `-1` | +| [`basis_index`](@ref)`(b, j)` | `j`, wrapped onto `1:N` where the basis is periodic | +| `eachindex(b)`, `axes(b)`, `eltype(b)` | `1:N`, `(Inclusion(domain), 1:N)`, the element type | + +!!! note "`order` here means `p + 1`" + That is the *spline* convention: a B-spline of order ``k`` is piecewise of degree + ``k - 1``. It is **not** the meaning `order` carries for the bases of + `CompactBasisFunctions`, where it is the number of basis functions. For a spline basis the + two differ, and `nbasis` is the one you want for a dimension. + +## Evaluating + +```julia +evaluate(b, j, x, d = 0) # basis function j +evaluate(b, û, x, d = 0) # the spline with coefficients û +b[x, j] b[x, :] b[X, j] b[X, :] b(x, j) +b' # a BSplineDerivative: first order only +``` + +`x` may be a number or a vector; with a vector the result is a vector (or, for `b[X, :]`, a +matrix of points × functions). + +```@example bases +xs = [0.1, 0.5, 0.9] +evaluate(b, 3, 0.3), evaluate(b, 3, 0.3, 2), size(b[xs, :]), b'[0.3, 3] +``` + +!!! warning "The argument order is not the same in the two spellings" + [`evaluate`](@ref) takes the **index first**, `evaluate(b, j, x, d)`, because the + derivative order comes last. `getindex` and the callable form take the **point first**, + `b[x, j]` and `b(x, j)`, to match a matrix of samples whose rows are points. Both spellings + are correct; mixing them up is silent whenever `j` and `x` are both plausible numbers. + +### Derivatives of arbitrary order + +The fourth argument of [`evaluate`](@ref) is the derivative order, any non-negative integer. +`b'` is a [`BSplineDerivative`](@ref), which is a lazy product and **first order only** — the +products do not compose, so a second derivative is not `b''`. + +```@example bases +[evaluate(b, 3, 0.3, d) for d in 0:4] +``` + +Orders above `p` are identically zero and are returned as such rather than raised, so a loop +over derivative orders needs no special case. A negative order is an error. + +### Outside the domain + +A bounded basis evaluates to **zero** outside its closed domain — not an error, and not an +extrapolated polynomial. That is the mathematically correct value for a compactly supported +function, and it is the behaviour a particle method wants; it is also silent, so it is worth +knowing. + +```@example bases +sum(b[1.7, :]), evaluate(b, 4, -0.2) +``` + +A [`PeriodicBSplineBasis`](@ref) instead accepts any real argument and reduces it onto the +domain, so the result is the periodic extension: + +```@example bases +c = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Periodic()) +sum(c[1.7, :]), evaluate(c, 4, 0.7) ≈ evaluate(c, 4, 0.7 + 3.0) +``` + +### At the right endpoint + +Knot spans are half-open, except the topmost, which is closed on the right so that the +clamped basis is interpolatory at ``b`` as it should be. The value there is exactly the limit +from the left. + +```@example bases +b[1.0, nbasis(b)], sum(b[1.0, :]) +``` + +## The local path + +A loop that needs *every* nonzero function at a point — a matrix assembly, a particle +deposition — should not call [`evaluate`](@ref) once per index. That path runs the recursion +literally, at ``O(2^p)`` per value; [`evaluate_all!`](@ref) fills a caller-supplied buffer with +the whole local block by de Boor's triangular scheme, in ``O(p^2)`` for all `p+1` together, and +allocates nothing: + +```@example bases +buf = zeros(local_width(b)) +j₀ = evaluate_all!(buf, b, 0.3) +j₀, buf, sum(buf) +``` + +```@example bases +all(buf[t] ≈ evaluate(b, basis_index(b, j₀ + t - 1), 0.3) for t in eachindex(buf)) +``` + +Three rules for using it: + +1. **The buffer must have exactly [`local_width`](@ref)`(b)` entries.** That is `p+1` for a + clamped or periodic basis, and possibly more for a recombined one, where a function near an + end spans the union of two parent supports. A wrong length is a `DimensionMismatch`. +2. **`j₀` is the first index *before* wrapping**, because the block is contiguous only in that + form. On a periodic basis it can fall outside `1:N`; put every index through + [`basis_index`](@ref) rather than writing `mod1` at the call site, since whether the wrap is + needed is a property of the basis and getting it wrong on a bounded basis silently folds the + two ends of the domain together. +3. **Outside a bounded domain the buffer is zeroed** and the *clamped* cell's `j₀` is + returned, so a deposition loop touches the same block it would for a point just inside and + adds nothing to it. + +```@example bases +buf2 = zeros(local_width(c)) +jc = evaluate_all!(buf2, c, 0.05) # near the periodic seam +jc, [basis_index(c, jc + t - 1) for t in eachindex(buf2)] +``` + +The allocating form [`evaluate_all`](@ref) returns `(j₀, values)` with a fresh buffer, which +is convenient at a call site that runs once and wasteful in a loop. + +A deposition therefore reads: + +```@example bases +coeffs = zeros(nbasis(b)) +particles = [(0.13, 1.0), (0.42, 2.0), (0.97, 0.5), (1.4, 3.0)] # the last is outside +scratch = zeros(local_width(b)) +for (x, w) in particles + j = evaluate_all!(scratch, b, x) + for t in eachindex(scratch) + i = basis_index(b, j + t - 1) + 1 ≤ i ≤ nbasis(b) || continue + coeffs[i] += w * scratch[t] + end +end +coeffs, sum(coeffs) +``` + +The sum is ``1.0 + 2.0 + 0.5 = 3.5`` and not ``6.5``: the partition of unity means each +particle deposits its whole weight, and the one outside the domain deposited nothing. + +## Cell lookup + +[`findcell`](@ref)`(b, x)` gives the cell index in `1:ncells`, reading the breakpoints cached +in the basis so that it allocates nothing; on a uniform mesh it is a division that does not +touch them at all. [`local_indices`](@ref)`(b, cell)` gives the indices of the functions +nonzero on a cell — `cell:(cell+p)` for a clamped basis, the same block wrapped for a periodic +one, and a precomputed range for a recombined one. + +```@example bases +findcell(b, 0.3), collect(local_indices(b, 3)), collect(local_indices(c, 1)) +``` + +Note the periodic case: cell `1` is reached by functions `6, 7, 8, 1`, wrapped. + +## Comparison + +`==` compares type, degree, mesh and boundary condition; `isequal` additionally requires the +same element type; `isapprox` compares the meshes approximately. + +```@example bases +b == BSplineBasis(UniformMesh(8, 0 .. 1), 3), +b == BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet()), +b == c +``` diff --git a/docs/src/usage/boundary.md b/docs/src/usage/boundary.md new file mode 100644 index 0000000..75957ac --- /dev/null +++ b/docs/src/usage/boundary.md @@ -0,0 +1,223 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Boundary Conditions](@id usage-boundary) + +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](@ref theory-boundary) +in the Theory section. + +```@example ubnd +using SimpleSplines +using LinearAlgebra +using SparseArrays +``` + +## The seven types + +| type | condition | order ``m`` | constructor | +|:--|:--|:--|:--| +| [`Free`](@ref) | none — the plain clamped basis | ``-1`` | `Free()` | +| [`Periodic`](@ref) | the periodic closure of the whole axis | ``-1`` | `Periodic()` | +| [`Dirichlet`](@ref) | ``u = 0`` | ``0`` | `Dirichlet()` | +| [`Neumann`](@ref) | ``u' = 0`` | ``1`` | `Neumann()` | +| [`Natural`](@ref) | ``u'' = 0`` | ``2`` | `Natural()` | +| [`Robin`](@ref) | ``\alpha u + \beta u' = 0`` | ``1``, or ``0`` if ``\beta = 0`` | `Robin(α, β)` | +| [`Constraint`](@ref) | ``\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. + +```@example ubnd +Robin(1, 2.0), Constraint(0, 0, 0, 1), constraint_order(Constraint(1, 0)) +``` + +## Three traits, one meaning + +Every condition is given its meaning in exactly one place, so that the recombination has a +single implementation: + +```@example ubnd +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 +``` + + - [`constraint_coefficients`](@ref) — the ``(c_0, c_1, \dots)``, or `nothing` for the two + conditions that impose none. + - [`constraint_order`](@ref) — the highest derivative appearing, or `-1`. Note + `constraint_order(Robin(1.0, 0.0)) == 0`: the trailing zero does not count. + - [`nconstraints`](@ref) — the degrees of freedom removed: `0` or `1`, never more, whatever + the order. + +## Writing a specification + +[`boundary_conditions`](@ref) 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 | + +```@example ubnd +(boundary_conditions(Dirichlet()), + boundary_conditions((Dirichlet(), Neumann())), + boundary_conditions(:periodic), + boundary_conditions((: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. + +```@example ubnd +try + boundary_conditions((Periodic(), Dirichlet())) +catch err + println(err.msg) +end +``` + +### 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. + +```@example ubnd +for s in (:Natural, :nothing, :Dirichlet, :quasiperiodic) + try + BoundaryCondition(s) + catch err + println(err.msg, "\n") + end +end +``` + +## Which basis a condition selects + +```@example ubnd +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 +``` + +Read it back with [`boundary`](@ref): + +```@example ubnd +(boundary(BSplineBasis(m, 3)), + boundary(BSplineBasis(m, 3, Periodic())), + boundary(BSplineBasis(m, 3, (Dirichlet(), Neumann())))) +``` + +## Working with a recombined basis + +A [`RecombinedBSplineBasis`](@ref) answers the whole [`AbstractBSplineBasis`](@ref) 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`](@ref) is the sparse ``R`` with +``\psi_j = \sum_i R_{ij}\varphi_i``. + +```@example ubnd +b = BSplineBasis(m, 3, Natural()) +bp = parent(b) +R = recombination_matrix(b) +nbasis(bp), nbasis(b), size(R), nnz(R) +``` + +Every assembly is the parent's conjugated by ``R``, which is worth knowing when comparing +against a hand-built reference: + +```@example ubnd +q, qp = SplineQuadrature(b), SplineQuadrature(bp) +maximum(abs, Matrix(mass_matrix(q)) - Matrix(R' * mass_matrix(qp) * R)) +``` + +**[`local_width`](@ref) 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!`](@ref) must have exactly that many entries, and how many depends on the +condition rather than on the degree alone: + +```@example ubnd +[(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))] +``` + +```@example ubnd +degree(b) + 1, local_width(b), +[length(collect(local_indices(b, k))) for k in 1:ncells(b)] +``` + +```@example ubnd +try + evaluate_all!(zeros(degree(b) + 1), b, 0.3) +catch err + println(err.msg) +end +``` + +So `local_width(b)`, never `degree(b) + 1`, is what a buffer should be sized by. + +**[`nodes`](@ref) 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. + +```@example ubnd +nodes(b) +``` + +**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. + +```@example ubnd +[(bc, nbasis(BSplineBasis(m, 4, bc))) + for bc in (Dirichlet(), Neumann(), Natural(), Constraint(0, 0, 0, 1))] +``` + +## 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`](@ref). 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`](@ref) 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`](@ref)` != mass_matrix(q) * ones(N)`. Check + [`polynomial_reproduction`](@ref) before assuming a conservation law survives. + - **An ill-scaled [`Robin`](@ref) 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. + +```@example ubnd +bd = BSplineBasis(m, 3, Dirichlet()) +qd = SplineQuadrature(bd) +maximum(abs, basis_integrals(qd) - mass_matrix(qd) * ones(nbasis(bd))) +``` + +That is the size of the discrepancy the partition of unity would have made zero. diff --git a/docs/src/usage/meshes.md b/docs/src/usage/meshes.md new file mode 100644 index 0000000..d4491bd --- /dev/null +++ b/docs/src/usage/meshes.md @@ -0,0 +1,187 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Meshes](@id usage-meshes) + +A [`Mesh`](@ref) is a subdivision of a closed interval ``[a,b]`` into `n` cells by `n+1` +breakpoints. It is geometry alone: it carries no degree and no boundary condition, and the +same mesh serves a clamped, a recombined and a periodic basis. + +```@example mesh +using SimpleSplines +``` + +## The four families + +| mesh | breakpoints | refining `n` gives | use it for | +|:--|:--|:--|:--| +| [`UniformMesh`](@ref) | ``y_i = a + (i-1)h`` | a uniform refinement | ordinary computation — but its periodic assemblies are *circulant*, so it confirms some identities for the wrong reason | +| [`GradedMesh`](@ref) | the image of a uniform mesh under one fixed smooth map | a genuine mesh family | measuring convergence rates | +| [`RandomMesh`](@ref) | perturbed cell widths from a seeded stream | an *unrelated* mesh | properties that must hold on **any** mesh | +| [`GeneralMesh`](@ref) | whatever you pass | whatever you pass | a subdivision none of the three describes | + +The distinction between the first three is about what is being *tested*, not about the +discretisation. A convergence rate measured across a family of `RandomMesh`es means nothing, +because each `n` gives a different mesh rather than a refinement of the previous one; that is +what `GradedMesh` is for. Conversely, an identity checked only on a `UniformMesh` may hold +there because the assemblies are circulant and fail everywhere else, which is what +`RandomMesh` catches. + +```@example mesh +[(nameof(typeof(m)), ncells(m), round(meshwidth(m); digits = 4)) + for m in (UniformMesh(8, 0 .. 1), GradedMesh(8, 0 .. 1), RandomMesh(8, 0 .. 1))] +``` + +## Constructors + +```julia +UniformMesh(n, domain) # UniformMesh(n) is [0, 2π] +GradedMesh(n, domain; amplitude = 0.12) +RandomMesh(n, domain; seed = 1, spread = 0.6) +GeneralMesh(breakpoints) +``` + +Each also has a `{T}` form — `UniformMesh{Float32}(8, 0 .. 1)` — for fixing the breakpoint +element type explicitly. + +### Three ways to write the domain + +A `ClosedInterval`, a two-tuple, or a bare number standing for ``[0, L]``. The last is kept +because `UniformMesh(16, 2π)` is the common periodic case. + +```@example mesh +UniformMesh(8, 0 .. 1) == UniformMesh(8, (0, 1)) == UniformMesh(8, 1) +``` + +```@example mesh +domain(UniformMesh(8)), domain(UniformMesh(4, -10 .. 10)) +``` + +`..` is re-exported, so `using SimpleSplines` is enough to write an interval; so are +`leftendpoint` and `rightendpoint`. + +```@example mesh +leftendpoint(domain(UniformMesh(4, -10 .. 10))), +rightendpoint(domain(UniformMesh(4, -10 .. 10))) +``` + +### Element types + +An **integer** domain is promoted to a float, because the breakpoints of a subdivision of it +are not integers — otherwise `UniformMesh(4, 0 .. 1)` would try to store `0.25` in an `Int` +and throw an `InexactError` from a long way away. `Rational` and extended-precision types are +left alone, since they are closed under the division `breakpoints` performs and narrowing them +would discard the exactness that is the reason for using them. + +```@example mesh +eltype(UniformMesh(4, 0 .. 1)), +eltype(UniformMesh(4, 0 .. 1//1)), +eltype(UniformMesh(4, big(0.0) .. big(1.0))) +``` + +## Accessors + +| call | returns | +|:--|:--| +| [`breakpoints`](@ref)`(m)` | the `n+1` cell boundaries, increasing, both endpoints included | +| [`ncells`](@ref)`(m)` | `n` | +| [`domain`](@ref)`(m)` | the `ClosedInterval` ``[a,b]`` | +| [`domainlength`](@ref)`(m)` | ``b - a``; named this way because it is the period ``L`` of a periodic basis | +| [`meshwidth`](@ref)`(m)` | the **largest** cell width, the ``h`` a convergence rate is measured against | +| `eltype(m)`, `length(m)`, `first(m)`, `last(m)` | element type, `ncells`, and the two endpoints | + +```@example mesh +m = GeneralMesh([-20.0; range(-10, 10; length = 5); 20.0]) +ncells(m), domain(m), meshwidth(m), diff(breakpoints(m)) +``` + +`meshwidth` is the largest cell, not the average: for that mesh it is `10.0`, not `40/6`. + +!!! warning "`breakpoints` must not be mutated" + For [`GradedMesh`](@ref), [`RandomMesh`](@ref) and [`GeneralMesh`](@ref) this returns the + mesh's *own* array rather than a copy, which is what keeps [`findcell`](@ref) off the + allocator on the innermost loop of a particle deposition. Writing to it corrupts the mesh, + and silently: the domain and cell count are unchanged so nothing rejects the result, while + `hash` and `==` now report a different mesh. [`UniformMesh`](@ref) computes its + breakpoints from a closed form and hands back a fresh vector, but that is an + implementation detail, not a licence. Take a `copy` if you need to modify it. + +## Equality is geometric, `isequal` is not + +`==` and `hash` compare the domain and the breakpoints. `isequal` also compares the *type*, +because the mesh type selects the assembly path — an equally spaced [`GeneralMesh`](@ref) +deliberately does not take the [`CirculantMass`](@ref) route that a [`UniformMesh`](@ref) with +the same breakpoints does. + +```@example mesh +u = UniformMesh(4, 0 .. 1) +g = GeneralMesh(collect(range(0, 1; length = 5))) +u == g, u ≈ g, isequal(u, g) +``` + +So the two are interchangeable as values and distinct as dictionary keys, and an equally +spaced `GeneralMesh` is a legitimate way to force the general assembly path in a test. + +## Keyword arguments worth knowing + +`GradedMesh`'s `amplitude` is the ``a`` in the map ``s \mapsto s + (a/2\pi)\sin 2\pi s``. It +must satisfy ``|a| < 1`` for the map to stay increasing. The derivative of the map varies +between ``1-a`` and ``1+a``, so the ratio of the widest to the narrowest cell approaches +``(1+a)/(1-a)`` as the mesh is refined — the two below differ because a cell width is a finite +difference of the map, not its derivative. + +```@example mesh +w = diff(breakpoints(GradedMesh(12, 0 .. 1; amplitude = 0.4))) +maximum(w) / minimum(w), (1 + 0.4) / (1 - 0.4) +``` + +`RandomMesh`'s widths are ``1 + \sigma r_i`` with ``r_i`` uniform on ``[0,1)`` and `spread` +``= \sigma``, rescaled to fill the domain. The stream is seeded, so the mesh is reproducible; +a different `seed` gives a different mesh. + +```@example mesh +breakpoints(RandomMesh(4, 0 .. 1; seed = 7)) == + breakpoints(RandomMesh(4, 0 .. 1; seed = 7)), +breakpoints(RandomMesh(4, 0 .. 1; seed = 7)) == + breakpoints(RandomMesh(4, 0 .. 1; seed = 8)) +``` + +## What is rejected + +```@example mesh +for bad in (() -> UniformMesh(0, 0 .. 1), + () -> UniformMesh(4, 1 .. 0), + () -> UniformMesh(4, 0 .. Inf), + () -> GradedMesh(4, 0 .. 1; amplitude = 1.5), + () -> RandomMesh(4, 0 .. 1; spread = -1), + () -> GeneralMesh([0.0]), + () -> GeneralMesh([0.0, 0.5, 0.5, 1.0]), + () -> GeneralMesh([0.0, 0.7, 0.3])) + try + bad() + catch err + println(err.msg) + end +end +``` + +A repeated breakpoint is rejected rather than tolerated: it is a cell of zero width, which the +Cox-de Boor recursion would silently skip, leaving a basis whose dimension does not match its +knot vector. Multiplicity belongs in the *knot vector*, where the basis puts it, not in the +mesh. + +## Finding a cell + +[`findcell`](@ref) returns the index of the cell containing a point, in `1:ncells`. The cells +are half-open except for the last, which is closed, so the right endpoint belongs to cell `n`. +A point outside the domain is clamped to the nearest cell. + +```@example mesh +u8 = UniformMesh(8, 0 .. 1) +findcell.(Ref(u8), [0.0, 0.06, 0.125, 0.99, 1.0, 5.0, -3.0]) +``` + +On a [`UniformMesh`](@ref) this is a division and needs no breakpoint vector at all; on the +others it is a `searchsortedlast`. Both are allocation-free when called on a *basis*, which +caches the breakpoints — see [Bases](@ref usage-bases). diff --git a/docs/src/usage/quadrature.md b/docs/src/usage/quadrature.md new file mode 100644 index 0000000..1ffbf52 --- /dev/null +++ b/docs/src/usage/quadrature.md @@ -0,0 +1,335 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Assembly](@id usage-assembly) + +Everything this package assembles comes out of one data structure. A +[`SplineQuadrature`](@ref) tabulates a basis and its derivatives at the Gauß-Legendre points +of every cell, and every matrix of the form + +```math +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, + +```math +A = \Phi_a \, \mathrm{diag}(f \odot w) \, \Phi_b^{\mathsf T} . +``` + +```@example asm +using SimpleSplines +using LinearAlgebra +using SparseArrays +``` + +## Construction + +```julia +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`](@ref)) and which representation the mass matrix takes +([`mass_operator`](@ref)); the contractions are identical. + +```@example asm +b = BSplineBasis(UniformMesh(16, 0 .. 2π), 3, Periodic()) +q = SplineQuadrature(b) +nbasis(q), degree(q), ncells(q), length(quadrature_nodes(q)) +``` + + - **`nq`** is the number of quadrature points per cell. The default, + [`quadrature_order`](@ref)`(p) = ⌈3p/2⌉`, integrates degree ``3p-1`` exactly — see + [Choosing `nq`](@ref 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`](@ref) for more is an error + naming the fix. + +```@example asm +quadrature_order.(0:5) +``` + +```@example asm +try + basis_values(SplineQuadrature(b; dmax = 1), 2) +catch err + println(err.msg) +end +``` + +## The table + +| call | returns | +|:--|:--| +| [`quadrature_nodes`](@ref)`(q)` | the `n*nq` global points, in cell order | +| [`quadrature_weights`](@ref)`(q)` | the matching weights, scaled by the cell widths, so `sum(w) == L` | +| [`basis_values`](@ref)`(q, d = 0)` | the sparse `N × n*nq` table ``\Phi_d[i,r] = D^d\phi_i(x_r)`` | +| [`basis`](@ref)`(q)` | the basis | +| [`mass_operator`](@ref)`(q)` | the [`MassOperator`](@ref) | +| `nbasis`, `degree`, `order`, `ncells`, `domainlength`, `meshwidth`, `eltype` | forwarded to the basis | + +```@example asm +sum(quadrature_weights(q)) ≈ domainlength(b) +``` + +The global index of the `r`-th point of cell `k` is `(k-1)*nq + r`, so a hand-written assembly +loop reads + +```@example asm +Φ = 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 +``` + +### 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. + +```@example asm +nnz(Φ), ncells(q) * (degree(q) + 1) * nq, nnz(Φ) / length(Φ) +``` + +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 + +| call | matrix | +|:--|:--| +| [`mass_matrix`](@ref)`(q)` | ``\mathbb{M}_{kl} = \int \phi_k \phi_l`` | +| [`stiffness_matrix`](@ref)`(q)` | ``\mathbb{K}_{kl} = \int \phi_k' \phi_l'`` | +| [`derivative_matrix`](@ref)`(q)` | ``S_{kl} = \int \phi_k \phi_l'`` | +| [`mixed_matrix`](@ref)`(q, a, b)` | ``\int D^a\phi_k \, D^b\phi_l`` | +| [`weighted_matrix`](@ref)`(q, f, a, b)` | ``\int f \, D^a\phi_k \, D^b\phi_l`` | +| [`basis_integrals`](@ref)`(q)` | the vector ``\int \phi_i`` | + +`stiffness_matrix` and `derivative_matrix` are `mixed_matrix(q, 1, 1)` and +`mixed_matrix(q, 0, 1)` — the same objects, not merely equal ones: + +```@example asm +stiffness_matrix(q) === mixed_matrix(q, 1, 1), +derivative_matrix(q) === mixed_matrix(q, 0, 1) +``` + +!!! note "`mixed_matrix` is memoised, and returns the cached object" + These are constants of the discretisation and are asked for inside every Newton iteration + of every step of a downstream integrator, so each `(a, b)` is assembled once and kept. + The consequence is that the result is **shared**: mutating it changes what every later + call returns. Take a `copy` before writing into one. + +```@example asm +mixed_matrix(q, 0, 2) === mixed_matrix(q, 0, 2) +``` + +The same holds for [`mass_matrix`](@ref) and [`basis_integrals`](@ref), 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: + +```@example asm +S = derivative_matrix(q) +maximum(abs, S + S') +``` + +The stiffness matrix is symmetric positive semi-definite with the constants in its kernel: + +```@example asm +K = stiffness_matrix(q) +(maximum(abs, K - K'), + minimum(eigvals(Symmetric(Matrix(K)))), + maximum(abs, K * ones(nbasis(q)))) +``` + +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: + +```@example asm +maximum(abs, mixed_matrix(q, 1, 2) + mixed_matrix(q, 0, 3)) +``` + +### A variable coefficient costs no more than a constant one + +[`weighted_matrix`](@ref) takes either a function of the coordinate or a vector already +sampled at [`quadrature_nodes`](@ref) — the second form is what a coefficient that is itself a +spline expansion becomes, and it avoids resampling a field already in hand. + +```@example asm +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)) +``` + +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`](@ref), which carries no +weight and is always `eltype(q)`: + +```@example asm +eltype(weighted_matrix(q, cis, 0, 0)), eltype(mixed_matrix(q, 0, 0)) +``` + +!!! warning "`weighted_matrix` allocates on every call, and there is no in-place form" + It is the one assembly that cannot be memoised — it depends on the field, so a time + integrator asks for a *different* one inside every Newton iteration. Measured at + ``N = 128``, ``p = 3``, ``n_q = 5`` that is about 260 kB per call, most of it the freshly + built `SparseMatrixCSC` result rather than a temporary. `weighted_matrix!` does not exist + yet; a caller in that position should hold the matrix across steps and accept the + allocation, or assemble the contraction by hand against a cached sparsity pattern. + +## Projection + +```julia +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. + +```@example asm +û = l2_projection(q, sin) +maximum(abs(evaluate(b, û, x) - sin(x)) for x in range(0, 2π; length = 401)) +``` + +```@example asm +v̂ = similar(û) +l2_projection!(v̂, q, sin) +v̂ == û +``` + +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`](@ref) and the [`BandedMass`](@ref) solve are themselves +allocation-free. The exception is a periodic basis on a non-uniform mesh — the +[`FactorizedMass`](@ref) case — where CHOLMOD has no in-place solve. + +```@example asm +fq = sin.(quadrature_nodes(q)) +l2_projection!(v̂, q, fq) # warm up +@allocated l2_projection!(v̂, q, fq) +``` + +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`](@ref), 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. + +| basis | mesh | representation | a solve is | +|:--|:--|:--|:--| +| [`PeriodicBSplineBasis`](@ref) | [`UniformMesh`](@ref) | [`CirculantMass`](@ref) | two planned transforms and a pointwise division, ``O(N\log N)``, no allocation | +| [`PeriodicBSplineBasis`](@ref) | any other | [`FactorizedMass`](@ref) | a sparse Cholesky; the only representation whose solve allocates | +| bounded — clamped or recombined | any | [`BandedMass`](@ref) | a banded Cholesky, ``O(Np)``, no allocation | + +```@example asm +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 +``` + +The circulant construction *verifies* circulance rather than assuming it, to `atol = 1e-10` by +default. 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. + +```@example asm +qg = SplineQuadrature(BSplineBasis(GradedMesh(16, 0 .. 2π), 3, Periodic())) +try + CirculantMass(mass_matrix(qg), nbasis(qg)) +catch err + println(err.msg) +end +``` + +### Solving + +```julia +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) +``` + +```@example asm +op = mass_operator(q) +rhs = basis_integrals(q) +x = op \ rhs +maximum(abs, x .- 1) # M⁻¹ ∫φᵢ = 1, by the partition of unity +``` + +```@example asm +y = similar(rhs) +mass_solve!(y, op, rhs) +y == x, (@allocated mass_solve!(y, op, rhs)) +``` + +[`mass_factorization`](@ref)`(q)` is the same operator, named for the case where the +factorisation rather than the matrix is what is wanted. + +## [Choosing `nq`](@id choosing-nq) + +An `nq`-point Gauß-Legendre rule is exact to degree ``2n_q - 1``. Three thresholds matter: + +| what needs to be exact | degree | requires | +|:--|:--|:--| +| the mass matrix ``\int \phi_k \phi_l`` | ``2p`` | ``n_q \ge p + 1`` | +| the antisymmetry of [`derivative_matrix`](@ref), i.e. ``\int (\phi_k\phi_l)'`` | ``2p - 1`` | ``n_q \ge p`` | +| a quadratic nonlinearity, ``\int u_h u_{h,x} \phi_i`` | ``3p - 1`` | ``n_q \ge \lceil 3p/2 \rceil``, the default | + +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: + +```@example asm +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 +``` + +Too coarse an `nq` makes the mass matrix **singular** rather than merely inexact, and the +constructor says so by name: + +```@example asm +try + SplineQuadrature(BSplineBasis(UniformMesh(8, 0 .. 1), 3); nq = 1) +catch err + println(err.msg) +end +``` + +## Two traps + +!!! warning "One quadrature per thread" + A `SplineQuadrature` carries mutable state behind an otherwise read-only interface: + [`mixed_matrix`](@ref) memoises into a `Dict`, and [`l2_projection!`](@ref) forms its + ``f \odot w`` product in a shared buffer. Neither is synchronised, so **one quadrature must + not be used from two threads at once.** Give each thread its own, or stay with the + allocating [`l2_projection`](@ref), which forms its own product. + +!!! warning "Several results are returned by reference" + [`mass_matrix`](@ref), [`basis_integrals`](@ref), [`basis_values`](@ref) and every + [`mixed_matrix`](@ref) hand back the quadrature's own storage. `copy` before mutating. diff --git a/docs/src/usage/splines.md b/docs/src/usage/splines.md new file mode 100644 index 0000000..b347ad9 --- /dev/null +++ b/docs/src/usage/splines.md @@ -0,0 +1,178 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Splines](@id usage-splines) + +A basis answers "what is ``\phi_j`` here". A [`Spline`](@ref) answers "what is ``u_h`` here": +it is a basis together with the coefficients of one element of its span, and it is callable. + +```@example spl +using SimpleSplines +``` + +## Construction + +```julia +Spline(basis, coefficients) # coefficients: a vector, or a D-array for a product basis +Spline(basis) # the zero spline, freshly allocated in the right shape +``` + +```@example spl +b = BSplineBasis(UniformMesh(12, 0 .. 1), 3) +q = SplineQuadrature(b) +s = Spline(b, l2_projection(q, x -> sinpi(2x))) +s, s(0.125), s(0.125, 1) +``` + +```@example spl +z = Spline(b) +size(z), all(iszero, coefficients(z)) +``` + +The shape is checked: a one-dimensional basis needs a vector of `nbasis(basis)` entries, a +product basis an array of `size(basis)`. + +```@example spl +try + Spline(b, zeros(nbasis(b) - 1)) +catch err + println(err.msg) +end +``` + +## Evaluation + +```julia +s(x) s(x, d) +evaluate(s, x) evaluate(s, x, d) +``` + +`d` is the derivative order — an integer on one axis, a per-axis tuple on a product basis. + +```@example spl +xs = range(0, 1; length = 5) +s.(xs), s.(xs, 1) +``` + +Evaluating a spline is ``O(p^2)`` and not ``O(N)``: only the local block of `p+1` functions +contributes, and the value is that block's weighted sum. Summing `evaluate(b, j, x)` over the +whole index range would give the same answer and cost ``O(N)``. + +```@example spl +s(0.3) ≈ sum(coefficients(s)[j] * b[0.3, j] for j in eachindex(b)) +``` + +Outside a bounded domain the value is zero, since every basis function vanishes there. + +```@example spl +s(-0.1), s(1.4) +``` + +## Derivatives as objects + +[`derivative`](@ref)`(s, d = 1)` gives a [`SplineDerivative`](@ref), which is callable and so +broadcasts the way `s` does. It holds no state of its own — it shares the coefficient array — +so a reprojection is visible through it with nothing rebuilt. + +```@example spl +ds = derivative(s) +ds, ds.(xs), maximum(abs(ds(x) - 2π * cospi(2x)) for x in range(0, 1; length = 201)) +``` + +```@example spl +derivative(s, 2)(0.3) ≈ s(0.3, 2) +``` + +This is the shape a right-hand side wants when the value and the derivative appear together in +one broadcast: `derivative(fs).(v) ./ fs.(v)` reads as it should. + +```@example spl +maximum(abs, derivative(s).(xs) ./ (1 .+ s.(xs) .^ 2)) +``` + +## Accessors + +```@example spl +(basis(s) === b, eltype(s), ndims(s), size(s), length(s), + nbasis(s), degree(s), order(s), domain(s)) +``` + +| call | returns | +|:--|:--| +| [`basis`](@ref)`(s)` | the basis | +| [`coefficients`](@ref)`(s)` | the coefficient array **itself**, not a copy | +| `eltype`, `ndims`, `size`, `length` | of the coefficient array | +| `nbasis`, `degree`, `order`, `domain` | forwarded to the basis | +| `similar(s)`, `similar(s, T)`, `copy(s)` | a spline on the same basis with fresh coefficients | + +!!! note "`coefficients` aliases, and that is the point" + It hands back the array the spline holds, so writing into it changes the spline. That is + what lets a projection be written into a spline already wired into a right-hand side, and + what lets a particle deposition rewrite the coefficients every step while the basis, the + quadrature and the mass factorisation stay put. + +```@example spl +t = Spline(b, copy(coefficients(s))) +coefficients(t) .= 0 +t(0.3) +``` + +For the same reason, a `Spline` built on an array you already hold shares it: + +```@example spl +û = zeros(nbasis(b)) +u = Spline(b, û) +û[5] = 1.0 +u(0.3) == b[0.3, 5] +``` + +Use `copy(s)` when you want an independent spline, and `similar(s)` when you want the same +shape with undefined contents. + +## Projecting into an existing spline + +[`l2_projection!`](@ref)`(s, q, f)` writes into the array `s` already holds and returns `s`, +so every reference to `s` sees the new function. + +```@example spl +l2_projection!(u, q, cospi) +u(0.25), coefficients(u) === û +``` + +## On a product basis + +Everything above carries over, with a coefficient array and tuple-valued derivative orders. +The integer form `derivative(s, k)` is the shorthand for ``\partial_k``. + +```@example spl +B = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(8, 0 .. 2), 3) +Q = TensorProductQuadrature(B) +g(x) = x[1]^2 * x[2] +S = Spline(B, l2_projection(Q, g)) +S, ndims(S), size(S) +``` + +```@example spl +(S((0.3, 1.4)), + S((0.3, 1.4), (1, 0)), + derivative(S, 1)((0.3, 1.4)), + derivative(S, 2)((0.3, 1.4))) +``` + +``\partial_1 g = 2x_1x_2`` and ``\partial_2 g = x_1^2``, so at ``(0.3, 1.4)`` those are +``0.84`` and ``0.09``. + +```@example spl +try + derivative(S, 3) +catch err + println(err.msg) +end +``` + +A mixed derivative is given directly as a tuple: + +```@example spl +S((0.3, 1.4), (1, 1)) # ∂₁∂₂ g = 2x₁ +``` diff --git a/docs/src/usage/tensorproduct.md b/docs/src/usage/tensorproduct.md new file mode 100644 index 0000000..2146089 --- /dev/null +++ b/docs/src/usage/tensorproduct.md @@ -0,0 +1,315 @@ +```@meta +CurrentModule = SimpleSplines +``` + +# [Tensor Products](@id usage-tensorproduct) + +A [`TensorProductBasis`](@ref) is ``D`` one-dimensional bases, one per axis, in any number of +dimensions. For the mathematics see [Tensor Products](@ref theory-tensorproduct). + +```@example utp +using SimpleSplines +using LinearAlgebra +using Random +Random.seed!(1234) +``` + +## Construction + +```julia +TensorProductBasis(b₁, b₂, …) TensorProductBasis((b₁, b₂, …)) +b₁ ⊗ b₂ ⊗ b₃ +``` + +`⊗` associates **flat**: three factors give one three-factor basis, not a nest of two. + +```@example utp +bx = BSplineBasis(UniformMesh(10, 0 .. 2π), 3, Periodic()) +bv = BSplineBasis(UniformMesh(8, -10 .. 10), 4) +bz = BSplineBasis(UniformMesh(6, 0 .. 1), 2, Dirichlet()) + +B = bx ⊗ bv ⊗ bz +ndims(B), size(B), nbasis(B) +``` + +```@example utp +bases(B) === (bx, bv, bz), (bx ⊗ bv) ⊗ bz == bx ⊗ (bv ⊗ bz) +``` + +A one-factor product is allowed, and is occasionally the point — it lets code written for `D` +dimensions run unchanged at ``D = 1``. Zero factors are rejected. + +```@example utp +B1 = TensorProductBasis(BSplineBasis(UniformMesh(7, 0 .. 1), 3)) +ndims(B1), size(B1) +``` + +## Accessors return per-axis tuples + +```@example utp +degree(B), order(B), ncells(B), size(B), local_width(B) +``` + +| call | returns | +|:--|:--| +| `ndims(B)`, `size(B)`, `size(B, d)`, `length(B)` | ``D``, the per-axis dimensions, one of them, and ``\prod_d N_d`` | +| [`nbasis`](@ref)`(B)` | `length(B)` — the total, **not** a tuple | +| [`degree`](@ref), [`order`](@ref), [`ncells`](@ref), [`mesh`](@ref), [`boundary`](@ref), [`local_width`](@ref) | per-axis tuples | +| [`meshwidth`](@ref)`(B)` | the **maximum** over the axes, a scalar | +| [`domain`](@ref)`(B)` | a `DomainSets` `ProductDomain`, so `x ∈ domain(B)` answers for a ``D``-vector | +| [`nodes`](@ref)`(B)` | the per-axis node vectors; the grid is their product and is not materialised | +| [`polynomial_reproduction`](@ref)`(B)` | the **minimum** over the axes | +| `CartesianIndices(B)`, `LinearIndices(B)`, `eachindex(B)` | index sets; `eachindex` is the Cartesian one | + +```@example utp +[1.0, 2.0, 0.5] ∈ domain(B), meshwidth(B), polynomial_reproduction(B) +``` + +Two of those are easy to misread. `nbasis(B)` is the *total* dimension, because that is the +length a flattened coefficient vector has; `size(B)` is the per-axis tuple. And +`boundary(B)` is a tuple of whatever each axis reports, so a periodic axis contributes a bare +`Periodic()` while a bounded one contributes a pair: + +```@example utp +boundary(B) +``` + +## Coefficients + +A coefficient array has size `size(B)`, one axis per factor. + +```@example utp +û = zeros(size(B)...) +û[3, 4, 2] = 1.0 +size(û) == size(B) +``` + +When it has to be flattened, the convention is that the **first axis varies fastest** — `vec`, +`LinearIndices(B)`, Julia's column-major order — and the assembled mass matrix is +`kron(M_D, …, M_1)` to match. + +```@example utp +LinearIndices(B)[3, 4, 2], findfirst(!iszero, vec(û)) +``` + +## Evaluating + +```julia +evaluate(B, I, x, d = ntuple(_ -> 0, D)) # I a CartesianIndex, a tuple, or a linear index +evaluate(B, û, x, d = ntuple(_ -> 0, D)) # the spline +B[x, I] B(x, I) +``` + +`x` is any indexable ``D``-vector: a tuple, an `SVector`, a `Vector`. As on one axis, +[`evaluate`](@ref) takes the **index first** and `getindex` the **point first**. + +```@example utp +Bp = BSplineBasis(UniformMesh(8, 0 .. 1), 3) ⊗ BSplineBasis(UniformMesh(8, 0 .. 1), 2) +b1, b2 = bases(Bp) +pt = (0.3, 0.4) +(evaluate(Bp, (3, 4), pt), + evaluate(Bp, CartesianIndex(3, 4), pt), + evaluate(Bp, LinearIndices(Bp)[3, 4], pt), + Bp[pt, (3, 4)], + evaluate(b1, 3, 0.3) * evaluate(b2, 4, 0.4)) +``` + +### Derivatives are a per-axis tuple + +There is no scalar derivative order. `d` is an `NTuple{D,Int}` of orders, giving +``\partial_1^{d_1}\cdots\partial_D^{d_D}``; the gradient component ``\partial_k`` is +`ntuple(i -> i == k ? 1 : 0, D)`. + +```@example utp +evaluate(Bp, (3, 4), pt, (1, 2)) ≈ + evaluate(b1, 3, 0.3, 1) * evaluate(b2, 4, 0.4, 2) +``` + +For a [`Spline`](@ref) on a product basis, `derivative(s, k::Integer)` builds that tuple for +you — see [Splines](@ref usage-splines). + +### The local path + +[`evaluate_all!`](@ref) takes a **tuple of buffers**, one per axis, each of length +`local_width(bases(B)[k])`, and returns the tuple of first indices. The ``D``-dimensional block +is deliberately not materialised: at ``D = 3`` with cubics that is 64 numbers per particle, +and the loop that consumes them can form each product as it goes. + +```@example utp +bufs = ntuple(k -> zeros(local_width(bases(Bp)[k])), ndims(Bp)) +j₀ = evaluate_all!(bufs, Bp, pt) +j₀, map(length, bufs) +``` + +```@example utp +all(bufs[k] ≈ evaluate_all(bases(Bp)[k], pt[k])[2] for k in 1:ndims(Bp)) +``` + +A deposition then reads: + +```@example utp +sz = size(Bp) +coeffs = zeros(sz) +for (pos, wt) in ((0.13, 0.22) => 1.0, (0.71, 0.44) => 2.5) + jj = evaluate_all!(bufs, Bp, pos) + for t in CartesianIndices(map(length, bufs)) + I = ntuple(k -> basis_index(bases(Bp)[k], jj[k] + t[k] - 1), ndims(Bp)) + all(k -> 1 ≤ I[k] ≤ sz[k], 1:ndims(Bp)) || continue + coeffs[I...] += wt * prod(ntuple(k -> bufs[k][t[k]], ndims(Bp))) + end +end +sum(coeffs) +``` + +The sum is ``3.5``, each particle having deposited its whole weight. As on one axis, the +returned indices are the ones *before* wrapping, so they go through [`basis_index`](@ref); and +an index outside `1:N_d` on a bounded axis is genuinely absent rather than zero, so it is +skipped. + +A point outside the box evaluates to exactly zero: + +```@example utp +BD = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet()) ⊗ + BSplineBasis(UniformMesh(8, 0 .. 1), 3) +v̂ = randn(size(BD)...) +evaluate(BD, v̂, (1.5, 0.5)), evaluate(BD, v̂, (0.5, -0.2)) +``` + +## Assembly + +```julia +TensorProductQuadrature(B; nq = map(quadrature_order, degree(B)), dmax = 3) +``` + +`nq` and `dmax` may each be a ``D``-tuple or a scalar broadcast to every axis. There is no +``D``-dimensional tabulation: the quadrature is ``D`` one-dimensional +[`SplineQuadrature`](@ref)s plus the [`KroneckerMass`](@ref) built from their mass operators. + +```@example utp +Q = TensorProductQuadrature(Bp; nq = (5, 7), dmax = (2, 3)) +map(length, quadrature_nodes(Q)), quadrature_grid_size(Q) +``` + +| call | returns | +|:--|:--| +| [`quadratures`](@ref)`(Q)` | the tuple of one-dimensional [`SplineQuadrature`](@ref)s | +| [`basis`](@ref)`(Q)` | the [`TensorProductBasis`](@ref) | +| [`quadrature_nodes`](@ref)`(Q)`, [`quadrature_weights`](@ref)`(Q)` | per-axis **tuples** of vectors | +| [`quadrature_grid_size`](@ref)`(Q)` | `(n_d * nq_d)` per axis | +| [`quadrature_sample`](@ref)`(Q, f)` | `f` evaluated on the grid, as an array of that size | +| [`mass_operator`](@ref)`(Q)`, [`mass_matrix`](@ref)`(Q)` | the [`KroneckerMass`](@ref), and its assembled Kronecker product | +| `ndims(Q)`, `size(Q)`, `nbasis(Q)`, `degree(Q)`, `eltype(Q)` | forwarded to the basis | + +The per-axis structure is why the one-dimensional accessors are reached through +`quadratures(Q)[k]`: + +```@example utp +q1, q2 = quadratures(Q) +size(basis_values(q1, 0)), size(basis_values(q2, 3)) +``` + +```@example utp +try + basis_values(q1, 3) # axis 1 was built with dmax = 2 +catch err + println(err.msg) +end +``` + +Note that a grid point is `(x[1][r₁], x[2][r₂], …)` and its weight the *product* of the +per-axis weights — [`quadrature_sample`](@ref) and [`contract`](@ref) apply that for you. + +```@example utp +X = quadrature_nodes(Q) +F = quadrature_sample(Q, x -> x[1]^2 * x[2]) +size(F) == quadrature_grid_size(Q), F[3, 4] ≈ X[1][3]^2 * X[2][4] +``` + +### Contraction and projection + +[`contract`](@ref)`(Q, F, d)` turns an array given on the quadrature grid into the load array + +```math +L_{i_1 \dots i_D} = \sum_{q_1 \dots q_D} F_{q_1 \dots q_D} + \prod_{k} D^{d_k} \phi^{(k)}_{i_k}(x_{q_k}) \, w_{q_k} , +``` + +as ``D`` sparse matrix products. `F` must already contain the whole integrand; the weights are +applied by `contract` itself, and always on a *copy*, so the caller's sample is not scaled. + +```@example utp +L = contract(Q, F) +size(L) == size(Bp), F[3, 4] ≈ X[1][3]^2 * X[2][4] # F is unchanged +``` + +[`l2_projection`](@ref) is then that contraction followed by the factored mass solve. `f` may +be a function of a ``D``-tuple or an array already sampled on the grid. + +```@example utp +f(x) = x[1]^2 * x[2] + 3x[1] - x[2]^2 + 1 +û2 = l2_projection(Q, f) +ŵ = similar(û2) +l2_projection!(ŵ, Q, f) +(size(û2), + maximum(abs, ŵ - û2), + maximum(abs, l2_projection(Q, quadrature_sample(Q, f)) - û2)) +``` + +```@example utp +maximum(abs(evaluate(Bp, û2, (x, y)) - f((x, y))) + for x in range(0, 1; length = 21), y in range(0, 1; length = 21)) +``` + +That is round-off rather than an approximation: `Bp` is cubic × quadratic, so +[`polynomial_reproduction`](@ref) is `2` and ``x^2 y`` is reproduced exactly. + +## The Kronecker mass operator + +[`KroneckerMass`](@ref) holds the ``D`` one-dimensional operators and is never assembled. Each +factor keeps whatever representation its own axis earned: + +```@example utp +op = mass_operator(TensorProductQuadrature(bx ⊗ bv)) +map(fac -> nameof(typeof(fac)), mass_factors(op)), size(op), ndims(op) +``` + +```julia +op * X op \ X ldiv!(Y, op, X) ldiv!(op, X) mass_solve!(Y, op, X) +mass_matrix(op) Matrix(op) mass_factors(op) +``` + +Arrays of size `size(B)` **and** flat vectors of length `length(B)` are both accepted; a +vector is reshaped on the first-axis-fastest convention. + +```@example utp +opp = mass_operator(Q) +Û = randn(size(Bp)...) +(maximum(abs, opp \ vec(Û) - vec(opp \ Û)), + maximum(abs, opp \ (opp * Û) - Û)) +``` + +[`mass_matrix`](@ref)`(op)` forms `kron(M_D, …, M_1)` on demand and does **not** store it — it +is what the representation exists to avoid, and at three dimensions it will not fit. It is +provided so that a test can check the factored solve against the dense one at a size where +both are possible. + +```@example utp +M1, M2 = map(mass_matrix, mass_factors(opp)) +maximum(abs, Matrix(opp) - kron(Matrix(M2), Matrix(M1))) +``` + +## What is not provided + +A matrix ``\int_\Omega f \, D^a\Phi_I \, D^b\Phi_J`` with a **non-separable** ``f`` does not +factorise, and there is no `weighted_matrix` for a tensor-product quadrature. Building one +means forming the Kronecker product of the one-dimensional tabulations explicitly and paying +the storage — at ``64^2`` cubic cells that is some 26 MB per derivative multi-index. The +pieces are exposed if you need it: + +```@example utp +Φ1, Φ2 = basis_values(quadratures(Q)[1], 0), basis_values(quadratures(Q)[2], 0) +size(Φ1), size(Φ2), size(kron(Φ2, Φ1)) +``` + +Note the factor order, `kron(Φ2, Φ1)`: the same first-axis-fastest convention as everywhere +else. diff --git a/scripts/evaluate_cost_scaling.jl b/scripts/evaluate_cost_scaling.jl new file mode 100644 index 0000000..cf47158 --- /dev/null +++ b/scripts/evaluate_cost_scaling.jl @@ -0,0 +1,84 @@ +# How the cost of `evaluate` and of `evaluate_all` grows with the degree, and therefore which +# complexity their docstrings should state. +# +# `_bspline` is the Cox-de Boor recursion written out as it stands, with no memoisation, so it +# splits into two subproblems at every level and the cost of one value is `O(2^p)`. +# `evaluate_all!` runs de Boor's triangular scheme instead, which is `O(p^2)` for all `p+1` +# nonzero functions *together* -- so the block form is not merely `p+1` times cheaper, it has a +# different exponent. At the cubics almost everything here uses, the difference hardly shows. +# +# This script is what the figures quoted in `evaluate`'s and `_bspline`'s docstrings, and in +# the cost table of `docs/src/theory/bsplines.md`, are read off. Absolute timings are a +# property of the machine, so the assertion at the end is on the *shape*: doubling the degree +# has to cost the single-function path far more than it costs the block. +# +# Run: julia --project=. scripts/evaluate_cost_scaling.jl + +using SimpleSplines +using Printf + +const DEGREES = (3, 6, 8, 12) +const CALLS = 20_000 +const REPEATS = 5 +const NCELLS = 40 +const X = 0.317 # interior, and not a breakpoint of the mesh below + +# The minimum over `REPEATS` runs, which is the estimator least polluted by whatever else the +# machine is doing. The value is accumulated so that the loop cannot be optimised away, and a +# zero total means the chosen function does not reach `X` and the timing is of nothing. +function nanoseconds_per_call(f) + f() # compile + best = Inf + for _ in 1:REPEATS + acc = 0.0 + t = @elapsed for _ in 1:CALLS + acc += f() + end + iszero(acc) && error("the accumulated value is zero; nothing was measured") + best = min(best, t) + end + return best / CALLS * 1e9 +end + +function measure(p) + b = BSplineBasis(UniformMesh(NCELLS, 0 .. 1), p) + # Any function nonzero at X. On a clamped basis those on cell c are c .. c+p. + j = first(local_indices(b, findcell(b, X))) + one_function = nanoseconds_per_call(() -> evaluate(b, j, X)) + # `evaluate_all`, not `evaluate_all!`: it is the form the docstrings compare against, and + # its freshly allocated buffer is part of what a caller pays. + whole_block = nanoseconds_per_call(() -> sum(last(evaluate_all(b, X)))) + return (; p, one_function, whole_block) +end + +rows = [measure(p) for p in DEGREES] + +@printf("%3s %14s %14s %8s\n", "p", "evaluate [ns]", "block [ns]", "ratio") +for r in rows + @printf("%3d %14.1f %14.1f %8.1f\n", + r.p, r.one_function, r.whole_block, r.one_function / r.whole_block) +end + +# The claim under test. Between p = 6 and p = 12 the degree doubles, so an `O(2^p)` path pays +# a factor of about 2^6 = 64 while an `O(p^2)` one pays about 4. The bounds are loose enough +# to survive a noisy machine and tight enough that the two cannot be confused: an `O(p^2)` +# single-function path would come nowhere near 20x. +lo = only(filter(r -> r.p == 6, rows)) +hi = only(filter(r -> r.p == 12, rows)) + +growth_one = hi.one_function / lo.one_function +growth_block = hi.whole_block / lo.whole_block + +@printf("\np = 6 -> 12: evaluate x%.1f block x%.1f\n", growth_one, growth_block) + +if growth_one < 20 + error("`evaluate` grew only $(round(growth_one; digits = 1))x from p = 6 to p = 12; " * + "the O(2^p) the docstrings claim would grow about 64x") +end +if growth_block > 8 + error("`evaluate_all` grew $(round(growth_block; digits = 1))x from p = 6 to p = 12; " * + "the O(p^2) the docstrings claim would grow about 4x") +end + +println("\nOK: the single-function path grows exponentially in the degree and the block " * + "path polynomially, as the docstrings state.") diff --git a/src/bspline.jl b/src/bspline.jl index 1ba1958..78eb458 100644 --- a/src/bspline.jl +++ b/src/bspline.jl @@ -34,6 +34,16 @@ built on them doubles as a check of the formulae themselves. It is the *referenc implementation: [`evaluate_all`](@ref) computes the same numbers in ``O(p^2)`` for the whole local block, and the test suite checks the two against each other. +Written as it stands the recursion is *unmemoised*, so it splits into two subproblems at every +level and the cost of one value is ``O(2^p)`` rather than ``O(p^2)``. That is deliberate — the +point of this function is to be the formula, not to be fast — but it is the reason +[`evaluate_all!`](@ref) exists rather than a convenience wrapper around a loop over this one. +Measured at 26 ns per call at ``p = 3`` against 21 ns for the *whole block* through +[`evaluate_all`](@ref), and 15.5 µs against 104 ns at ``p = 12``; doubling the degree from +``6`` to ``12`` costs the single value 64 times more and the block 2.6 times more, which is +``2^6`` against ``(12/6)^2``. `scripts/evaluate_cost_scaling.jl` measures it and fails if the +two exponents stop being distinguishable. + A knot span of zero length — a repeated knot — contributes nothing, which is what makes the recursion well defined at a knot of multiplicity greater than one. The guard is on the span being positive rather than on the numerator, because it is the division that fails. @@ -167,8 +177,22 @@ Base.eltype(::AbstractBSplineBasis{T}) where {T} = T Base.eltype(::Type{<:AbstractBSplineBasis{T}}) where {T} = T Base.eachindex(b::AbstractBSplineBasis) = Base.OneTo(nbasis(b)) Base.axes(b::AbstractBSplineBasis) = (Inclusion(domain(b)), eachindex(b)) +""" + grid(b::AbstractBSplineBasis) + +The [`nodes`](@ref) 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. +""" ContinuumArrays.grid(b::AbstractBSplineBasis) = nodes(b) +""" + nnodes(b::AbstractBSplineBasis) + +The number of [`nodes`](@ref), which for a spline basis is [`nbasis`](@ref)`(b)`: there is one +Greville abscissa per basis function. + +The generic belongs to `GeometricBase` and is extended here rather than redefined. +""" nnodes(b::AbstractBSplineBasis) = nbasis(b) @doc raw""" @@ -552,10 +576,12 @@ gives zero — which is the mathematically correct value for a compactly support 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. A loop that needs every nonzero -function at a point — a particle deposition, or a matrix assembly — should use -[`evaluate_all`](@ref), which is ``O(p^2)`` for the whole block of `p+1` rather than -``O(p^2)`` for each one. +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`](@ref) 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`. """ function evaluate(b::BSplineBasis{T}, j::Integer, x::Number, d::Integer = 0) where {T} @boundscheck (1 ≤ j ≤ nbasis(b)) || throw(BoundsError(b, j)) @@ -661,6 +687,16 @@ end ## Local evaluation ## --------------------------------------------------------------------------------------- +# The shared implementation. `findcell` itself carries the docstring, below: a `@doc` block +# here would document *this* function instead, leaving the exported name undocumented. +function _findcell(y::AbstractVector, n::Integer, x::Number) + x ≤ y[begin] && return 1 + x ≥ y[end] && return n + # searchsortedlast gives the k with y[k] ≤ x < y[k+1]; both ends are handled above, so + # the result is already in 1:n and needs no clamping. + return searchsortedlast(y, x) +end + @doc raw""" findcell(b::AbstractBSplineBasis, x) findcell(m::Mesh, x) @@ -671,19 +707,15 @@ The right endpoint `b` belongs to the last cell, not to a cell of its own — th 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`](@ref) reduce `x` onto the domain first. -""" -function _findcell(y::AbstractVector, n::Integer, x::Number) - x ≤ y[begin] && return 1 - x ≥ y[end] && return n - # searchsortedlast gives the k with y[k] ≤ x < y[k+1]; both ends are handled above, so - # the result is already in 1:n and needs no clamping. - return searchsortedlast(y, x) -end +It allocates nothing: a basis and a non-uniform mesh both hand over a stored breakpoint +vector, and on a [`UniformMesh`](@ref) the cell index is a division that does not read the +breakpoints at all. +""" findcell(m::Mesh, x::Number) = _findcell(breakpoints(m), ncells(m), x) -# On a uniform mesh the cell index is a division, with no need for the breakpoint vector at -# all. Kept as its own method because it is the one a particle loop actually takes. +# Kept as its own method, rather than left to the generic one above, because it is the method +# a particle loop actually takes. function findcell(m::UniformMesh, x::Number) n = ncells(m) t = (x - m.a) / (m.b - m.a) @@ -692,7 +724,6 @@ function findcell(m::UniformMesh, x::Number) return min(floor(Int, t * n) + 1, n) end -# Reads the breakpoints cached in the basis, so this allocates nothing. findcell(b::AbstractBSplineBasis, x::Number) = _findcell(breakpoints(b), ncells(b), x) findcell(b::BSplineBasis{T, <:UniformMesh}, x::Number) where {T} = findcell(mesh(b), x) function findcell(b::PeriodicBSplineBasis{T, <:UniformMesh}, x::Number) where {T} @@ -703,11 +734,12 @@ end evaluate_all(b::AbstractBSplineBasis, x, d = 0) evaluate_all!(values, b::AbstractBSplineBasis, x, d = 0) -The `d`-th derivatives at `x` of the `p+1` basis functions that do not vanish there, +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 `p+1` entries; `values[t]` is the derivative of basis function +`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)`. ```jldoctest @@ -725,10 +757,11 @@ julia> sum(v) ≈ 1 # the other N - p - 1 functions vanish a 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)`` through [`evaluate`](@ref) one index +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`](@ref) 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. +particle and a loop over the whole basis, and, within each term of it, between the triangular +scheme and the unmemoised recursion. !!! note "The index may need wrapping" `j₀ + t - 1` is the index *before* wrapping, and the block is contiguous only in that @@ -753,7 +786,7 @@ D^{m+1} \phi_j^r = r \left( \frac{D^m \phi_j^{r-1}}{x_{j+r} - x_j} ``` 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. [`_bspline`](@ref) computes the same numbers from the +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. """ @@ -920,5 +953,12 @@ end Base.adjoint(b::AbstractBSplineBasis) = Derivative(axes(b, 1)) * b -# Retained so that code written against the periodic-only version keeps resolving. +""" + PeriodicBSplineDerivative + +The [`BSplineDerivative`](@ref) of a [`PeriodicBSplineBasis`](@ref) specifically. + +A compatibility alias. Dispatch on [`BSplineDerivative`](@ref) instead, which covers all three +bases. +""" const PeriodicBSplineDerivative = QMul2{<:Derivative, <:PeriodicBSplineBasis} diff --git a/src/quadrature.jl b/src/quadrature.jl index e241a51..adda9ec 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -237,6 +237,11 @@ for a bounded basis — clamped or recombined — whose mass matrix has no seam """ mass_operator(q::SplineQuadrature) = q.mass +""" + basis(q::SplineQuadrature) + +The [`AbstractBSplineBasis`](@ref) the quadrature was built for. +""" basis(q::SplineQuadrature) = q.basis nbasis(q::SplineQuadrature) = nbasis(q.basis) degree(q::SplineQuadrature) = degree(q.basis) diff --git a/src/tensorproduct.jl b/src/tensorproduct.jl index d878d7d..b45a642 100644 --- a/src/tensorproduct.jl +++ b/src/tensorproduct.jl @@ -129,6 +129,12 @@ it is not materialised, since at three dimensions it is the whole coefficient ar """ nodes(B::TensorProductBasis) = map(nodes, B.bases) +""" + nnodes(B::TensorProductBasis) + +The number of points of the tensor-product node grid, ``\\prod_d N_d`` — a scalar, unlike +[`nodes`](@ref)`(B)`, which returns the per-axis vectors whose product that grid is. +""" nnodes(B::TensorProductBasis) = length(B) """ @@ -532,9 +538,36 @@ struct TensorProductQuadrature{T, D, BT <: TensorProductBasis{T, D}, QS <: Tuple end end +""" + basis(q::TensorProductQuadrature) + +The [`TensorProductBasis`](@ref) the quadrature was built for. +""" basis(q::TensorProductQuadrature) = q.basis + +""" + quadratures(q::TensorProductQuadrature) + +The tuple of one-dimensional [`SplineQuadrature`](@ref)s, 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. +""" quadratures(q::TensorProductQuadrature) = q.quadratures + +""" + mass_operator(q::TensorProductQuadrature) + +The [`KroneckerMass`](@ref) built from the per-axis mass operators, never assembled. +""" mass_operator(q::TensorProductQuadrature) = q.mass + +""" + mass_matrix(q::TensorProductQuadrature) + +The assembled Kronecker product `kron(M_D, …, M_1)`, formed on demand — see +[`mass_matrix`](@ref)`(op::KroneckerMass)` for why it is not stored. +""" mass_matrix(q::TensorProductQuadrature) = mass_matrix(q.mass) nbasis(q::TensorProductQuadrature) = nbasis(q.basis) degree(q::TensorProductQuadrature) = degree(q.basis)