Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow creating Galerkin matrices with custom functions in integrand #82

Merged
merged 3 commits into from
Dec 1, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "BSplineKit"
uuid = "093aae92-e908-43d7-9660-e50ee39d5a0a"
authors = ["Juan Ignacio Polanco <juan-ignacio.polanco@cnrs.fr>"]
version = "0.16.5"
version = "0.17.0"

[deps]
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
Expand Down
2 changes: 1 addition & 1 deletion src/Galerkin/Galerkin.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ using ..Recombinations: num_constraints
using ..Collocation: collocation_matrix # for Documenter

using BandedMatrices
using FastGaussQuadrature: gausslegendre
using FastGaussQuadrature: FastGaussQuadrature
using LinearAlgebra: Hermitian, ⋅
using SparseArrays
using StaticArrays
Expand Down
86 changes: 64 additions & 22 deletions src/Galerkin/linear.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
"""
galerkin_matrix(
[f::Function],
B::AbstractBSplineBasis,
[deriv = (Derivative(0), Derivative(0))],
[MatrixType = BandedMatrix{Float64}],
[derivatives = (Derivative(0), Derivative(0))],
[MatrixType = BandedMatrix{Float64}];
[quadrature = default_quadrature((B, B))],
)

Compute Galerkin mass or stiffness matrix, as well as more general variants
Expand Down Expand Up @@ -55,8 +57,8 @@ See [`collocation_matrix`](@ref) for a discussion on matrix types.
## Derivatives of basis functions

Galerkin matrices associated to the derivatives of basis functions may be
constructed using the optional `deriv` parameter.
For instance, if `deriv = (Derivative(0), Derivative(2))`, the matrix
constructed using the optional `derivatives` parameter.
For instance, if `derivatives = (Derivative(0), Derivative(2))`, the matrix
``⟨ ϕ_i, ϕ_j'' ⟩`` is constructed, where primes denote derivatives.
Note that, if the derivative orders are different, the resulting matrix is not
symmetric, and a `Hermitian` view is not returned in those cases.
Expand All @@ -77,35 +79,66 @@ will not even be square if the bases have different lengths.

In practice, this feature may be used to combine a B-spline basis `B`, with a
recombined basis `R` generated from `B` (see [Basis recombination](@ref basis-recombination-api)).

## Integrating more general functions

Say we wanted to integrate the more general term

```math
L_{ij} = ⟨ f(x) ϕ_i(x), ϕ_j(x) ⟩ = ∫_{a}^{b} f(x) ϕ_i(x) ϕ_j(x) \\, \\mathrm{d}x
```

To obtain an approximation of this matrix, one would pass the ``f(x)`` function as a first
positional argument to `galerkin_matrix` (or [`galerkin_matrix!`](@ref)).
This can also be combined with the `derivatives` argument if one wants to consider
derivatives of ``ϕ_i`` or ``ϕ_j``.

Note that, in this case, the computation of the integrals is not guaranteed to be exact,
since Gauss--Legendre quadratures are only "exact" when the integrand is a polynomial (the
product of two B-splines is a piecewise polynomial).
To improve accuracy, one may want to increase the number `n` of quadrature nodes.
For this, pass `quadrature = Galerkin.gausslegendre(Val(n))` as a keyword argument.

"""
function galerkin_matrix end

function galerkin_matrix(
f::F,
Bs::NTuple{2,AbstractBSplineBasis},
deriv::DerivativeCombination{2} = Derivative.((0, 0)),
::Type{M} = _default_matrix_type(first(Bs)),
) where {M <: AbstractMatrix}
::Type{M} = _default_matrix_type(first(Bs));
kwargs...,
) where {F <: Function, M <: AbstractMatrix}
B1, B2 = Bs
symmetry = M <: Hermitian
if symmetry && (deriv[1] !== deriv[2] || B1 !== B2)
throw(ArgumentError(
"input matrix type incorrectly assumes symmetry"))
end
A = allocate_galerkin_matrix(M, Bs)
galerkin_matrix!(A, Bs, deriv)
galerkin_matrix!(f, A, Bs, deriv; kwargs...)
end

function galerkin_matrix(
f::F,
B::AbstractBSplineBasis,
deriv::DerivativeCombination{2} = Derivative.((0, 0)),
::Type{Min} = _default_matrix_type(B),
) where {Min <: AbstractMatrix}
::Type{Min} = _default_matrix_type(B);
kwargs...,
) where {F <: Function, Min <: AbstractMatrix}
symmetry = deriv[1] === deriv[2]
T = eltype(Min)
M = symmetry ? Hermitian{T,Min} : Min
galerkin_matrix((B, B), deriv, M)
galerkin_matrix(f, (B, B), deriv, M; kwargs...)
end

galerkin_matrix(B, ::Type{M}) where {M <: AbstractMatrix} =
galerkin_matrix(B, Derivative.((0, 0)), M)
galerkin_matrix(f::F, B, ::Type{M}; kwargs...) where {F <: Function, M <: AbstractMatrix} =
galerkin_matrix(f, B, Derivative.((0, 0)), M; kwargs...)

galerkin_matrix(B::AbstractBSplineBasis, args...; kws...) =
galerkin_matrix(Returns(true), B, args...; kws...)
galerkin_matrix(B::NTuple{2, AbstractBSplineBasis}, args...; kws...) =
galerkin_matrix(Returns(true), B, args...; kws...)

_default_matrix_type(::Type{<:AbstractBSplineBasis}) = BandedMatrix{Float64}
_default_matrix_type(::Type{<:PeriodicBSplineBasis}) = SparseMatrixCSC{Float64}
Expand Down Expand Up @@ -158,8 +191,10 @@ function allocate_galerkin_matrix(::Type{M}, Bs, symmetry) where {M <: BandedMat
end

"""
galerkin_matrix!(A::AbstractMatrix, B::AbstractBSplineBasis,
deriv = (Derivative(0), Derivative(0)))
galerkin_matrix!(
[f::Function], A::AbstractMatrix, B::AbstractBSplineBasis, [deriv = (Derivative(0), Derivative(0))];
[quadrature],
)

Fill preallocated Galerkin matrix.

Expand All @@ -170,17 +205,22 @@ in `deriv` must be the same.
More generally, it is possible to combine different functional bases by passing
a tuple of `AbstractBSplineBasis` as `B`.

See [`galerkin_matrix`](@ref) for details.
See [`galerkin_matrix`](@ref) for details on arguments.
"""
function galerkin_matrix! end

galerkin_matrix!(M, B::AbstractBSplineBasis, args...) =
galerkin_matrix!(M, (B, B), args...)
galerkin_matrix!(M::AbstractMatrix, args...; kwargs...) =
galerkin_matrix!(Returns(true), M, args...; kwargs...)

galerkin_matrix!(f::F, M::AbstractMatrix, B::AbstractBSplineBasis, args...; kwargs...) where {F <: Function} =
galerkin_matrix!(f, M, (B, B), args...; kwargs...)

function galerkin_matrix!(
f::F,
S::AbstractMatrix, Bs::Tuple{Vararg{AbstractBSplineBasis,2}},
deriv = Derivative.((0, 0)),
)
deriv = Derivative.((0, 0));
quadrature = default_quadrature(Bs),
) where {F <: Function}
_check_bases(Bs)
B1, B2 = Bs
same_ij = B1 == B2 && deriv[1] == deriv[2]
Expand All @@ -198,8 +238,10 @@ function galerkin_matrix!(
@assert k == order(B2) && ts === knots(B2)

# Quadrature information (nodes, weights).
quadx, quadw = _quadrature_prod(Val(2k - 2))
@assert length(quadx) == k # we need k quadrature points per knot segment
quadx, quadw = quadrature
if length(quadx) < k
@warn lazy"quadrature: it is recommended use at least k = $k quadrature nodes (using n = $(length(quadx)))"
end

A = if S isa Hermitian
deriv[1] === deriv[2] ||
Expand Down Expand Up @@ -262,7 +304,7 @@ function galerkin_matrix!(
@assert iszero(bj)
continue
end
@inbounds A[i, j] += y * bi * bj
@inbounds A[i, j] += y * bi * bj * f(x) # if `f` was not passed, the default is f(x) = true (= 1)
end
end
end
Expand Down
12 changes: 9 additions & 3 deletions src/Galerkin/quadratures.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,19 @@
# Here, p is the polynomial order (p = 2k - 2 for the product of two B-splines).
function _quadrature_prod(::Val{p}) where {p}
n = cld(p + 1, 2)
_gausslegendre(Val(n))
gausslegendre(Val(n))
end

function default_quadrature(Bs::Tuple{Vararg{AbstractBSplineBasis}})
# Polynomial order of integrand (the value is usually inferred).
polynomial_order = sum(B -> order(B) - 1, Bs)
_quadrature_prod(Val(polynomial_order))
end

# Compile-time computation of Gauss--Legendre nodes and weights.
# If the @generated branch is taken, then the computation is done at compile
# time, making sure that no runtime allocations are performed.
function _gausslegendre(::Val{n}) where {n}
function gausslegendre(::Val{n}) where {n}
if @generated
vecs = _gausslegendre_impl(Val(n))
:( $vecs )
Expand All @@ -35,7 +41,7 @@ function _gausslegendre(::Val{n}) where {n}
end

function _gausslegendre_impl(::Val{n}) where {n}
data = gausslegendre(n)
data = FastGaussQuadrature.gausslegendre(n)
map(SVector{n}, data)
end

Expand Down
40 changes: 37 additions & 3 deletions test/galerkin.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@ function test_galerkin_analytical()
t = knots(B)
G = galerkin_matrix(B)

@inferred Galerkin.default_quadrature((B, B))

let i = 8
@test G[i, i] ≈ (t[i + 2] - t[i]) / 3
@test G[i - 1, i] ≈ (t[i + 1] - t[i]) / 6
Expand All @@ -31,6 +33,40 @@ function test_galerkin_analytical()
nothing
end

function test_galerkin_with_function()
N = 32
xs = [-cospi(n / N) for n = 0:N]
f(x) = x^2
B = BSplineBasis(BSplineOrder(4), copy(xs))
G_base = @inferred galerkin_matrix(B)
G = @inferred galerkin_matrix(f, B)

# In all cases, multiplying the integrand by x^2 should result in a smaller integral
# value, since x ∈ [-1, 1].
@test all(zip(G_base.data, G.data)) do (u, v)
iszero(u) || u > v
end

# By default quadratures use `k` quadrature nodes, which is enough to "exactly"
# integrate a product of two B-splines in a knot segment (polynomial of degree 2k - 2).
k = order(B)
quad_base = Galerkin.default_quadrature((B, B))
n_base = length(quad_base[1])
@test n_base == k

# However this quadrature is no longer exact when we multiply by an extra term.
# If we multiply by x^2, then we have a polynomial of degree 2k, and we need (k + 1)
# quadrature nodes.
G_improved = @inferred galerkin_matrix(f, B; quadrature = Galerkin.gausslegendre(Val(k + 1)))
@test maximum(abs, G - G_improved) > 1e-8 # values slightly changed!

# Further increasing the number of quadrature nodes is unnecessary...
G2 = @inferred galerkin_matrix(f, B; quadrature = Galerkin.gausslegendre(Val(k + 2)))
@test maximum(abs, G2 - G_improved) < 1e-16 # basically nothing changed!

nothing
end

function test_galerkin_projection(B)
# Idea: project a spline with all coefficients set to c_i = 1.
# The resulting projection should be equal to the sum of the rows (or
Expand Down Expand Up @@ -260,12 +296,10 @@ end

@testset "Galerkin" begin
test_galerkin_analytical()

test_galerkin_with_function()
gauss_lobatto_points(N) = [-cos(π * n / N) for n = 0:N]

knots_base = gauss_lobatto_points(41)
B = BSplineBasis(5, copy(knots_base))
test_galerkin_projection(B)

test_galerkin_recombined()
end