From 6ba74ebee66b6f6c0d9e6dd6c8691987e9b906be Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 22 May 2026 14:03:32 +0100 Subject: [PATCH 1/3] Add abstract type `Abstract1DQuadrature` and function `quadrature1D()` to create quadrature options chosen by type dispatch. --- src/FiniteElementMatrices.jl | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl index 9d4a9fc..2993930 100644 --- a/src/FiniteElementMatrices.jl +++ b/src/FiniteElementMatrices.jl @@ -18,6 +18,9 @@ export lagrange_x, d_lagrange_dx end +abstract type Abstract1DQuadrature end +struct Default1DQuadrature <: Abstract1DQuadrature end + struct ElementCoordinates # precomputed data for calculating # the Lagrange polynomials, including @@ -132,8 +135,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 @@ -146,7 +150,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) # compute integral # int P_i(z) Q_j(z) poly(z) s d z # with poly(z) = (s z + c)^power @@ -193,8 +197,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 @@ -208,7 +213,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) # compute integral # int P_i(z) Q_j(z) S_k(z) poly(z) d z # with poly(z) = (s z + c)^power @@ -270,8 +275,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) @@ -293,9 +302,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) 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) # 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 @@ -333,5 +342,11 @@ function _finite_element_matrix( return matrix end +function quadrature1D(::Default1DQuadrature, nquad) + # default quadrature running from [-1,1] + zz, wz = gausslegendre(nquad) + return zz, wz +end + end \ No newline at end of file From 9464c8340fad9b58852c026adc06a8ae879b6af1 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 22 May 2026 16:40:00 +0100 Subject: [PATCH 2/3] Extend the available quadrature options to include `GLSpecifiedLimits` where the integration range on the element is permitted to be different from the normal `[-1,1]` (reference coords) `[c - s, c + s]` (physical coordinates). We supply the minimum and maximum limits to be used through `GLSpecifiedLimits(v_min,v_max)` in the physical coordinate `v = s x + c` with `x` the reference coordinate and `s` and `c` the scale and shift factors, respectively. Test this feature by computing the mass matrix for 1D linear and nonlinear operators, and 2D linear operators, using the default method, and by adding together results of the integration over subdomains. --- src/FiniteElementMatrices.jl | 54 ++++++++++++++++-- test/runtests.jl | 103 ++++++++++++++++++++++++++++------- 2 files changed, 131 insertions(+), 26 deletions(-) diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl index 2993930..7c23f4a 100644 --- a/src/FiniteElementMatrices.jl +++ b/src/FiniteElementMatrices.jl @@ -11,7 +11,8 @@ using FastGaussQuadrature: gausslegendre export lagrange_x, d_lagrange_dx, finite_element_matrix, - ElementCoordinates + ElementCoordinates, + GLSpecifiedLimits @enum LagrangeFunctionType begin lagrange_x @@ -20,6 +21,14 @@ 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 @@ -150,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 = quadrature1D(quadrature_option,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 @@ -213,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 = quadrature1D(quadrature_option,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 @@ -302,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 = quadrature1D(quadrature_option_x1,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 = quadrature1D(quadrature_option_x2,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 @@ -342,11 +351,44 @@ function _finite_element_matrix( return matrix end -function quadrature1D(::Default1DQuadrature, nquad) +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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 9299d0a..3f3aec3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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! @@ -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) @@ -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) @@ -135,14 +136,14 @@ 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 @@ -150,19 +151,19 @@ function test_second_derivative(;nodes::node_type=GLL, 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 @@ -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, @@ -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) @@ -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)), @@ -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 @@ -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 From 3f8df15488d49bd8a0f4e2be64d26e847669db04 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 17 Jun 2026 09:24:28 +0100 Subject: [PATCH 3/3] Include information on how to use `GLSpecifiedLimits()` in the README.md. --- README.md | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/README.md b/README.md index c8d3088..0aacaee 100644 --- a/README.md +++ b/README.md @@ -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