Commit 87476a16 authored by Juan Ignacio Polanco's avatar Juan Ignacio Polanco

Define RecombinedBSplineBasis

parent f23d5b70
......@@ -15,7 +15,7 @@ module BasisSplines
export Collocation
export AbstractBSplineBasis, BSplineBasis
export AbstractBSplineBasis, BSplineBasis, RecombinedBSplineBasis
export Spline, BSpline, BSplineOrder, Derivative
export knots, order, coefficients
export augment_knots
......@@ -46,6 +46,7 @@ BSplineOrder(k::Integer) = BSplineOrder{k}()
include("knots.jl")
include("basis.jl")
include("recombined.jl")
const AnyBSplineBasis = Union{BSplineBasis, BSplines.BSplineBasis}
......
......@@ -79,17 +79,20 @@ function collocation_points!(
collocation_points!(method, x, B)
end
# TODO make this work with RecombinedBSplineBasis (remove boundaries!)
function collocation_points!(::AvgKnots, x, B)
N = length(B)
function collocation_points!(::AvgKnots, x, B::AbstractBSplineBasis)
N = length(B) # number of functions in basis
@assert length(x) == N
k = order(B)
t = knots(B)
j = 0
T = eltype(x)
j = 0
if B isa RecombinedBSplineBasis
j += 1 # skip boundaries
end
v::T = inv(k - 1)
a::T, b::T = t[k], t[N + 1] # domain boundaries
a::T, b::T = t[k], t[end - k + 1] # domain boundaries
for i in eachindex(x)
j += 1 # generally i = j, unless x has weird indexation
......@@ -119,7 +122,6 @@ end
x::AbstractVector,
[deriv::Derivative = Derivative(0)],
[MatrixType = SparseMatrixCSC{Float64}];
bc = nothing,
clip_threshold = eps(eltype(MatrixType)),
)
......@@ -178,52 +180,17 @@ bandwidth of the matrix may be larger than the expected bandwidth.
- `Matrix{T}`: a regular dense matrix. Generally performs worse than the
alternatives, especially for large problems.
## Boundary conditions
Basis recombination (see Boyd 2000, ch. 6) is a common technique for imposing
boundary conditions (BCs) in Galerkin methods. In this approach, the basis is
"recombined" so that each basis function individually satisfies the BCs.
The new basis, `{ϕⱼ(x), j ∈ 1:(N-2)}`, has two fewer functions than the original
B-spline basis, `{bⱼ(x), j ∈ 1:N}`. Due to this, the number of collocation
points needed to obtain a square collocation matrix is `N - 2`. In particular,
for the matrix to be invertible, there must be **no** collocation points at the
boundaries.
Thanks to the local support of B-splines, basis recombination involves only a
little portion of the original B-spline basis. For instance, since there is only
one B-spline that is non-zero at each boundary, removing that function from the
basis is enough to apply homogeneous Dirichlet BCs. Imposing BCs for derivatives
is a bit more complex, but not much.
The optional keyword argument `bc` allows specifying homogeneous BCs. It accepts
a `Derivative` object, which sets the order of the derivative to be imposed with
homogeneous BCs. Some typical choices are:
- `bc = Derivative(0)` sets homogeneous Dirichlet BCs (`u = 0` at the
boundaries) by removing the first and last B-splines, i.e. ϕ₁ = b₂;
- `bc = Derivative(1)` sets homogeneous Neumann BCs (`du/dx = 0` at the
boundaries) by adding the two first (and two last) B-splines,
i.e. ϕ₁ = b₁ + b₂.
For now, the two boundaries are given the same BC (but this could be easily
extended...). And actually, only the two above choices are available for now!
See also [`collocation_matrix!`](@ref).
"""
function collocation_matrix(
B::AbstractBSplineBasis, x::AbstractVector,
deriv::Derivative = Derivative(0),
::Type{MatrixType} = SparseMatrixCSC{Float64};
bc=nothing, kwargs...) where {MatrixType}
kwargs...) where {MatrixType}
Nx = length(x)
Nb = length(B)
if bc !== nothing
Nb -= 2 # remove two basis functions if BCs are applied
end
C = allocate_collocation_matrix(MatrixType, (Nx, Nb), order(B); kwargs...)
collocation_matrix!(C, B, x, deriv; bc=bc, kwargs...)
collocation_matrix!(C, B, x, deriv; kwargs...)
end
collocation_matrix(B, x, ::Type{M}; kwargs...) where {M} =
......@@ -257,8 +224,7 @@ end
"""
collocation_matrix!(
C::AbstractMatrix{T}, B::AbstractBSplineBasis, x::AbstractVector,
[deriv::Derivative = Derivative(0)];
bc = nothing, clip_threshold = eps(T))
[deriv::Derivative = Derivative(0)]; clip_threshold = eps(T))
Fill preallocated collocation matrix.
......@@ -267,12 +233,10 @@ See [`collocation_matrix`](@ref) for details.
function collocation_matrix!(
C::AbstractMatrix{T}, B::AbstractBSplineBasis, x::AbstractVector,
deriv::Derivative = Derivative(0);
bc = nothing,
clip_threshold = eps(T)) where {T}
Nx, Nb = size(C)
with_bc = bc !== nothing
if Nx != length(x) || Nb != length(B) - 2 * with_bc
if Nx != length(x) || Nb != length(B)
throw(ArgumentError("wrong dimensions of collocation matrix"))
end
......@@ -280,7 +244,7 @@ function collocation_matrix!(
b_lo, b_hi = bandwidths(C)
for j = 1:Nb, i = 1:Nx
b = recombine_bspline(bc, B, j, x[i], deriv, T)
b = evaluate_bspline(B, j, x[i], deriv, T)
# Skip very small values (and zeros).
# This is important for SparseMatrixCSC, which also stores explicit
......@@ -300,27 +264,4 @@ function collocation_matrix!(
C
end
# TODO define RecombinedBSplineBasis?
# Perform basis recombination for applying homogeneous BCs.
# No basis recombination.
recombine_bspline(::Nothing, B, j, args...) =
evaluate_bspline(B, j, args...)
# For homogeneous Dirichlet BCs: just shift the B-spline basis (removing b₁).
recombine_bspline(::Derivative{0}, B, j, args...) =
evaluate_bspline(B, j + 1, args...)
# Homogeneous Neumann BCs.
function recombine_bspline(::Derivative{1}, B, j, args...)
Nb = length(B) - 2 # length of recombined basis
if j == 1
evaluate_bspline(B, 1, args...) + evaluate_bspline(B, 2, args...)
elseif j == Nb
evaluate_bspline(B, Nb + 1, args...) + evaluate_bspline(B, Nb + 2, args...)
else
# Same as for Dirichlet.
recombine_bspline(Derivative(0), B, j, args...)
end
end
end
end # module
......@@ -43,7 +43,7 @@ end
BSplineBasis(BSplineOrder(k), args...; kwargs...)
# Make BSplineBasis behave as scalar when broadcasting.
Broadcast.broadcastable(B::BSplineBasis) = Ref(B)
Broadcast.broadcastable(B::AbstractBSplineBasis) = Ref(B)
"""
length(g::BSplineBasis)
......@@ -51,7 +51,8 @@ Broadcast.broadcastable(B::BSplineBasis) = Ref(B)
Returns the number of B-splines composing a spline.
"""
Base.length(g::BSplineBasis) = g.N
Base.size(g::BSplineBasis) = (g.N, )
Base.size(g::AbstractBSplineBasis) = (length(g), )
Base.parent(g::BSplineBasis) = g
"""
knots(g::BSplineBasis)
......@@ -138,7 +139,7 @@ differentiation order.
"""
evaluate_bspline(
B::BSplineBasis, i::Integer, x, [deriv=Derivative(0)], [T=Float64]
B::AbstractBSplineBasis, 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
......@@ -195,7 +196,7 @@ 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,
function evaluate_bspline!(b::AbstractVector{T}, B::AbstractBSplineBasis, i,
x::AbstractVector, args...) where {T}
broadcast!(x -> evaluate_bspline(B, i, x, args..., T), b, x)
end
......
"""
RecombinedBSplineBasis{D, k, T}
Functional basis defined from the recombination of a [`BSplineBasis`](@ref)
in order to satisfy certain homogeneous boundary conditions (BCs).
The basis recombination technique is a common way of applying BCs in Galerkin
methods. It is described for instance in Boyd 2000 (ch. 6), in the context of
a Chebyshev basis. In this approach, the original basis is "recombined" so that
each basis function individually satisfies the BCs.
The new basis, `{ϕⱼ(x), j ∈ 1:(N-2)}`, has two fewer functions than the original
B-spline basis, `{bⱼ(x), j ∈ 1:N}`. Due to this, the number of collocation
points needed to obtain a square collocation matrix is `N - 2`. In particular,
for the matrix to be invertible, there must be **no** collocation points at the
boundaries.
Thanks to the local support of B-splines, basis recombination involves only a
little portion of the original B-spline basis. For instance, since there is only
one B-spline that is non-zero at each boundary, removing that function from the
basis is enough to apply homogeneous Dirichlet BCs. Imposing BCs for derivatives
is a bit more complex, but not much.
# Order of boundary condition
The type parameter `D` represents the order of the BC. The recombined basis
requires the specification of a `Derivative` object determining the order of the
homogeneous BCs to be applied at the two boundaries.
Some typical choices are:
- `Derivative(0)` sets homogeneous Dirichlet BCs (`u = 0` at the
boundaries) by removing the first and last B-splines, i.e. ϕ₁ = b₂;
- `Derivative(1)` sets homogeneous Neumann BCs (`du/dx = 0` at the
boundaries) by adding the two first (and two last) B-splines,
i.e. ϕ₁ = b₁ + b₂.
For now, the two boundaries are given the same BC (but this could be easily
extended...). And actually, only the two above choices are available for now.
---
RecombinedBSplineBasis(::Derivative{D}, B::BSplineBasis)
Construct `RecombinedBSplineBasis` from B-spline basis `B`, satisfying
homogeneous boundary conditions of order `D >= 0`.
"""
struct RecombinedBSplineBasis{
D, k, T, Parent <: BSplineBasis{k,T},
} <: AbstractBSplineBasis{k,T}
B :: Parent # original B-spline basis
function RecombinedBSplineBasis(
::Derivative{D}, B::BSplineBasis{k,T}) where {k,T,D}
Parent = typeof(B)
new{D,k,T,Parent}(B)
end
end
"""
RecombinedBSplineBasis(::Derivative{D}, args...; kwargs...)
Construct [`RecombinedBSplineBasis`](@ref) from B-spline basis, satisfying
homogeneous boundary conditions of order `D >= 0`.
This variant does not require a previously constructed [`BSplineBasis`](@ref).
Arguments are passed to the `BSplineBasis` constructor.
"""
RecombinedBSplineBasis(order::Derivative, args...; kwargs...) =
RecombinedBSplineBasis(order, BSplineBasis(args...; kwargs...))
"""
parent(R::RecombinedBSplineBasis)
Get original B-spline basis.
"""
Base.parent(R::RecombinedBSplineBasis) = R.B
"""
length(R::RecombinedBSplineBasis)
Returns the number of functions in the recombined basis.
"""
Base.length(R::RecombinedBSplineBasis) = length(parent(R)) - 2
knots(R::RecombinedBSplineBasis) = knots(parent(R))
order(R::RecombinedBSplineBasis{D,k}) where {D,k} = k
Base.eltype(::Type{RecombinedBSplineBasis{D,k,T}}) where {D,k,T} = T
# For homogeneous Dirichlet BCs: just shift the B-spline basis (removing b₁).
evaluate_bspline(R::RecombinedBSplineBasis{0}, j, args...) =
evaluate_bspline(parent(R), j + 1, args...)
# Homogeneous Neumann BCs.
function evaluate_bspline(R::RecombinedBSplineBasis, j, args...)
Nb = length(R) # length of recombined basis
B = parent(R)
if j == 1
evaluate_bspline(B, 1, args...) + evaluate_bspline(B, 2, args...)
elseif j == Nb
evaluate_bspline(B, Nb + 1, args...) + evaluate_bspline(B, Nb + 2, args...)
else
# Same as for Dirichlet.
evaluate_bspline(B, j + 1, args...)
end
end
......@@ -21,13 +21,19 @@ function test_collocation(B::BSplineBasis, xcol, ::Type{T} = Float64) where {T}
C = collocation_matrix(B, xcol)
N = length(B)
# Recombined basis for Dirichlet BCs (u = 0)
@inferred RecombinedBSplineBasis(Derivative(0), B)
R0 = RecombinedBSplineBasis(Derivative(0), B)
@test length(R0) == N - 2
# Collocation points for recombined basis, for implicit boundary
# conditions (points at the boundaries are removed!).
x = @view xcol[2:end-1]
x = collocation_points(R0)
@test length(x) == N - 2
@test x == @view xcol[2:end-1]
# Dirichlet BCs (u = 0)
@inferred collocation_matrix(B, x, bc=Derivative(0))
C0 = collocation_matrix(B, x, bc=Derivative(0))
@inferred collocation_matrix(R0, x)
C0 = collocation_matrix(R0, x)
@test size(C0) == (N - 2, N - 2)
# C0 is simply the interior elements of C.
......@@ -36,7 +42,9 @@ function test_collocation(B::BSplineBasis, xcol, ::Type{T} = Float64) where {T}
end
# Neumann BCs (du/dx = 0)
C1 = collocation_matrix(B, x, bc=Derivative(1))
R1 = RecombinedBSplineBasis(Derivative(1), B)
@test collocation_points(R1) == x # same as for Dirichlet
C1 = collocation_matrix(R1, x)
@test size(C1) == (N - 2, N - 2)
let r = 2:(N - 1)
......@@ -87,6 +95,8 @@ function test_splines(B::BSplineBasis, knots_in)
end
xcol = collocation_points(B, method=Collocation.AvgKnots())
@test xcol[1] == knots_in[1]
@test xcol[end] == knots_in[end]
@testset "Collocation (k = $k)" begin
test_collocation(B, xcol)
......
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