Galerkin arrays

galerkin_matrix(
    B::AbstractBSplineBasis,
    [deriv = (Derivative(0), Derivative(0))],
    [MatrixType = BandedMatrix{Float64}],
)

Compute Galerkin mass or stiffness matrix, as well as more general variants of these.

Extended help

The Galerkin mass matrix is defined as

\[M_{ij} = ⟨ ϕ_i, ϕ_j ⟩ \quad \text{for} \quad i ∈ [1, N] \text{ and } j ∈ [1, N],\]

where $ϕ_i(x)$ is the $i$-th basis function and N = length(B) is the number of functions in the basis B. Here, $⟨f, g⟩$ is the $L^2$ inner product between functions $f$ and $g$.

Since products of B-splines are themselves piecewise polynomials, integrals can be computed exactly using Gaussian quadrature rules. To do this, we use Gauss–Legendre quadratures via the FastGaussQuadrature package.

Matrix layout and types

The mass matrix is banded with $2k - 1$ bands. Moreover, the matrix is symmetric and positive definite, and only $k$ bands are needed to fully describe the matrix. Hence, a Hermitian view of an underlying matrix is returned.

By default, the underlying matrix holding the data is a BandedMatrix that defines the upper part of the symmetric matrix. Other types of container are also supported, including regular sparse matrices (SparseMatrixCSC) and dense arrays (Matrix). See collocation_matrix 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 $⟨ ϕ_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.

Combining different bases

More generally, it is possible to compute matrices of the form $⟨ ψ_i^{(n)}, ϕ_j^{(m)} ⟩$, where n and m are derivative orders, and $ψ_i$ and $ϕ_j$ belong to two different (but related) bases B₁ and B₂. For this, instead of the B parameter, one must pass a tuple of bases (B₁, B₂). The restriction is that the bases must have the same parent B-spline basis. That is, they must share the same set of B-spline knots and be of equal polynomial order.

Note that, if both bases are different, the matrix will not be symmetric, and 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).

galerkin_matrix!(A::AbstractMatrix, B::AbstractBSplineBasis,
                 deriv = (Derivative(0), Derivative(0)))

Fill preallocated Galerkin matrix.

The matrix may be a Hermitian view, in which case only one half of the matrix will be filled. Note that, for the matrix to be symmetric, both derivative orders 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 for details.

galerkin_tensor(
    B::AbstractBSplineBasis,
    (D₁::Derivative, D₂::Derivative, D₃::Derivative),
    [T = Float64],
)

Compute 3D banded tensor appearing from quadratic terms in Galerkin method.

As with galerkin_matrix, it is also possible to combine different functional bases by passing, instead of B, a tuple (B₁, B₂, B₃) of three AbstractBSplineBasis. For now, the first two bases, B₁ and B₂, must have the same length.

The tensor is efficiently stored in a BandedTensor3D object.

galerkin_tensor!(
    A::BandedTensor3D,
    B::AbstractBSplineBasis,
    (D₁::Derivative, D₂::Derivative, D₃::Derivative),
)

Compute 3D Galerkin tensor in-place.

See galerkin_tensor for details.