Commit 019f0ea8 authored by Juan Ignacio Polanco's avatar Juan Ignacio Polanco

Define and use BSplineOrder type

parent 67b66f70
......@@ -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!
......@@ -32,7 +32,15 @@ using StaticArrays: MVector
# We're transitioning to using the registered BSplines package...
import BSplines
using BSplines: Derivative
export Derivative
"""
BSplineOrder(k::Integer)
Specifies the B-spline order `k`.
"""
struct BSplineOrder{k} end
BSplineOrder(k::Integer) = BSplineOrder{k}()
include("knots.jl")
include("basis.jl")
......
......@@ -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,8 +52,8 @@ 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
Returns order of B-splines.
"""
......@@ -145,15 +145,16 @@ function evaluate_bspline(B::BSplineBasis, i::Integer, x::Real,
end
k = order(B)
t = knots(B)
evaluate_bspline_diff(deriv, Val(k), t, i, x, T)
evaluate_bspline_diff(deriv, BSplineOrder(k), t, i, x, T)
end
# No derivative
evaluate_bspline_diff(::Derivative{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(::Derivative{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)
......@@ -161,12 +162,12 @@ function evaluate_bspline_diff(::Derivative{N}, ::Val{k}, t, i, x,
if !iszero(dt)
# Recursively evaluate derivative `N - 1` of B-spline of order `k - 1`.
y += evaluate_bspline_diff(
Derivative(N - 1), Val(k - 1), t, i, x, T) / dt
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(
Derivative(N - 1), Val(k - 1), t, i + 1, x, T) / dt
Derivative(N - 1), BSplineOrder(k - 1), t, i + 1, x, T) / dt
end
y * (k - 1)
end
......@@ -188,8 +189,8 @@ function evaluate_bspline!(b::AbstractVector{T}, B::BSplineBasis, i,
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]
......@@ -198,8 +199,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
......@@ -218,12 +219,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
......
"""
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)
......
......@@ -113,7 +113,7 @@ function Base.diff(S::Spline,
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
......@@ -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