Allow creating Galerkin matrices with custom functions in integrand (#82

)

* Allow passing custom function to galerkin_matrix

* Add tests

* Bump version

Project.toml

@@ -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"

src/Galerkin/Galerkin.jl

@@ -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

src/Galerkin/linear.jl

@@ -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
@@ -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.
@@ -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}
@@ -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.
 
@@ -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]
@@ -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] ||
@@ -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

src/Galerkin/quadratures.jl

@@ -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 )
@@ -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
 

test/galerkin.jl

@@ -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
@@ -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
@@ -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