Commit cbf5586e authored by Juan Ignacio Polanco's avatar Juan Ignacio Polanco

Merge branch 'BSplines'

parents 89256cf1 d0cfb49e
......@@ -4,6 +4,7 @@ authors = ["Juan Ignacio Polanco <jipolanc@gmail.com>"]
version = "0.1.0"
[deps]
BSplines = "488c2830-172b-11e9-1591-253b8a7df96d"
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
......
......@@ -15,7 +15,7 @@ module BasisSplines
export Collocation
export BSplineBasis, Spline, BSpline
export BSplineBasis, Spline, BSpline, BSplineOrder, Derivative
export knots, order, coefficients
export augment_knots
export evaluate_bspline, evaluate_bspline!
......@@ -29,8 +29,25 @@ using LinearAlgebra: Symmetric
using SparseArrays
using StaticArrays: MVector
# We're transitioning to using the registered BSplines package...
import BSplines
using BSplines: Derivative
import BSplines: order, knots
"""
BSplineOrder(k::Integer)
Specifies the B-spline order `k`.
"""
struct BSplineOrder{k} end
BSplineOrder(k::Integer) = BSplineOrder{k}()
include("knots.jl")
include("basis.jl")
const AnyBSplineBasis = Union{BSplineBasis, BSplines.BSplineBasis}
include("spline.jl")
include("galerkin.jl")
......
module Collocation
using ..BasisSplines
import BSplines
using BSplines: Derivative
using BandedMatrices
using SparseArrays
......@@ -114,8 +116,8 @@ end
collocation_matrix(
B::BSplineBasis,
x::AbstractVector,
[deriv::Derivative = Derivative(0)],
[MatrixType = SparseMatrixCSC{Float64}];
Ndiff::Val = Val(0),
clip_threshold = eps(eltype(MatrixType)),
)
......@@ -130,7 +132,7 @@ points:
where `Nx = length(x)` is the number of collocation points, and
`Nb = length(B)` is the number of B-splines in `B`.
To obtain a matrix associated to the B-spline derivatives, set the `Ndiff`
To obtain a matrix associated to the B-spline derivatives, set the `deriv`
argument to the order of the derivative.
Given the B-spline coefficients `{u[j], j = 1:Nb}`, one can recover the values
......@@ -176,15 +178,20 @@ alternatives, especially for large problems.
See also [`collocation_matrix!`](@ref).
"""
function collocation_matrix(B::BSplineBasis, x::AbstractVector,
function collocation_matrix(
B::BSplineBasis, x::AbstractVector,
deriv::Derivative = Derivative(0),
::Type{MatrixType} = SparseMatrixCSC{Float64};
kwargs...) where {MatrixType}
Nx = length(x)
Nb = length(B)
C = allocate_collocation_matrix(MatrixType, (Nx, Nb), order(B))
collocation_matrix!(C, B, x; kwargs...)
collocation_matrix!(C, B, x, deriv; kwargs...)
end
collocation_matrix(B, x, ::Type{M}; kwargs...) where {M} =
collocation_matrix(B, x, Derivative(0), M; kwargs...)
allocate_collocation_matrix(::Type{M}, dims, k) where {M <: AbstractMatrix} =
M(undef, dims...)
......@@ -212,16 +219,18 @@ end
"""
collocation_matrix!(
C::AbstractMatrix{T}, B::BSplineBasis, x::AbstractVector;
Ndiff::Val = Val(0), clip_threshold = eps(T))
C::AbstractMatrix{T}, B::BSplineBasis, x::AbstractVector,
[deriv::Derivative = Derivative(0)];
clip_threshold = eps(T))
Fill preallocated collocation matrix.
See [`collocation_matrix`](@ref) for details.
"""
function collocation_matrix!(
C::AbstractMatrix{T}, B::BSplineBasis, x::AbstractVector;
Ndiff::Val = Val(0), clip_threshold = eps(T)) where {T}
C::AbstractMatrix{T}, B::BSplineBasis, x::AbstractVector,
deriv::Derivative = Derivative(0),
clip_threshold = eps(T)) where {T}
Nx, Nb = size(C)
if Nx != length(x) || Nb != length(B)
......@@ -232,7 +241,7 @@ function collocation_matrix!(
b_lo, b_hi = bandwidths(C)
for j = 1:Nb, i = 1:Nx
b = evaluate_bspline(B, j, x[i], T, Ndiff=Ndiff)
b = evaluate_bspline(B, j, x[i], deriv, T)
# Skip very small values (and zeros).
# This is important for SparseMatrixCSC, which also stores explicit
......
......@@ -7,7 +7,7 @@ The basis is defined by a set of knots and by the B-spline order.
---
BSplineBasis(k::Union{Val,Integer}, knots::Vector; augment=true)
BSplineBasis(k::Union{BSplineOrder,Integer}, knots::Vector; augment=true)
Create B-spline basis of order `k` with the given knots.
......@@ -17,7 +17,7 @@ have multiplicity `k`. See also [`augment_knots`](@ref).
struct BSplineBasis{k, T} # k: B-spline order
N :: Int # number of B-splines ("resolution")
t :: Vector{T} # knots (length = N + k)
function BSplineBasis(::Val{k}, knots::AbstractVector{T};
function BSplineBasis(::BSplineOrder{k}, knots::AbstractVector{T};
augment=true) where {k,T}
k :: Integer
if k <= 0
......@@ -30,7 +30,7 @@ struct BSplineBasis{k, T} # k: B-spline order
end
@inline BSplineBasis(k::Integer, args...; kwargs...) =
BSplineBasis(Val(k), args...; kwargs...)
BSplineBasis(BSplineOrder(k), args...; kwargs...)
# Make BSplineBasis behave as scalar when broadcasting.
Broadcast.broadcastable(B::BSplineBasis) = Ref(B)
......@@ -52,13 +52,15 @@ Returns the knots of the spline grid.
knots(g::BSplineBasis) = g.t
"""
order(::Type{BSplineBasis})
order(::Type{Spline})
order(::Type{BSplineBasis}) -> Int
order(::Type{Spline}) -> Int
order(::BSplineOrder) -> Int
Returns order of B-splines.
Returns order of B-splines as an integer.
"""
order(::Type{<:BSplineBasis{k}}) where {k} = k
order(b::BSplineBasis) = order(typeof(b))
order(::BSplineOrder{k}) where {k} = k
"""
BSpline{B <: BSplineBasis}
......@@ -114,85 +116,83 @@ using `isempty`, which returns `true` for such a range.
common_support(bs::Vararg{BSpline}) = (support.(bs)...)
"""
(b::BSpline)(x, [Ndiff = Val(0)])
(b::BSpline)(x, [deriv = Derivative(0)])
Evaluate B-spline at coordinate `x`.
To evaluate a derivative, pass the `Ndiff` parameter with the wanted
To evaluate a derivative, pass the `deriv` parameter with the wanted
differentiation order.
"""
(b::BSpline)(x, Ndiff=Val(0)) =
evaluate_bspline(basis(b), b.i, x, eltype(b), Ndiff=Ndiff)
(b::BSpline)(x, deriv=Derivative(0)) =
evaluate_bspline(basis(b), b.i, x, deriv, eltype(b))
"""
evaluate_bspline(
B::BSplineBasis, i::Integer, x, [T=Float64];
Ndiff::{Integer, Val} = Val(0),
B::BSplineBasis, i::Integer, x, [deriv=Derivative(0)], [T=Float64]
)
Evaluate i-th B-spline in the given basis at `x` (can be a coordinate or a
vector of coordinates).
The `N`-th derivative of bᵢ(x) may be evaluated by passing `Ndiff = Val(N)`.
The `N`-th derivative of bᵢ(x) may be evaluated by passing `Derivative(N)`.
See also [`evaluate_bspline!`](@ref).
"""
function evaluate_bspline(B::BSplineBasis, i::Integer, x::Real,
::Type{T} = Float64;
Ndiff::Union{Val,Integer} = Val(0)) where {T,D}
deriv::Derivative = Derivative(0),
::Type{T} = Float64) where {T}
N = length(B)
if !(1 <= i <= N)
throw(DomainError(i, "B-spline index must be in 1:$N"))
end
k = order(B)
t = knots(B)
Ndiff_val = _make_val(Ndiff)
evaluate_bspline_diff(Ndiff_val, Val(k), t, i, x, T)
evaluate_bspline_diff(deriv, BSplineOrder(k), t, i, x, T)
end
@inline _make_val(x) = Val(x)
@inline _make_val(x::Val) = x
# No derivative
evaluate_bspline_diff(::Val{0}, ::Val{k}, t, i, x, ::Type{T}) where {k,T} =
_evaluate_bspline(Val(k), t, i, x, T)
evaluate_bspline_diff(::Derivative{0}, ::BSplineOrder{k}, t, i, x,
::Type{T}) where {k,T} =
_evaluate_bspline(BSplineOrder(k), t, i, x, T)
# N-th derivative
function evaluate_bspline_diff(::Val{N}, ::Val{k}, t, i, x,
function evaluate_bspline_diff(::Derivative{N}, ::BSplineOrder{k}, t, i, x,
::Type{T}) where {N,k,T}
@assert N > 0
y = zero(T)
dt = t[i + k - 1] - t[i]
if !iszero(dt)
# Recursively evaluate derivative `N - 1` of B-spline of order `k - 1`.
y += evaluate_bspline_diff(Val(N - 1), Val(k - 1), t, i, x, T) / dt
y += evaluate_bspline_diff(
Derivative(N - 1), BSplineOrder(k - 1), t, i, x, T) / dt
end
dt = t[i + k] - t[i + 1]
if !iszero(dt)
y -= evaluate_bspline_diff(Val(N - 1), Val(k - 1), t, i + 1, x, T) / dt
y -= evaluate_bspline_diff(
Derivative(N - 1), BSplineOrder(k - 1), t, i + 1, x, T) / dt
end
y * (k - 1)
end
evaluate_bspline(B::BSplineBasis, i, x::AbstractVector, args...; kwargs...) =
evaluate_bspline.(B, i, x, args...; kwargs...)
evaluate_bspline(B::BSplineBasis, i, x::AbstractVector, args...) =
evaluate_bspline.(B, i, x, args...)
"""
evaluate_bspline!(b::AbstractVector, B::BSplineBasis, i::Integer,
x::AbstractVector; kwargs...)
x::AbstractVector, args...)
Evaluate i-th B-spline at positions `x` and write result to `b`.
See also [`evaluate_bspline`](@ref).
"""
function evaluate_bspline!(b::AbstractVector{T}, B::BSplineBasis, i,
x::AbstractVector; kwargs...) where {T}
broadcast!(x -> evaluate_bspline(B, i, x, T; kwargs...), b, x)
x::AbstractVector, args...) where {T}
broadcast!(x -> evaluate_bspline(B, i, x, args..., T), b, x)
end
# Specialisation for first order B-splines.
function _evaluate_bspline(::Val{1}, t::AbstractVector, i::Integer, x::Real,
::Type{T}) where {T}
function _evaluate_bspline(::BSplineOrder{1}, t::AbstractVector, i::Integer,
x::Real, ::Type{T}) where {T}
# Local support of the B-spline.
@inbounds ta = t[i]
@inbounds tb = t[i + 1]
......@@ -201,8 +201,8 @@ function _evaluate_bspline(::Val{1}, t::AbstractVector, i::Integer, x::Real,
end
# General case of order k >= 2.
function _evaluate_bspline(::Val{k}, t::AbstractVector, i::Integer, x::Real,
::Type{T}) where {T,k}
function _evaluate_bspline(::BSplineOrder{k}, t::AbstractVector, i::Integer,
x::Real, ::Type{T}) where {T,k}
k::Int
@assert k >= 2
......@@ -221,12 +221,12 @@ function _evaluate_bspline(::Val{k}, t::AbstractVector, i::Integer, x::Real,
y = zero(T)
if tb1 != ta
y += _evaluate_bspline(Val(k - 1), t, i, x, T) *
y += _evaluate_bspline(BSplineOrder(k - 1), t, i, x, T) *
(x - ta) / (tb1 - ta)
end
if ta1 != tb
y += _evaluate_bspline(Val(k - 1), t, i + 1, x, T) *
y += _evaluate_bspline(BSplineOrder(k - 1), t, i + 1, x, T) *
(tb - x) / (tb - ta1)
end
......
"""
galerkin_matrix(
B::BSplineBasis, [MatrixType = BandedMatrix{Float64}];
Ndiff::Val = (Val(0), Val(0)),
B::BSplineBasis,
[deriv = (Derivative(0), Derivative(0))],
[MatrixType = BandedMatrix{Float64}]
)
Compute Galerkin mass or stiffness matrix.
......@@ -16,10 +17,10 @@ Here, ⟨⋅,⋅⟩ is the [L² inner
product](https://en.wikipedia.org/wiki/Square-integrable_function#Properties)
between functions.
To obtain a matrix associated to the B-spline derivatives, set the `Ndiff`
argument to the order of the derivative.
For instance, if `Ndiff = (Val(0), Val(2))`, this returns the matrix
`⟨ bᵢ, bⱼ'' ⟩`.
To obtain a matrix associated to the B-spline derivatives, set the `deriv`
argument to the order of the derivatives.
For instance, if `deriv = (Derivative(0), Derivative(2))`, this returns the
matrix `⟨ bᵢ, bⱼ'' ⟩`.
Note that the Galerkin matrix is banded,
with `k + 1` and `k + 2` for `k` even and odd, respectively.
......@@ -35,12 +36,11 @@ Other types of container are also supported, including regular sparse matrices
(`SparseMatrixCSC`) and dense arrays (`Matrix`).
"""
function galerkin_matrix(
B::BSplineBasis,
::Type{M} = BandedMatrix{Float64};
Ndiff = (Val(0), Val(0)),
B::AnyBSplineBasis,
deriv = Derivative.((0, 0)),
::Type{M} = BandedMatrix{Float64},
) where {M <: AbstractMatrix}
N = length(B)
deriv = _galerkin_make_Ndiff(Ndiff)
symmetry = deriv[1] === deriv[2]
A = allocate_galerkin_matrix(M, N, order(B), symmetry)
......@@ -48,9 +48,12 @@ function galerkin_matrix(
# Make the matrix symmetric if possible.
S = symmetry ? Symmetric(A) : A
galerkin_matrix!(S, B, Ndiff=deriv)
galerkin_matrix!(S, B, deriv)
end
galerkin_matrix(B, ::Type{M}) where {M} =
galerkin_matrix(B, Derivative.((0, 0)), M)
allocate_galerkin_matrix(::Type{M}, N, etc...) where {M <: AbstractMatrix} =
M(undef, N, N)
......@@ -59,30 +62,28 @@ allocate_galerkin_matrix(::Type{SparseMatrixCSC{T}}, N, etc...) where {T} =
function allocate_galerkin_matrix(::Type{M}, N, k,
symmetry) where {M <: BandedMatrix}
# The upper/lower bandwidths are:
# - for even k: Nb = k / 2 (total = k + 1 bands)
# - for odd k: Nb = (k + 1) / 2 (total = k + 2 bands)
# The upper/lower bandwidths are Nb = k (total = 2k + 1 bands).
# Note that if the matrix is also symmetric, then we only need the upper
# band.
Nb = (k + 1) >> 1
Nb = k
bands = symmetry ? (0, Nb) : (Nb, Nb)
M(undef, (N, N), bands)
end
"""
galerkin_matrix!(A::AbstractMatrix, B::BSplineBasis;
Ndiff::Val = (Val(0), Val(0)))
galerkin_matrix!(A::AbstractMatrix, B::BSplineBasis,
deriv = (Derivative(0), Derivative(0)))
Fill preallocated Galerkin matrix.
The matrix may be a `Symmetric` view, in which case only one half of the matrix
will be filled. Note that, for the matrix to be symmetric, both derivative orders
in `Ndiff` must be the same.
in `deriv` must be the same.
See also [`galerkin_matrix`](@ref).
"""
function galerkin_matrix!(S::AbstractMatrix, B::BSplineBasis;
Ndiff = (Val(0), Val(0)))
function galerkin_matrix!(S::AbstractMatrix, B::BSplineBasis,
deriv = Derivative.((0, 0)))
N = size(S, 1)
if N != length(B)
......@@ -93,16 +94,14 @@ function galerkin_matrix!(S::AbstractMatrix, B::BSplineBasis;
k = order(B)
t = knots(B)
h = (k + 1) ÷ 2 # k/2 if k is even
h = k - 1
T = eltype(S)
# Quadrature information (weights, nodes).
quad = _quadrature_prod(k)
deriv = _galerkin_make_Ndiff(Ndiff)
if S isa Symmetric
deriv[1] === deriv[2] || error("matrix will not be symmetric with Ndiff = $Ndiff")
deriv[1] === deriv[2] || error("matrix will not be symmetric with deriv = $deriv")
fill_upper = S.uplo === 'U'
fill_lower = S.uplo === 'L'
A = parent(S)
......@@ -134,8 +133,75 @@ function galerkin_matrix!(S::AbstractMatrix, B::BSplineBasis;
S
end
_galerkin_make_Ndiff(v::Val) = (v, v)
_galerkin_make_Ndiff(v::Tuple{Vararg{<:Val,2}}) = v
# TODO merge this with other variant!
function galerkin_matrix!(S::AbstractMatrix, B::BSplines.BSplineBasis,
deriv = Derivative.((0, 0)))
N = size(S, 1)
if N != length(B)
throw(ArgumentError("wrong dimensions of Galerkin matrix"))
end
fill!(S, 0)
k = order(B)
t = knots(B)
T = eltype(S)
# Quadrature information (weights, nodes).
quad = _quadrature_prod(k)
same_deriv = deriv[1] === deriv[2]
if S isa Symmetric
same_deriv || error("matrix will not be symmetric with deriv = $deriv")
fill_upper = S.uplo === 'U'
fill_lower = S.uplo === 'L'
A = parent(S)
else
fill_upper = true
fill_lower = true
A = S
end
# B-splines evaluated at a given point.
bi = zeros(T, k)
bj = copy(bi)
# Integrate over each segment between knots.
for l = k:N
a, b = t[l], t[l + 1] # segment (t[l], t[l + 1])
α = (b - a) / 2
β = (a + b) / 2
x, w = quad
for n in eachindex(w)
y = α * x[n] + β
# Evaluate B-splines at quadrature point.
# This evaluates B-splines (i, j) ∈ [(off + 1):(off + k)].
off = BSplines.bsplines!(bi, B, y, deriv[1], leftknot=l) :: Int
if same_deriv
copy!(bj, bi)
else
δ = BSplines.bsplines!(bj, B, y, deriv[2], leftknot=l) :: Int
@assert δ === off
end
@assert off === l - k
for j = 1:k
istart = fill_upper ? 1 : j
iend = fill_lower ? k : j
for i = istart:iend
A[off + i, off + j] += α * w[n] * bi[i] * bj[j]
end
end
end
end
S
end
# Generate quadrature information for B-spline product.
# Returns weights and nodes for integration in [-1, 1].
......
"""
augment_knots(knots::AbstractVector, k::Integer)
augment_knots(knots::AbstractVector, k::Union{Integer,BSplineOrder})
Modifies the input knots to make sure that the first and last knot have
the multiplicity `k` for splines of order `k`.
......@@ -33,6 +33,8 @@ function augment_knots(knots::AbstractVector{T},
t
end
augment_knots(knots, ::BSplineOrder{k}) where {k} = augment_knots(knots, k)
"""
multiplicity(knots, i)
......
......@@ -74,11 +74,12 @@ function (S::Spline)(x)
end
"""
diff(S::Spline, [N::Union{Val, Integer} = Val(1)]) -> Spline
diff(S::Spline, [deriv::Derivative = Derivative(1)]) -> Spline
Return `N`-th derivative of spline `S` as a new spline.
"""
function Base.diff(S::Spline, ::Val{Ndiff} = Val(1)) where {Ndiff}
function Base.diff(S::Spline,
::Derivative{Ndiff} = Derivative(1)) where {Ndiff}
Ndiff :: Integer
@assert Ndiff >= 1
......@@ -112,14 +113,13 @@ function Base.diff(S::Spline, ::Val{Ndiff} = Val(1)) where {Ndiff}
N = length(u)
Nt = length(t)
t_new = view(t, (1 + Ndiff):(Nt - Ndiff))
B = BSplineBasis(Val(k - Ndiff), t_new, augment=false)
B = BSplineBasis(BSplineOrder(k - Ndiff), t_new, augment=false)
Spline(B, view(du, (1 + Ndiff):N))
end
# Zeroth derivative: return S itself.
@inline Base.diff(S::Spline, ::Val{0}) = S
@inline Base.diff(S::Spline, Ndiff::Integer) = diff(S, Val(Ndiff))
@inline Base.diff(S::Spline, ::Derivative{0}) = S
"""
integral(S::Spline)
......@@ -154,7 +154,7 @@ function integral(S::Spline)
end
end
B = BSplineBasis(Val(k + 1), t_int, augment=false)
B = BSplineBasis(BSplineOrder(k + 1), t_int, augment=false)
Spline(B, β)
end
......
......@@ -78,12 +78,12 @@ function test_splines(B::BSplineBasis, knots_in)
nothing
end
function test_splines(::Val{k}) where {k}
function test_splines(::BSplineOrder{k}) where {k}
knots_in = let N = 10 + k
[-cos(n * π / N) for n = 0:N]
end
@inferred BSplineBasis(Val(k), knots_in)
@inferred BSplineBasis(BSplineOrder(k), knots_in)
@inferred (() -> BSplineBasis(k, knots_in))()
g = BSplineBasis(k, knots_in)
......@@ -94,8 +94,8 @@ function test_splines(::Val{k}) where {k}
end
function main()
test_splines(Val(3))
test_splines(Val(4))
test_splines(BSplineOrder(3))
test_splines(BSplineOrder(4))
@testset "Galerkin" begin
test_galerkin()
end
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment