Skip to content

Quadrature options using type dispatch #8

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.

Sign up for GitHub

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

Merged
mrhardman merged 3 commits into from
Jun 17, 2026
Merged

Quadrature options using type dispatch #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6ba74eb
Add abstract type `Abstract1DQuadrature` and function `quadrature1D()…
mrhardman May 22, 2026
9464c83
Extend the available quadrature options to include `GLSpecifiedLimits…
mrhardman May 22, 2026
3f8df15
Include information on how to use `GLSpecifiedLimits()` in the README…
mrhardman Jun 17, 2026
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
33 changes: 33 additions & 0 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,39 @@ K_2D = finite_element_matrix(l_i,l_j,coordinate_x1,
```
where we have two coordinates, `x1` (with polynomials `l_i` and `l_j`) and `x2` (with polynomials `l_k` and `l_n`). Note that the `@enum` variables used to select the Lagrange polynomials are the same between the two coordinates.
The matrix created above is indexed in the order `K_2D[i,k,j,n]`, indexing over the first polynomial `l_i` of `x1`, then the first polynomial `l_k` of `x2`, then the second polynomial `l_j` of `x1` with the final index representing the polynomial `l_n` of `x2`. The same adaptive quadrature keyword arguments are used as in the 1D case.

# Integration limits

For some function inputs $`\rho(x)`$, the domain over which $`\rho(x)`$ is nonzero may not coincide with the element limits $`[x_l, x_u]`$. To handle these cases efficiently with Gauss-Legendre quadrature points, we provide a 1D quadrature option `GLSpecifiedLimits(x_min, x_max)` which can limit the range of integration where necessary without introducing the need to define $`\rho(x)`$ where it is zero, e.g., with a Heaviside function. Note that the values `x_min` and `x_max` should be supplied in the physical coordinate normalisation for the coordinate $`x`$=`coordinate`, rather than in the reference normalisation on the reference domain $`[-1,1]`$ defined for each element. For example, for a $`\rho(x)`$ which is nonzero between $`[x_{-},x_{+}]`$ = `[x_min, x_max]`, we compute the integral
```math
A_{ij} = \int^{\textrm{min}(x_u,x_{+})}_{\textrm{max}(x_l,x_-)} P_i(x) Q_j(x) \rho(x) d x.
```
For the integral
```math
M_{ij} = \int^{\textrm{min}(x_u,x_{+})}_{\textrm{max}(x_l,x_-)} \Phi_i(x) \Phi_j(x) \rho(x) d x,
```
the corresponding syntax for the matrix element $`M_{ij}`$ is as follows.
```
finite_element_matrix(lagrange_x,lagrange_x,coordinate,
kernel_function=(x -> rho(x)),
quadrature_option=GLSpecifiedLimits(x_min,x_max))
```
For matrices of integrals in two dimensions, we use a similar syntax, as follows.
```
finite_element_matrix(lagrange_x,lagrange_x,x1_coordinate,
lagrange_x,lagrange_x,x2_coordinate,
kernel_function=((x1,x2)->rho(x1,x2)),
quadrature_option_x1=GLSpecifiedLimits(x1_min,x1_max),
quadrature_option_x2=GLSpecifiedLimits(x2_min,x2_max))
```
The default quadrature options for integrals in the range $`[x_{l},x_{u}]`$ are as follows.
```
quadrature_option=Default1DQuadrature()
quadrature_option_x1=Default1DQuadrature()
quadrature_option_x2=Default1DQuadrature()
```


# Examples

There are several examples of taking first and second derivatives
Expand Down
79 changes: 68 additions & 11 deletions src/FiniteElementMatrices.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,25 @@ using FastGaussQuadrature: gausslegendre
export lagrange_x,
d_lagrange_dx,
finite_element_matrix,
ElementCoordinates
ElementCoordinates,
GLSpecifiedLimits

@enum LagrangeFunctionType begin
lagrange_x
d_lagrange_dx
end

abstract type Abstract1DQuadrature end
struct Default1DQuadrature <: Abstract1DQuadrature end
struct GLSpecifiedLimits <: Abstract1DQuadrature
# minimum of range of integration
# in physical grid v = s x + c
# with x in [-1, 1]
v_min::Float64
# maximum of range of integration
v_max::Float64
end

struct ElementCoordinates
# precomputed data for calculating
# the Lagrange polynomials, including
Expand Down Expand Up @@ -132,8 +144,9 @@ function _finite_element_matrix(
kernel_function::TFunction=((v -> 1.0)),
additional_quadrature_points::Int64=0,
# argument to permit single implementation of
adaptive_quadrature_points::Int64=0
) where TFunction
adaptive_quadrature_points::Int64=0,
quadrature_option::TQuad=Default1DQuadrature()
) where {TFunction, TQuad <: Abstract1DQuadrature}
lpoly_data = coordinate.lpoly_data
ngrid = length(coordinate.lpoly_data.x_nodes)
scale = coordinate.scale
Expand All @@ -146,7 +159,7 @@ function _finite_element_matrix(
# nquad chosen for exact results for default inputs
# with kernel = 1.0 and zero additional quadrature points
nquad = ngrid + additional_quadrature_points + adaptive_quadrature_points
zz, wz = gausslegendre(nquad)
zz, wz = quadrature1D(quadrature_option,nquad,coordinate)
# compute integral
# int P_i(z) Q_j(z) poly(z) s d z
# with poly(z) = (s z + c)^power
Expand Down Expand Up @@ -193,8 +206,9 @@ function _finite_element_matrix(
# rather than of the reference coordinate z on [-1,1]
kernel_function::TFunction=((v -> 1.0)),
additional_quadrature_points::Int64=0,
adaptive_quadrature_points::Int64=0
) where TFunction
adaptive_quadrature_points::Int64=0,
quadrature_option::TQuad=Default1DQuadrature(),
) where {TFunction, TQuad <: Abstract1DQuadrature}
lpoly_data = coordinate.lpoly_data
ngrid = length(coordinate.lpoly_data.x_nodes)
scale = coordinate.scale
Expand All @@ -208,7 +222,7 @@ function _finite_element_matrix(
# nquad chosen for exact results for default inputs
# with kernel = 1.0 and zero additional quadrature points
nquad = 2*ngrid + additional_quadrature_points + adaptive_quadrature_points
zz, wz = gausslegendre(nquad)
zz, wz = quadrature1D(quadrature_option,nquad,coordinate)
# compute integral
# int P_i(z) Q_j(z) S_k(z) poly(z) d z
# with poly(z) = (s z + c)^power
Expand Down Expand Up @@ -270,8 +284,12 @@ function _finite_element_matrix(
kernel_function::TFunction=((v1,v2) -> 1.0),
additional_quadrature_points_x1::Int64=0,
additional_quadrature_points_x2::Int64=0,
adaptive_quadrature_points::Int64=0
) where TFunction
adaptive_quadrature_points::Int64=0,
quadrature_option_x1::TQuad1=Default1DQuadrature(),
quadrature_option_x2::TQuad2=Default1DQuadrature(),
) where {TFunction,
TQuad1 <: Abstract1DQuadrature,
TQuad2 <: Abstract1DQuadrature}
# coordinate x1 data
lpoly_data_x1 = coordinate_x1.lpoly_data
ngrid_x1 = length(coordinate_x1.lpoly_data.x_nodes)
Expand All @@ -293,9 +311,9 @@ function _finite_element_matrix(
# nquad chosen for exact results for default inputs
# with kernel = 1.0 and zero additional quadrature points
nquad_x1 = ngrid_x1 + additional_quadrature_points_x1 + adaptive_quadrature_points
zz_x1, wz_x1 = gausslegendre(nquad_x1)
zz_x1, wz_x1 = quadrature1D(quadrature_option_x1,nquad_x1,coordinate_x1)
nquad_x2 = ngrid_x2 + additional_quadrature_points_x2 + adaptive_quadrature_points
zz_x2, wz_x2 = gausslegendre(nquad_x2)
zz_x2, wz_x2 = quadrature1D(quadrature_option_x2,nquad_x2,coordinate_x2)
# compute integral
# \int \int \left(P1_i(z_1) Q1_j(z_1) P2_i(z_1) Q2_j(z_2)
# kernel(s_1 z_1 + c_1, s_2 z_2 + c_2) \right) s_1 s_2 d z_1 d z_2
Expand Down Expand Up @@ -333,5 +351,44 @@ function _finite_element_matrix(
return matrix
end

function quadrature1D(::Default1DQuadrature, nquad::Int64, coordinate::ElementCoordinates)
# default quadrature running from [-1,1]
zz, wz = gausslegendre(nquad)
return zz, wz
end
function quadrature1D(qopt::GLSpecifiedLimits, nquad::Int64, coordinate::ElementCoordinates)
# limits specified for integration
v_min = qopt.v_min
v_max = qopt.v_max
# limits based on the range supported by the Lagrange polynomials
v_lower_limit = coordinate.shift - coordinate.scale
v_upper_limit = coordinate.shift + coordinate.scale
# some checks on the inputs
if !(v_min < v_upper_limit)
error("invalid integration range: v_min >= v_upper_limit")
elseif !(v_max > v_lower_limit)
error("invalid integration range: v_max =< v_lower_limit")
elseif !(v_min < v_max)
error("invalid integration range: v_min >= v_max")
end
# obtain valid limits of integration
# in the physical range [v_lower_limit,v_upper_limit]
v_min = max(v_min, v_lower_limit)
v_max = min(v_max, v_upper_limit)
# in the reference range [-1,1]
x_min = (v_min - coordinate.shift)/coordinate.scale
x_max = (v_max - coordinate.shift)/coordinate.scale
# default quadrature running from [-1,1]
zz, wz = gausslegendre(nquad)
# quadrature running from [x_min, x_max]
zx = zeros(Float64, nquad)
wx = zeros(Float64, nquad)
x_scale = 0.5*(x_max - x_min)
x_shift = 0.5*(x_max + x_min)
@. zx = x_scale*zz + x_shift
@. wx = x_scale*wz
return zx, wx
end


end
103 changes: 83 additions & 20 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@ using Test: @testset, @test
using FiniteElementMatrices: lagrange_x,
d_lagrange_dx,
finite_element_matrix,
ElementCoordinates
ElementCoordinates,
GLSpecifiedLimits
using FastGaussQuadrature: gausslobatto, gaussradau, gausslegendre
using LinearAlgebra: lu, ldiv!, mul!

Expand Down Expand Up @@ -52,9 +53,9 @@ function test_first_derivative(;nodes::node_type=GLL,
# check S_ij = P_ji to cover tests where
# derivative acts on test function
@test isapprox(S,transpose(P), atol=2.0e-15)
result = Array{Float64,1}(undef,ngrid)
dummy = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
result = Array{Float64,1}(undef,ngrid)
dummy = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
# function to test
f = Array{Float64,1}(undef,ngrid)
df = Array{Float64,1}(undef,ngrid)
Expand Down Expand Up @@ -96,9 +97,9 @@ function test_first_derivative_nonlinear_operators(
# check Y100_ijk = Y001_kji
@test isapprox(permutedims(Y100, [3,2,1]),Y001,atol=2.0e-15)
# check that a first derivative can be carried out with Y100
result = Array{Float64,1}(undef,ngrid)
result = Array{Float64,1}(undef,ngrid)
dummy = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
# function to test
f = Array{Float64,1}(undef,ngrid)
df = Array{Float64,1}(undef,ngrid)
Expand Down Expand Up @@ -135,34 +136,34 @@ function test_second_derivative(;nodes::node_type=GLL,
M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
S = -finite_element_matrix(d_lagrange_dx,lagrange_x,0,coordinate)
K = -finite_element_matrix(d_lagrange_dx,d_lagrange_dx,0,coordinate)
result = Array{Float64,1}(undef,ngrid)
dummy = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
result = Array{Float64,1}(undef,ngrid)
dummy = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
# function to test
f = Array{Float64,1}(undef,ngrid)
df = Array{Float64,1}(undef,ngrid)
d2f = Array{Float64,1}(undef,ngrid)
# to test the weak first derivative,
# to test the weak first derivative,
# choose a fn that vanishes at x = +-1
# to avoid including boundary terms
for i in 1:ngrid
y = scale*x[i] + shift
f[i] = sin(pi*(y-shift)/scale)
df[i] = (pi/scale)*cos(pi*(y-shift)/scale)
end

# test the performance of a first derivative
# taken using the "weak" methods via
# taken using the "weak" methods via
# integration by parts
luM = lu(M)
mul!(dummy,S,f)
ldiv!(result,luM,dummy)

@. err = result - df
maxerr = maximum(abs.(err))
@test maxerr < atol
# to test the weak second derivative,

# to test the weak second derivative,
# choose a fn that where the derivative vanishes at x = +-1
# to avoid including boundary terms
for i in 1:ngrid
Expand All @@ -171,14 +172,14 @@ function test_second_derivative(;nodes::node_type=GLL,
d2f[i] = -((pi/scale)^2)*cos(pi*(y-shift)/scale)
end
# test the performance of a second derivative
# taken using the "weak" methods via
# taken using the "weak" methods via
# integration by parts
mul!(dummy,K,f)
ldiv!(result,luM,dummy)
@. err = result - d2f
maxerr = maximum(abs.(err))
@test maxerr < atol

end

function test_nonlinear_operators(;nodes::node_type=GLL,
Expand All @@ -195,9 +196,9 @@ function test_nonlinear_operators(;nodes::node_type=GLL,
Y100 = finite_element_matrix(d_lagrange_dx,lagrange_x,lagrange_x,0,coordinate)
Y010 = finite_element_matrix(lagrange_x,d_lagrange_dx,lagrange_x,0,coordinate)
Y001 = finite_element_matrix(lagrange_x,lagrange_x,d_lagrange_dx,0,coordinate)
result = Array{Float64,1}(undef,ngrid)
result = Array{Float64,1}(undef,ngrid)
dummy = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
err = Array{Float64,1}(undef,ngrid)
# function to test
f = Array{Float64,1}(undef,ngrid)
u = Array{Float64,1}(undef,ngrid)
Expand Down Expand Up @@ -330,6 +331,62 @@ function test_nonpolynomial_kernel(;nodes::node_type=GLL,
#println(maxerr)
end

function test_limited_integration_range(;nodes::node_type=GLL,
ngrid::Int64=5,
y0::Float64=-1.1,
y1::Float64=1.3,
atol::Float64=2.0e-13)
x = reference_nodes(nodes,ngrid)
scale = 0.5*(y1-y0)
shift = 0.5*(y0+y1)
coordinate = ElementCoordinates(x,scale,shift)
# midpoint splitting the integral
ymid = y0 + (1.0/3.0)*(y1 - y0)
# test of operators with two polynomials
M_total = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
M_lower_range = finite_element_matrix(lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(y0,ymid))
M_upper_range = finite_element_matrix(lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(ymid,y1))
M_err = deepcopy(M_total)
# test the total matches the exact result
@. M_err = M_total - M_lower_range - M_upper_range
maxerr = maximum(abs.(M_err))
#println(maxerr)
@test maxerr < atol

# test of operators with three polynomials
M_nl_total = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,0,coordinate)
M_nl_lower_range = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(y0,ymid))
M_nl_upper_range = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,coordinate,quadrature_option=GLSpecifiedLimits(ymid,y1))
M_nl_err = deepcopy(M_nl_total)
# test the total matches the exact result
@. M_nl_err = M_nl_total - M_nl_lower_range - M_nl_upper_range
maxerr_nl = maximum(abs.(M_nl_err))
#println(maxerr_nl)
@test maxerr_nl < atol

# test of 2D operators with two polynomials
M2D_total = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate,
lagrange_x,lagrange_x,0,coordinate)
M2D_range1 = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
lagrange_x,lagrange_x,coordinate,
quadrature_option_x1=GLSpecifiedLimits(y0,ymid), quadrature_option_x2=GLSpecifiedLimits(y0,ymid))
M2D_range2 = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
lagrange_x,lagrange_x,coordinate,
quadrature_option_x1=GLSpecifiedLimits(ymid,y1), quadrature_option_x2=GLSpecifiedLimits(y0,ymid))
M2D_range3 = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
lagrange_x,lagrange_x,coordinate,
quadrature_option_x1=GLSpecifiedLimits(y0,ymid), quadrature_option_x2=GLSpecifiedLimits(ymid,y1))
M2D_range4 = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
lagrange_x,lagrange_x,coordinate,
quadrature_option_x1=GLSpecifiedLimits(ymid,y1), quadrature_option_x2=GLSpecifiedLimits(ymid,y1))
M2D_err = deepcopy(M2D_total)
# test the total matches the exact result
@. M2D_err = M2D_total - M2D_range1 - M2D_range2 - M2D_range3 - M2D_range4
maxerr_2D = maximum(abs.(M2D_err))
#println(maxerr_2D)
@test maxerr_2D < atol
end

function test_nonlinear_operator_nonpolynomial_kernel(;nodes::node_type=GLL,
ngrid::Int64=20,
func1::Function=(x -> sin(x)),
Expand Down Expand Up @@ -514,7 +571,7 @@ function runtests()
y0 = -0.5, y1=1.0,
ngrid=n)
end

println("test nonlinear operators:")
for nodes in nodes_list
n = 25
Expand All @@ -541,6 +598,12 @@ function runtests()
test_2D_linear_operators(;nodes=nodes,
ngrid_x1=4, ngrid_x2=5, atol=1.0e-13)
end

println("test limited integration range:")
for nodes in nodes_list
println(" -- test: $(nodes) ngrid=5 poly")
test_limited_integration_range(;nodes = nodes, ngrid=5)
end
end
end

Expand Down
Loading