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

Start work on Galerkin matrices

parent 294b3be6
......@@ -5,6 +5,7 @@ version = "0.1.0"
[deps]
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
......@@ -21,12 +21,15 @@ export augment_knots
export evaluate_bspline, evaluate_bspline!
export integral
using BandedMatrices
using Reexport
using LinearAlgebra: Symmetric
using StaticArrays: MVector
include("knots.jl")
include("basis.jl")
include("spline.jl")
include("galerkin.jl")
include("Collocation.jl")
......
......@@ -112,8 +112,12 @@ end
"""
collocation_matrix(
B::BSplineBasis, x::AbstractVector, [MatrixType=SparseMatrixCSC{Float64}];
Ndiff::Val = Val(0), clip_threshold = eps(eltype(MatrixType)))
B::BSplineBasis,
x::AbstractVector,
[MatrixType=SparseMatrixCSC{Float64}];
Ndiff::Val = Val(0),
clip_threshold = eps(eltype(MatrixType)),
)
Return banded collocation matrix mapping B-spline coefficients to spline values
at the collocation points `x`.
......@@ -219,9 +223,11 @@ function collocation_matrix!(
C::AbstractMatrix{T}, B::BSplineBasis, x::AbstractVector;
Ndiff::Val = Val(0), clip_threshold = eps(T)) where {T}
Nx, Nb = size(C)
if Nx != length(x) || Nb != length(B)
throw(ArgumentError("wrong dimensions of collocation matrix"))
end
fill!(C, 0)
b_lo, b_hi = bandwidths(C)
......
function galerkin_matrix(
B::BSplineBasis,
::Type{M} = BandedMatrix{Float64}
) where {M <: AbstractMatrix}
N = length(B)
A = allocate_galerkin_matrix(M, N, order(B))
galerkin_matrix!(A, B)
end
allocate_galerkin_matrix(::Type{M}, N, k) where {M <: AbstractMatrix} =
Symmetric(M(undef, N, N))
function allocate_galerkin_matrix(::Type{M}, N, k) 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)
# Note that the matrix is also symmetric, so we only need the upper band.
Nb = (k + 1) >> 1
A = M(undef, (N, N), (0, Nb))
Symmetric(A)
end
function galerkin_matrix!(S::Symmetric, B::BSplineBasis)
N = size(S, 1)
if N != length(B)
throw(ArgumentError("wrong dimensions of Galerkin matrix"))
end
fill!(S, 0)
# The matrix is symmetric, so we fill only the upper part.
# For now we assume that S uses the upper part of its parent.
@assert S.uplo === 'U'
A = parent(S)
k = order(B)
h = (k + 1) >> 2 # k/2 if k is even
# Upper part: j >= i
for j = 1:N
# We're only visiting the elements that have non-zero values.
# In other words, we know that S[i, j] = 0 outside the chosen interval.
istart = clamp(j - h, 1, N)
bj = BSpline(B, j)
for i = istart:j
bi = BSpline(B, i)
end
end
S
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