From 5a9906bb38e69a7abc0f7fc8ffa55759c7dc03b5 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 24 Nov 2025 14:11:13 +0000
Subject: [PATCH 01/12] Addition of basic 2D finite element matrix that can be
 tested against the product of two 1D matrices. Additional test script
 provided to test this basic functionality.

---
 src/FiniteElementMatrices.jl | 75 ++++++++++++++++++++++++++++++++++++
 test/2DMatricesTest.jl       | 64 ++++++++++++++++++++++++++++++
 test/runtests.jl             |  3 +-
 3 files changed, 141 insertions(+), 1 deletion(-)
 create mode 100644 test/2DMatricesTest.jl

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index e52fb5e..f49e008 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -172,4 +172,79 @@ function finite_element_matrix(
     return matrix
 end
 
+function finite_element_matrix(
+    fn1_x1_type::LagrangeFunctionType,
+    fn2_x1_type::LagrangeFunctionType,
+    power_x1::Int64,
+    coordinate_x1::ElementCoordinates,
+    fn1_x2_type::LagrangeFunctionType,
+    fn2_x2_type::LagrangeFunctionType,
+    power_x2::Int64,
+    coordinate_x2::ElementCoordinates
+    )
+    # coordinate x1 data
+    lpoly_data_x1 = coordinate_x1.lpoly_data
+    ngrid_x1 = length(coordinate_x1.lpoly_data.x_nodes)
+    scale_x1 = coordinate_x1.scale
+    shift_x1 = coordinate_x1.shift
+    # coordinate x2 data
+    lpoly_data_x2 = coordinate_x2.lpoly_data
+    ngrid_x2 = length(coordinate_x2.lpoly_data.x_nodes)
+    scale_x2 = coordinate_x2.scale
+    shift_x2 = coordinate_x2.shift
+    # the finite element array to be returned
+    matrix = zeros(Float64,ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    # the function objects for the required polynomials
+    lagrange11 = select_lagrange_function(fn1_x1_type)
+    lagrange12 = select_lagrange_function(fn1_x2_type)
+    lagrange21 = select_lagrange_function(fn2_x1_type)
+    lagrange22 = select_lagrange_function(fn2_x2_type)
+    # the normalisation factors due to v = scale * x + shift
+    prefactor11 = select_lagrange_prefactor(fn1_x1_type,scale_x1)
+    prefactor12 = select_lagrange_prefactor(fn1_x2_type,scale_x2)
+    prefactor21 = select_lagrange_prefactor(fn2_x1_type,scale_x1)
+    prefactor22 = select_lagrange_prefactor(fn2_x2_type,scale_x2)
+    # nquad chosen for exact results
+    nquad_x1 = ngrid_x1 + power_x1
+    zz_x1, wz_x1 = gausslegendre(nquad_x1)
+    nquad_x2 = ngrid_x2 + power_x2
+    zz_x2, wz_x2 = gausslegendre(nquad_x2)
+    # compute integral
+    # int P_i(z) Q_j(z) poly(z) s d z
+    # with poly(z) = (s z + c)^power
+    # s = scale, c = shift and
+    # P_i(z), Q_i(z) in [l_i(z), (1/s) d l_i(z) / d z]
+    for jx2 in 1:ngrid_x2
+        jx2_lpoly_data = lpoly_data_x2.lpoly_data[jx2]
+        for jx1 in 1:ngrid_x1
+            jx1_lpoly_data = lpoly_data_x1.lpoly_data[jx1]
+            for ix2 in 1:ngrid_x2
+                ix2_lpoly_data = lpoly_data_x2.lpoly_data[ix2]
+                for ix1 in 1:ngrid_x1
+                    ix1_lpoly_data = lpoly_data_x1.lpoly_data[ix1]
+                    for lx2 in 1:nquad_x2
+                        zzl2 = zz_x2[lx2]
+                        polyz2 = get_polyz(power_x2,scale_x2,shift_x2,zzl2)
+                        wgt2 = scale_x2*wz_x2[lx2]
+                        lagrange_factor2 = (prefactor12*lagrange12(ix2_lpoly_data,zzl2)*
+                                               prefactor22*lagrange22(jx2_lpoly_data,zzl2))
+                        for lx1 in 1:nquad_x1
+                            zzl1 = zz_x1[lx1]
+                            polyz1 = get_polyz(power_x1,scale_x1,shift_x1,zzl1)
+                            wgt1 = scale_x1*wz_x1[lx1]
+                            lagrange_factor1 = (prefactor11*lagrange11(ix1_lpoly_data,zzl1)*
+                                               prefactor21*lagrange21(jx1_lpoly_data,zzl1))
+                            matrix[ix1,ix2,jx1,jx2] += (wgt1*polyz1*wgt2*polyz2*
+                                                        lagrange_factor1*
+                                                        lagrange_factor2)
+                        end
+                    end
+                end
+            end
+        end
+    end
+    return matrix
+end
+
+
 end
\ No newline at end of file
diff --git a/test/2DMatricesTest.jl b/test/2DMatricesTest.jl
new file mode 100644
index 0000000..b2733a7
--- /dev/null
+++ b/test/2DMatricesTest.jl
@@ -0,0 +1,64 @@
+using FiniteElementMatrices
+
+using FastGaussQuadrature: gausslobatto, gaussradau, gausslegendre
+using Test
+
+# x1
+# endpoints of grid
+y1 = 1.2
+y0 = -1.3
+ngrid_x1 = 5
+x, w = gausslegendre(ngrid_x1)
+scale_x1 = 0.5*(y1-y0)
+shift_x1 = 0.5*(y0+y1)
+coordinate_x1 = ElementCoordinates(x,scale_x1,shift_x1)
+# x2
+# endpoints of grid
+y1 = 0.9
+y0 = -0.78
+ngrid_x2 = 4
+x, w = gausslegendre(ngrid_x2)
+scale_x2 = 0.5*(y1-y0)
+shift_x2 = 0.5*(y0+y1)
+coordinate_x2 = ElementCoordinates(x,scale_x2,shift_x2)
+
+M_x1 = zeros(ngrid_x1,ngrid_x1)
+M_x2 = zeros(ngrid_x2,ngrid_x2)
+M_2D_err = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+M_2D_exact = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+
+
+@testset "FiniteElementMatrices2Ddev" begin
+    println("test 2D matrices:")
+    for power1 in 1:3
+        for power2 in 1:3
+            for l1A in (lagrange_x, d_lagrange_dx)
+                for l1B in (lagrange_x, d_lagrange_dx)
+                    for l2A in (lagrange_x, d_lagrange_dx)
+                        for l2B in (lagrange_x, d_lagrange_dx)
+                            #1D matrices
+                            M_x1 .= finite_element_matrix(l1A,l1B,power1,coordinate_x1)
+                            M_x2 .= finite_element_matrix(l2A,l2B,power2,coordinate_x2)
+                            #2D matrices
+                            M_2D .= finite_element_matrix(l1A,l1B,power1,coordinate_x1,
+                                                        l2A,l2B,power2,coordinate_x2)
+                            for i2 in 1:ngrid_x2
+                                for i1 in 1:ngrid_x1
+                                    for j2 in 1:ngrid_x2
+                                        for j1 in 1:ngrid_x1
+                                            M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+                                        end
+                                    end
+                                end
+                            end
+                            @. M_2D_err = abs(M_2D-M_2D_exact)
+                            @test maximum(M_2D_err) < 1.0e-13
+                            #println(maximum(M_2D_err))
+                            #println(maximum(M_2D))
+                        end
+                    end
+                end
+            end
+        end
+    end
+end # FiniteElementMatrices2Ddev
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 894d219..2a1a513 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -349,4 +349,5 @@ end
 end
 
 using .FiniteElementMatricesTests
-FiniteElementMatricesTests.runtests()
\ No newline at end of file
+FiniteElementMatricesTests.runtests()
+include(joinpath(@__DIR__,"2DMatricesTest.jl"))
\ No newline at end of file

From e750eaa562d38bb18f84026cd7e1436afe319a32 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 26 Jan 2026 10:52:57 +0000
Subject: [PATCH 02/12] Refactor 2D matrix function for kernel function inputs.

---
 src/FiniteElementMatrices.jl | 35 ++++++++++---
 test/2DMatricesTest.jl       | 96 +++++++++++++++++++-----------------
 2 files changed, 77 insertions(+), 54 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index dda9525..e9d2a80 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -175,7 +175,7 @@ function finite_element_matrix(
                 ith_lpoly_data = lpoly_data.lpoly_data[i]
                 for l in 1:nquad
                     zzl = zz[l]
-                    funcz = kernel_function(scale*zzl+shift)
+                    funcz = kernel_function(scale*zzl+shift)::Float64
                     matrix[i,j,k] += (scale*wz[l]*funcz*
                             prefactor1*lagrange1(ith_lpoly_data,zzl)*
                             prefactor2*lagrange2(jth_lpoly_data,zzl)*
@@ -197,6 +197,25 @@ function finite_element_matrix(
     power_x2::Int64,
     coordinate_x2::ElementCoordinates
     )
+    return finite_element_matrix(fn1_x1_type, fn2_x1_type, coordinate_x1,
+                fn1_x2_type, fn2_x2_type, coordinate_x2;
+                kernel_function=((v1,v2) -> (v1^power_x1)*(v2^power_x2)),
+                additional_quadrature_points_x1=power_x1,
+                additional_quadrature_points_x2=power_x2)
+end
+function finite_element_matrix(
+    fn1_x1_type::LagrangeFunctionType,
+    fn2_x1_type::LagrangeFunctionType,
+    coordinate_x1::ElementCoordinates,
+    fn1_x2_type::LagrangeFunctionType,
+    fn2_x2_type::LagrangeFunctionType,
+    coordinate_x2::ElementCoordinates;
+    # function of the "physical" coord v = s z + c
+    # rather than of the reference coordinate z on [-1,1]
+    kernel_function=((v1,v2) -> 1.0)::Function,
+    additional_quadrature_points_x1=0::Int64,
+    additional_quadrature_points_x2=0::Int64,
+    )
     # coordinate x1 data
     lpoly_data_x1 = coordinate_x1.lpoly_data
     ngrid_x1 = length(coordinate_x1.lpoly_data.x_nodes)
@@ -220,13 +239,13 @@ function finite_element_matrix(
     prefactor21 = select_lagrange_prefactor(fn2_x1_type,scale_x1)
     prefactor22 = select_lagrange_prefactor(fn2_x2_type,scale_x2)
     # nquad chosen for exact results
-    nquad_x1 = ngrid_x1 + power_x1
+    nquad_x1 = ngrid_x1 + additional_quadrature_points_x1
     zz_x1, wz_x1 = gausslegendre(nquad_x1)
-    nquad_x2 = ngrid_x2 + power_x2
+    nquad_x2 = ngrid_x2 + additional_quadrature_points_x2
     zz_x2, wz_x2 = gausslegendre(nquad_x2)
     # compute integral
-    # int P_i(z) Q_j(z) poly(z) s d z
-    # with poly(z) = (s z + c)^power
+    # \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
     # s = scale, c = shift and
     # P_i(z), Q_i(z) in [l_i(z), (1/s) d l_i(z) / d z]
     for jx2 in 1:ngrid_x2
@@ -239,17 +258,17 @@ function finite_element_matrix(
                     ix1_lpoly_data = lpoly_data_x1.lpoly_data[ix1]
                     for lx2 in 1:nquad_x2
                         zzl2 = zz_x2[lx2]
-                        polyz2 = get_polyz(power_x2,scale_x2,shift_x2,zzl2)
                         wgt2 = scale_x2*wz_x2[lx2]
                         lagrange_factor2 = (prefactor12*lagrange12(ix2_lpoly_data,zzl2)*
                                                prefactor22*lagrange22(jx2_lpoly_data,zzl2))
                         for lx1 in 1:nquad_x1
                             zzl1 = zz_x1[lx1]
-                            polyz1 = get_polyz(power_x1,scale_x1,shift_x1,zzl1)
                             wgt1 = scale_x1*wz_x1[lx1]
                             lagrange_factor1 = (prefactor11*lagrange11(ix1_lpoly_data,zzl1)*
                                                prefactor21*lagrange21(jx1_lpoly_data,zzl1))
-                            matrix[ix1,ix2,jx1,jx2] += (wgt1*polyz1*wgt2*polyz2*
+                            kernel_x1_x2 = kernel_function(scale_x1*zzl1+shift_x1,
+                                                        scale_x2*zzl2+shift_x2)::Float64
+                            matrix[ix1,ix2,jx1,jx2] += (wgt1*wgt2*kernel_x1_x2*
                                                         lagrange_factor1*
                                                         lagrange_factor2)
                         end
diff --git a/test/2DMatricesTest.jl b/test/2DMatricesTest.jl
index b2733a7..78d6480 100644
--- a/test/2DMatricesTest.jl
+++ b/test/2DMatricesTest.jl
@@ -3,62 +3,66 @@ using FiniteElementMatrices
 using FastGaussQuadrature: gausslobatto, gaussradau, gausslegendre
 using Test
 
-# x1
-# endpoints of grid
-y1 = 1.2
-y0 = -1.3
-ngrid_x1 = 5
-x, w = gausslegendre(ngrid_x1)
-scale_x1 = 0.5*(y1-y0)
-shift_x1 = 0.5*(y0+y1)
-coordinate_x1 = ElementCoordinates(x,scale_x1,shift_x1)
-# x2
-# endpoints of grid
-y1 = 0.9
-y0 = -0.78
-ngrid_x2 = 4
-x, w = gausslegendre(ngrid_x2)
-scale_x2 = 0.5*(y1-y0)
-shift_x2 = 0.5*(y0+y1)
-coordinate_x2 = ElementCoordinates(x,scale_x2,shift_x2)
+function runtests_2D(;ngrid_x1=5,ngrid_x2=4)
+    # x1
+    # endpoints of grid
+    y1 = 1.2
+    y0 = -1.3
+    x1, w1 = gausslegendre(ngrid_x1)
+    scale_x1 = 0.5*(y1-y0)
+    shift_x1 = 0.5*(y0+y1)
+    coordinate_x1 = ElementCoordinates(x1,scale_x1,shift_x1)
+    # x2
+    # endpoints of grid
+    y1 = 0.9
+    y0 = -0.78
+    x2, w2 = gausslegendre(ngrid_x2)
+    scale_x2 = 0.5*(y1-y0)
+    shift_x2 = 0.5*(y0+y1)
+    coordinate_x2 = ElementCoordinates(x2,scale_x2,shift_x2)
 
-M_x1 = zeros(ngrid_x1,ngrid_x1)
-M_x2 = zeros(ngrid_x2,ngrid_x2)
-M_2D_err = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
-M_2D_exact = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_x1 = zeros(ngrid_x1,ngrid_x1)
+    M_x2 = zeros(ngrid_x2,ngrid_x2)
+    M_2D = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_2D_err = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_2D_exact = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
 
 
-@testset "FiniteElementMatrices2Ddev" begin
-    println("test 2D matrices:")
-    for power1 in 1:3
-        for power2 in 1:3
-            for l1A in (lagrange_x, d_lagrange_dx)
-                for l1B in (lagrange_x, d_lagrange_dx)
-                    for l2A in (lagrange_x, d_lagrange_dx)
-                        for l2B in (lagrange_x, d_lagrange_dx)
-                            #1D matrices
-                            M_x1 .= finite_element_matrix(l1A,l1B,power1,coordinate_x1)
-                            M_x2 .= finite_element_matrix(l2A,l2B,power2,coordinate_x2)
-                            #2D matrices
-                            M_2D .= finite_element_matrix(l1A,l1B,power1,coordinate_x1,
-                                                        l2A,l2B,power2,coordinate_x2)
-                            for i2 in 1:ngrid_x2
-                                for i1 in 1:ngrid_x1
-                                    for j2 in 1:ngrid_x2
-                                        for j1 in 1:ngrid_x1
-                                            M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+    @testset "FiniteElementMatrices2Ddev" begin
+        println("test 2D matrices:")
+        for power1 in 0:3
+            for power2 in 0:3
+                for l1A in (lagrange_x, d_lagrange_dx)
+                    for l1B in (lagrange_x, d_lagrange_dx)
+                        for l2A in (lagrange_x, d_lagrange_dx)
+                            for l2B in (lagrange_x, d_lagrange_dx)
+                                #1D matrices
+                                M_x1 .= finite_element_matrix(l1A,l1B,power1,coordinate_x1)
+                                M_x2 .= finite_element_matrix(l2A,l2B,power2,coordinate_x2)
+                                #2D matrices
+                                M_2D .= finite_element_matrix(l1A,l1B,power1,coordinate_x1,
+                                                            l2A,l2B,power2,coordinate_x2)
+                                for i2 in 1:ngrid_x2
+                                    for i1 in 1:ngrid_x1
+                                        for j2 in 1:ngrid_x2
+                                            for j1 in 1:ngrid_x1
+                                                M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+                                            end
                                         end
                                     end
                                 end
+                                @. M_2D_err = abs(M_2D-M_2D_exact)
+                                @test maximum(M_2D_err) < 1.0e-13
+                                #println(maximum(M_2D_err))
+                                #println(maximum(M_2D))
                             end
-                            @. M_2D_err = abs(M_2D-M_2D_exact)
-                            @test maximum(M_2D_err) < 1.0e-13
-                            #println(maximum(M_2D_err))
-                            #println(maximum(M_2D))
                         end
                     end
                 end
             end
         end
     end
-end # FiniteElementMatrices2Ddev
\ No newline at end of file
+    return nothing
+end
+runtests_2D(ngrid_x1=5,ngrid_x2=4)
+runtests_2D(ngrid_x1=4,ngrid_x2=5)

From 30bd714bff067df938e7dacede2f6918408b9452 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 26 Jan 2026 12:56:15 +0000
Subject: [PATCH 03/12] Refactor to remove the separate function that handles
 the Lagrange polynomial prefactor.

---
 src/FiniteElementMatrices.jl | 68 +++++++++++++-----------------------
 1 file changed, 25 insertions(+), 43 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index e9d2a80..5903e6e 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -4,7 +4,8 @@ module FiniteElementMatrices
 
 using LagrangePolynomials: lagrange_poly,
                            lagrange_poly_derivative,
-                           LagrangePolyData
+                           LagrangePolyData,
+                           LagrangePolyDataBase
 using FastGaussQuadrature: gausslegendre
 
 export lagrange_x,
@@ -49,25 +50,17 @@ struct ElementCoordinates
     end
 end
 
-function select_lagrange_function(fn_type::LagrangeFunctionType)
+function select_lagrange_function(fn_type::LagrangeFunctionType,scale::Float64)
+    # including the normalisation factors due to v = scale * x + shift
     if fn_type == lagrange_x
-        func = lagrange_poly
+        func = ((lpoly_data::LagrangePolyDataBase,z::Float64) -> lagrange_poly(lpoly_data,z))
     elseif fn_type == d_lagrange_dx
-        func = lagrange_poly_derivative
+        # 1/scale due to scale in x in d l(z(x)) /d x
+        func = ((lpoly_data::LagrangePolyDataBase,z::Float64) -> (1.0/scale)*lagrange_poly_derivative(lpoly_data,z))
     end
     return func
 end
 
-function select_lagrange_prefactor(fn_type::LagrangeFunctionType,
-                                    scale::Float64)
-    if fn_type == lagrange_x
-        prefactor = 1.0
-    elseif fn_type == d_lagrange_dx
-        prefactor = 1.0/scale
-    end
-    return prefactor
-end
-
 function finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
@@ -95,11 +88,8 @@ function finite_element_matrix(
     # the finite element array to be returned
     matrix = zeros(Float64,ngrid,ngrid)
     # the function objects for the required polynomials
-    lagrange1 = select_lagrange_function(fn1_type)
-    lagrange2 = select_lagrange_function(fn2_type)
-    # the normalisation factors due to v = scale * x + shift
-    prefactor1 = select_lagrange_prefactor(fn1_type,scale)
-    prefactor2 = select_lagrange_prefactor(fn2_type,scale)
+    lagrange1 = select_lagrange_function(fn1_type,scale)
+    lagrange2 = select_lagrange_function(fn2_type,scale)
     # nquad chosen for exact results for all power, ngrid
     nquad = ngrid + additional_quadrature_points
     zz, wz = gausslegendre(nquad)
@@ -116,8 +106,8 @@ function finite_element_matrix(
                 zzl = zz[l]
                 funcz = kernel_function(scale*zzl+shift)::Float64
                 matrix[i,j] += (scale*wz[l]*funcz*
-                           prefactor1*lagrange1(ith_lpoly_data,zzl)*
-                           prefactor2*lagrange2(jth_lpoly_data,zzl))
+                           lagrange1(ith_lpoly_data,zzl)*
+                           lagrange2(jth_lpoly_data,zzl))
             end
         end
     end
@@ -152,13 +142,9 @@ function finite_element_matrix(
     # the finite element array to be returned
     matrix = zeros(Float64,ngrid,ngrid,ngrid)
     # the function objects for the required polynomials
-    lagrange1 = select_lagrange_function(fn1_type)
-    lagrange2 = select_lagrange_function(fn2_type)
-    lagrange3 = select_lagrange_function(fn3_type)
-    # the normalisation factors due to v = scale * x + shift
-    prefactor1 = select_lagrange_prefactor(fn1_type,scale)
-    prefactor2 = select_lagrange_prefactor(fn2_type,scale)
-    prefactor3 = select_lagrange_prefactor(fn3_type,scale)
+    lagrange1 = select_lagrange_function(fn1_type,scale)
+    lagrange2 = select_lagrange_function(fn2_type,scale)
+    lagrange3 = select_lagrange_function(fn3_type,scale)
     # nquad chosen for exact results for all power, ngrid
     nquad = 2*ngrid + additional_quadrature_points
     zz, wz = gausslegendre(nquad)
@@ -177,9 +163,9 @@ function finite_element_matrix(
                     zzl = zz[l]
                     funcz = kernel_function(scale*zzl+shift)::Float64
                     matrix[i,j,k] += (scale*wz[l]*funcz*
-                            prefactor1*lagrange1(ith_lpoly_data,zzl)*
-                            prefactor2*lagrange2(jth_lpoly_data,zzl)*
-                            prefactor3*lagrange3(kth_lpoly_data,zzl))
+                            lagrange1(ith_lpoly_data,zzl)*
+                            lagrange2(jth_lpoly_data,zzl)*
+                            lagrange3(kth_lpoly_data,zzl))
                 end
             end
         end
@@ -229,15 +215,11 @@ function finite_element_matrix(
     # the finite element array to be returned
     matrix = zeros(Float64,ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
     # the function objects for the required polynomials
-    lagrange11 = select_lagrange_function(fn1_x1_type)
-    lagrange12 = select_lagrange_function(fn1_x2_type)
-    lagrange21 = select_lagrange_function(fn2_x1_type)
-    lagrange22 = select_lagrange_function(fn2_x2_type)
+    lagrange11 = select_lagrange_function(fn1_x1_type,scale_x1)
+    lagrange12 = select_lagrange_function(fn1_x2_type,scale_x2)
+    lagrange21 = select_lagrange_function(fn2_x1_type,scale_x1)
+    lagrange22 = select_lagrange_function(fn2_x2_type,scale_x2)
     # the normalisation factors due to v = scale * x + shift
-    prefactor11 = select_lagrange_prefactor(fn1_x1_type,scale_x1)
-    prefactor12 = select_lagrange_prefactor(fn1_x2_type,scale_x2)
-    prefactor21 = select_lagrange_prefactor(fn2_x1_type,scale_x1)
-    prefactor22 = select_lagrange_prefactor(fn2_x2_type,scale_x2)
     # nquad chosen for exact results
     nquad_x1 = ngrid_x1 + additional_quadrature_points_x1
     zz_x1, wz_x1 = gausslegendre(nquad_x1)
@@ -259,13 +241,13 @@ function finite_element_matrix(
                     for lx2 in 1:nquad_x2
                         zzl2 = zz_x2[lx2]
                         wgt2 = scale_x2*wz_x2[lx2]
-                        lagrange_factor2 = (prefactor12*lagrange12(ix2_lpoly_data,zzl2)*
-                                               prefactor22*lagrange22(jx2_lpoly_data,zzl2))
+                        lagrange_factor2 = (lagrange12(ix2_lpoly_data,zzl2)*
+                                               lagrange22(jx2_lpoly_data,zzl2))
                         for lx1 in 1:nquad_x1
                             zzl1 = zz_x1[lx1]
                             wgt1 = scale_x1*wz_x1[lx1]
-                            lagrange_factor1 = (prefactor11*lagrange11(ix1_lpoly_data,zzl1)*
-                                               prefactor21*lagrange21(jx1_lpoly_data,zzl1))
+                            lagrange_factor1 = (lagrange11(ix1_lpoly_data,zzl1)*
+                                               lagrange21(jx1_lpoly_data,zzl1))
                             kernel_x1_x2 = kernel_function(scale_x1*zzl1+shift_x1,
                                                         scale_x2*zzl2+shift_x2)::Float64
                             matrix[ix1,ix2,jx1,jx2] += (wgt1*wgt2*kernel_x1_x2*

From 150bea7e5a61831be03d53f3e6ee119b3f525c5b Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 26 Jan 2026 16:05:18 +0000
Subject: [PATCH 04/12] Move 2D matrices test into main test/runtests.jl file.

---
 test/2DMatricesTest.jl | 68 -------------------------------------
 test/runtests.jl       | 76 +++++++++++++++++++++++++++++++++++++++---
 2 files changed, 71 insertions(+), 73 deletions(-)
 delete mode 100644 test/2DMatricesTest.jl

diff --git a/test/2DMatricesTest.jl b/test/2DMatricesTest.jl
deleted file mode 100644
index 78d6480..0000000
--- a/test/2DMatricesTest.jl
+++ /dev/null
@@ -1,68 +0,0 @@
-using FiniteElementMatrices
-
-using FastGaussQuadrature: gausslobatto, gaussradau, gausslegendre
-using Test
-
-function runtests_2D(;ngrid_x1=5,ngrid_x2=4)
-    # x1
-    # endpoints of grid
-    y1 = 1.2
-    y0 = -1.3
-    x1, w1 = gausslegendre(ngrid_x1)
-    scale_x1 = 0.5*(y1-y0)
-    shift_x1 = 0.5*(y0+y1)
-    coordinate_x1 = ElementCoordinates(x1,scale_x1,shift_x1)
-    # x2
-    # endpoints of grid
-    y1 = 0.9
-    y0 = -0.78
-    x2, w2 = gausslegendre(ngrid_x2)
-    scale_x2 = 0.5*(y1-y0)
-    shift_x2 = 0.5*(y0+y1)
-    coordinate_x2 = ElementCoordinates(x2,scale_x2,shift_x2)
-
-    M_x1 = zeros(ngrid_x1,ngrid_x1)
-    M_x2 = zeros(ngrid_x2,ngrid_x2)
-    M_2D = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
-    M_2D_err = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
-    M_2D_exact = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
-
-
-    @testset "FiniteElementMatrices2Ddev" begin
-        println("test 2D matrices:")
-        for power1 in 0:3
-            for power2 in 0:3
-                for l1A in (lagrange_x, d_lagrange_dx)
-                    for l1B in (lagrange_x, d_lagrange_dx)
-                        for l2A in (lagrange_x, d_lagrange_dx)
-                            for l2B in (lagrange_x, d_lagrange_dx)
-                                #1D matrices
-                                M_x1 .= finite_element_matrix(l1A,l1B,power1,coordinate_x1)
-                                M_x2 .= finite_element_matrix(l2A,l2B,power2,coordinate_x2)
-                                #2D matrices
-                                M_2D .= finite_element_matrix(l1A,l1B,power1,coordinate_x1,
-                                                            l2A,l2B,power2,coordinate_x2)
-                                for i2 in 1:ngrid_x2
-                                    for i1 in 1:ngrid_x1
-                                        for j2 in 1:ngrid_x2
-                                            for j1 in 1:ngrid_x1
-                                                M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
-                                            end
-                                        end
-                                    end
-                                end
-                                @. M_2D_err = abs(M_2D-M_2D_exact)
-                                @test maximum(M_2D_err) < 1.0e-13
-                                #println(maximum(M_2D_err))
-                                #println(maximum(M_2D))
-                            end
-                        end
-                    end
-                end
-            end
-        end
-    end
-    return nothing
-end
-runtests_2D(ngrid_x1=5,ngrid_x2=4)
-runtests_2D(ngrid_x1=4,ngrid_x2=5)
diff --git a/test/runtests.jl b/test/runtests.jl
index fee42a5..c239592 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -379,6 +379,66 @@ function test_nonlinear_operator_nonpolynomial_kernel(;nodes::node_type=GLL,
     #println(maxerr)
 end
 
+function test_2D_linear_operators(;nodes::node_type=GLL,
+            ngrid_x1=5::Int64,ngrid_x2=4::Int64,
+            atol=1.0e-13::Float64)
+    # x1
+    # endpoints of grid
+    y1 = 1.2
+    y0 = -1.3
+    x1 = reference_nodes(nodes,ngrid_x1)
+    scale_x1 = 0.5*(y1-y0)
+    shift_x1 = 0.5*(y0+y1)
+    coordinate_x1 = ElementCoordinates(x1,scale_x1,shift_x1)
+    # x2
+    # endpoints of grid
+    y1 = 0.9
+    y0 = -0.78
+    x2 = reference_nodes(nodes,ngrid_x2)
+    scale_x2 = 0.5*(y1-y0)
+    shift_x2 = 0.5*(y0+y1)
+    coordinate_x2 = ElementCoordinates(x2,scale_x2,shift_x2)
+
+    M_x1 = zeros(ngrid_x1,ngrid_x1)
+    M_x2 = zeros(ngrid_x2,ngrid_x2)
+    M_2D = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_2D_err = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    M_2D_exact = zeros(ngrid_x1,ngrid_x2,ngrid_x1,ngrid_x2)
+    # test that a 2D matrix with polynomial kernel can be
+    # calculated with an identical result to the
+    # product of 1D matrices with the same kernels
+    for power1 in 0:3
+        for power2 in 0:3
+            for l1A in (lagrange_x, d_lagrange_dx)
+                for l1B in (lagrange_x, d_lagrange_dx)
+                    for l2A in (lagrange_x, d_lagrange_dx)
+                        for l2B in (lagrange_x, d_lagrange_dx)
+                            #1D matrices
+                            M_x1 .= finite_element_matrix(l1A,l1B,power1,coordinate_x1)
+                            M_x2 .= finite_element_matrix(l2A,l2B,power2,coordinate_x2)
+                            #2D matrices
+                            M_2D .= finite_element_matrix(l1A,l1B,power1,coordinate_x1,
+                                                        l2A,l2B,power2,coordinate_x2)
+                            for i2 in 1:ngrid_x2
+                                for i1 in 1:ngrid_x1
+                                    for j2 in 1:ngrid_x2
+                                        for j1 in 1:ngrid_x1
+                                            M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+                                        end
+                                    end
+                                end
+                            end
+                            @. M_2D_err = abs(M_2D-M_2D_exact)
+                            @test maximum(M_2D_err) < atol
+                        end
+                    end
+                end
+            end
+        end
+    end
+    return nothing
+end
+
 function runtests()
     @testset "FiniteElementMatrices" begin
         println("test FiniteElementMatrices")
@@ -414,7 +474,6 @@ function runtests()
             end
         end
         println("test second derivative:")
-        nodes_list = [GLL, GLR, GLe]
         for nodes in nodes_list
             n = 20
             println("    -- test: $(nodes) ngrid=$(n) trig")
@@ -424,7 +483,6 @@ function runtests()
         end
         
         println("test nonlinear operators:")
-        nodes_list = [GLL, GLR, GLe]
         for nodes in nodes_list
             n = 25
             println("    -- test: $(nodes) ngrid=$(n) trig")
@@ -434,18 +492,26 @@ function runtests()
         end
 
         println("test nonpolynomial kernels:")
-        nodes_list = [GLL, GLR, GLe]
         for nodes in nodes_list
             n = 20
             println("    -- test: $(nodes) ngrid=$(n) trig")
             test_nonpolynomial_kernel(; nodes=nodes, ngrid=n)
             test_nonlinear_operator_nonpolynomial_kernel(; nodes=nodes, ngrid=n)
         end
+
+        println("test 2D operators:")
+        for nodes in nodes_list
+            println("    -- test: $(nodes) ngrid_x1=5 ngrid_x2=4 poly")
+            test_2D_linear_operators(;nodes=nodes,
+                ngrid_x1=5, ngrid_x2=4, atol=1.0e-13)
+            println("    -- test: $(nodes) ngrid_x1=4 ngrid_x2=5 poly")
+            test_2D_linear_operators(;nodes=nodes,
+                ngrid_x1=4, ngrid_x2=5, atol=1.0e-13)
+        end
     end
 end
 
 end
 
 using .FiniteElementMatricesTests
-FiniteElementMatricesTests.runtests()
-include(joinpath(@__DIR__,"2DMatricesTest.jl"))
\ No newline at end of file
+FiniteElementMatricesTests.runtests()
\ No newline at end of file

From 36b0fcc72da333c1c777fad79ee75782e15c3009 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 28 Jan 2026 08:43:04 +0000
Subject: [PATCH 05/12] Update comments.

---
 src/FiniteElementMatrices.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index 5903e6e..4e227ed 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -90,7 +90,8 @@ function finite_element_matrix(
     # the function objects for the required polynomials
     lagrange1 = select_lagrange_function(fn1_type,scale)
     lagrange2 = select_lagrange_function(fn2_type,scale)
-    # nquad chosen for exact results for all power, ngrid
+    # nquad chosen for exact results for default inputs
+    # with kernel = 1.0 and zero additional quadrature points
     nquad = ngrid + additional_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
@@ -145,7 +146,8 @@ function finite_element_matrix(
     lagrange1 = select_lagrange_function(fn1_type,scale)
     lagrange2 = select_lagrange_function(fn2_type,scale)
     lagrange3 = select_lagrange_function(fn3_type,scale)
-    # nquad chosen for exact results for all power, ngrid
+    # nquad chosen for exact results for default inputs
+    # with kernel = 1.0 and zero additional quadrature points
     nquad = 2*ngrid + additional_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
@@ -220,7 +222,8 @@ function finite_element_matrix(
     lagrange21 = select_lagrange_function(fn2_x1_type,scale_x1)
     lagrange22 = select_lagrange_function(fn2_x2_type,scale_x2)
     # the normalisation factors due to v = scale * x + shift
-    # nquad chosen for exact results
+    # 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
     zz_x1, wz_x1 = gausslegendre(nquad_x1)
     nquad_x2 = ngrid_x2 + additional_quadrature_points_x2

From ded1a1f8e12eb48c48b37ccba8e301ce9a1add66 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 28 Jan 2026 10:58:29 +0000
Subject: [PATCH 06/12] Add p-adaptive quadrature method to permit the user to
 supply complicated kernels and quickly obtain accurate results.

---
 src/FiniteElementMatrices.jl | 69 ++++++++++++++++++++++++++++++------
 1 file changed, 59 insertions(+), 10 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index 4e227ed..5022548 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -61,18 +61,53 @@ function select_lagrange_function(fn_type::LagrangeFunctionType,scale::Float64)
     return func
 end
 
+# Method for p-adaptive quadrature in 1D and 2D.
+# Note that the the rate at which the number of quadrature
+# points is increased is the same for each coordinate.
+function finite_element_matrix(args...;
+    atol=1.0e-12::Float64, rtol=1.0e-12::Float64, max_iterations=1::Int64,
+    quadrature_increment=5::Int64, kwargs...)
+    matrix = _finite_element_matrix(args...;
+                adaptive_quadrature_points=0,
+                kwargs...)
+    iteration = 0
+    maxvalue = maximum(abs.(matrix))
+    maxerr = maxvalue
+    while iteration < max_iterations && maxerr > atol + rtol*maxvalue
+        # increment the iteration
+        iteration += 1
+        matrixnew = _finite_element_matrix(args...;
+                        adaptive_quadrature_points=quadrature_increment*iteration,
+                        kwargs...)
+        maxvalue = maximum(abs.(matrixnew))
+        maxerr = maximum(abs.(matrixnew .- matrix))
+        # update the matrix
+        matrix .= matrixnew
+    end
+    if maxerr > atol + rtol*maxvalue && max_iterations > 0
+        error("Unable to achieve specified error bounds: \n
+              maxerr < atol + rtol*maxvalue \n
+              after max_iterations = $max_iterations with atol=$atol and rtol=$rtol we obtain \n
+              maxerr = $maxerr > $atol + $rtol*$maxvalue")
+    end
+    return matrix
+end
+
+# Method for computing polynomial integrals exactly
+# for matrices with two Lagrange polymomials in 1D
 function finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     power::Int64,
     coordinate::ElementCoordinates
     )
-    return finite_element_matrix(fn1_type, fn2_type, coordinate;
+    return _finite_element_matrix(fn1_type, fn2_type, coordinate;
                 kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
 end
 
-function finite_element_matrix(
+# The basic method for matrices with two Lagrange polymomials in 1D
+function _finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     coordinate::ElementCoordinates;
@@ -80,6 +115,8 @@ function finite_element_matrix(
     # rather than of the reference coordinate z on [-1,1]
     kernel_function=((v -> 1.0))::Function,
     additional_quadrature_points=0::Int64,
+    # argument to permit single implementation of
+    adaptive_quadrature_points=0::Int64
     )
     lpoly_data = coordinate.lpoly_data
     ngrid = length(coordinate.lpoly_data.x_nodes)
@@ -92,7 +129,7 @@ function finite_element_matrix(
     lagrange2 = select_lagrange_function(fn2_type,scale)
     # nquad chosen for exact results for default inputs
     # with kernel = 1.0 and zero additional quadrature points
-    nquad = ngrid + additional_quadrature_points
+    nquad = ngrid + additional_quadrature_points + adaptive_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
     # int P_i(z) Q_j(z) poly(z) s d z
@@ -115,6 +152,8 @@ function finite_element_matrix(
     return matrix
 end
 
+# Method for computing polynomial integrals exactly
+# for matrices with three Lagrange polymomials in 1D
 function finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
@@ -122,11 +161,13 @@ function finite_element_matrix(
     power::Int64,
     coordinate::ElementCoordinates
     )
-    return finite_element_matrix(fn1_type, fn2_type, fn3_type,
+    return _finite_element_matrix(fn1_type, fn2_type, fn3_type,
                 coordinate; kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
 end
-function finite_element_matrix(
+
+# The basic method for matrices with three Lagrange polymomials in 2D
+function _finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     fn3_type::LagrangeFunctionType,
@@ -135,6 +176,7 @@ function finite_element_matrix(
     # rather than of the reference coordinate z on [-1,1]
     kernel_function=((v -> 1.0))::Function,
     additional_quadrature_points=0::Int64,
+    adaptive_quadrature_points=0::Int64
     )
     lpoly_data = coordinate.lpoly_data
     ngrid = length(coordinate.lpoly_data.x_nodes)
@@ -148,7 +190,7 @@ function finite_element_matrix(
     lagrange3 = select_lagrange_function(fn3_type,scale)
     # nquad chosen for exact results for default inputs
     # with kernel = 1.0 and zero additional quadrature points
-    nquad = 2*ngrid + additional_quadrature_points
+    nquad = 2*ngrid + additional_quadrature_points + adaptive_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
     # int P_i(z) Q_j(z) S_k(z) poly(z) d z
@@ -175,6 +217,9 @@ function finite_element_matrix(
     return matrix
 end
 
+# Method for computing polynomial integrals exactly
+# for matrices with two 2D basis functions
+# (two Lagrange polymomials per coordinate) in 2D
 function finite_element_matrix(
     fn1_x1_type::LagrangeFunctionType,
     fn2_x1_type::LagrangeFunctionType,
@@ -185,13 +230,16 @@ function finite_element_matrix(
     power_x2::Int64,
     coordinate_x2::ElementCoordinates
     )
-    return finite_element_matrix(fn1_x1_type, fn2_x1_type, coordinate_x1,
+    return _finite_element_matrix(fn1_x1_type, fn2_x1_type, coordinate_x1,
                 fn1_x2_type, fn2_x2_type, coordinate_x2;
                 kernel_function=((v1,v2) -> (v1^power_x1)*(v2^power_x2)),
                 additional_quadrature_points_x1=power_x1,
                 additional_quadrature_points_x2=power_x2)
 end
-function finite_element_matrix(
+
+# Basic method for computing matrices with two 2D basis functions
+# (two Lagrange polymomials per coordinate) in 2D
+function _finite_element_matrix(
     fn1_x1_type::LagrangeFunctionType,
     fn2_x1_type::LagrangeFunctionType,
     coordinate_x1::ElementCoordinates,
@@ -203,6 +251,7 @@ function finite_element_matrix(
     kernel_function=((v1,v2) -> 1.0)::Function,
     additional_quadrature_points_x1=0::Int64,
     additional_quadrature_points_x2=0::Int64,
+    adaptive_quadrature_points=0::Int64
     )
     # coordinate x1 data
     lpoly_data_x1 = coordinate_x1.lpoly_data
@@ -224,9 +273,9 @@ function finite_element_matrix(
     # the normalisation factors due to v = scale * x + shift
     # 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
+    nquad_x1 = ngrid_x1 + additional_quadrature_points_x1 + adaptive_quadrature_points
     zz_x1, wz_x1 = gausslegendre(nquad_x1)
-    nquad_x2 = ngrid_x2 + additional_quadrature_points_x2
+    nquad_x2 = ngrid_x2 + additional_quadrature_points_x2 + adaptive_quadrature_points
     zz_x2, wz_x2 = gausslegendre(nquad_x2)
     # compute integral
     # \int \int \left(P1_i(z_1) Q1_j(z_1) P2_i(z_1) Q2_j(z_2)

From f5ef90934f854cd231594f51e439df13e466fbe4 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 28 Jan 2026 11:08:32 +0000
Subject: [PATCH 07/12] Add a check to prevent calculation of matrices with
 power < 0.

---
 src/FiniteElementMatrices.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index 5022548..304fad3 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -61,6 +61,12 @@ function select_lagrange_function(fn_type::LagrangeFunctionType,scale::Float64)
     return func
 end
 
+function check_power(power::Int64, varname::String)
+    if power < 0
+        error("Require power > 0, whereas $varname=$power was passed")
+    end
+    return nothing
+end
 # Method for p-adaptive quadrature in 1D and 2D.
 # Note that the the rate at which the number of quadrature
 # points is increased is the same for each coordinate.
@@ -101,6 +107,7 @@ function finite_element_matrix(
     power::Int64,
     coordinate::ElementCoordinates
     )
+    check_power(power, "power")
     return _finite_element_matrix(fn1_type, fn2_type, coordinate;
                 kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
@@ -161,6 +168,7 @@ function finite_element_matrix(
     power::Int64,
     coordinate::ElementCoordinates
     )
+    check_power(power, "power")
     return _finite_element_matrix(fn1_type, fn2_type, fn3_type,
                 coordinate; kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
@@ -230,6 +238,8 @@ function finite_element_matrix(
     power_x2::Int64,
     coordinate_x2::ElementCoordinates
     )
+    check_power(power_x1, "power_x1")
+    check_power(power_x2, "power_x2")
     return _finite_element_matrix(fn1_x1_type, fn2_x1_type, coordinate_x1,
                 fn1_x2_type, fn2_x2_type, coordinate_x2;
                 kernel_function=((v1,v2) -> (v1^power_x1)*(v2^power_x2)),

From 51f57389782fea56b912dff2d34a70371832f270 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 28 Jan 2026 11:25:06 +0000
Subject: [PATCH 08/12] Add a test of the p-adaptive interface for 2D matrices,
 comparing the result to a produce of 1D matrices.

---
 test/runtests.jl | 33 +++++++++++++++++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index c239592..b7d36e4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -436,6 +436,39 @@ function test_2D_linear_operators(;nodes::node_type=GLL,
             end
         end
     end
+    # test produces of 1D matrices for a simple separable Jacobian
+    for l1A in (lagrange_x, d_lagrange_dx)
+        for l1B in (lagrange_x, d_lagrange_dx)
+            for l2A in (lagrange_x, d_lagrange_dx)
+                for l2B in (lagrange_x, d_lagrange_dx)
+                    #1D matrices
+                    M_x1 .= finite_element_matrix(l1A,l1B,coordinate_x1;
+                                kernel_function=(x1 -> exp(-x1^2)), quadrature_increment=10,
+                                max_iterations=10, atol=1.0e-13, rtol=1.0e-13)
+                    M_x2 .= finite_element_matrix(l2A,l2B,coordinate_x2;
+                                kernel_function=(x2 -> exp(-x2^2)), quadrature_increment=10,
+                                max_iterations=10, atol=1.0e-13, rtol=1.0e-13)
+                    #2D matrices
+                    M_2D .= finite_element_matrix(l1A,l1B,coordinate_x1,
+                                                l2A,l2B,coordinate_x2;
+                                                kernel_function=((x1,x2)->exp(-x1^2 - x2^2)),
+                                                max_iterations=10, quadrature_increment=10,
+                                                atol=1.0e-13, rtol=1.0e-13)
+                    for i2 in 1:ngrid_x2
+                        for i1 in 1:ngrid_x1
+                            for j2 in 1:ngrid_x2
+                                for j1 in 1:ngrid_x1
+                                    M_2D_exact[j1,j2,i1,i2] = M_x1[j1,i1]*M_x2[j2,i2]
+                                end
+                            end
+                        end
+                    end
+                    @. M_2D_err = abs(M_2D-M_2D_exact)
+                    @test maximum(M_2D_err) < atol
+                end
+            end
+        end
+    end
     return nothing
 end
 

From 5f7fe1daf8b3692b0589989b3659efd75c8d4a6a Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 28 Jan 2026 11:27:47 +0000
Subject: [PATCH 09/12] Improve comment.

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index b7d36e4..9299d0a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -436,7 +436,7 @@ function test_2D_linear_operators(;nodes::node_type=GLL,
             end
         end
     end
-    # test produces of 1D matrices for a simple separable Jacobian
+    # tests 2D matrices against products of 1D matrices for a simple separable Jacobian
     for l1A in (lagrange_x, d_lagrange_dx)
         for l1B in (lagrange_x, d_lagrange_dx)
             for l2A in (lagrange_x, d_lagrange_dx)

From 229adb35f94c051f4b65ffd336a58919e19d319f Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Thu, 29 Jan 2026 11:52:26 +0000
Subject: [PATCH 10/12] Print message if p-adaptive integration converges.

---
 src/FiniteElementMatrices.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index 304fad3..500fb1e 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -72,7 +72,7 @@ end
 # points is increased is the same for each coordinate.
 function finite_element_matrix(args...;
     atol=1.0e-12::Float64, rtol=1.0e-12::Float64, max_iterations=1::Int64,
-    quadrature_increment=5::Int64, kwargs...)
+    quadrature_increment=5::Int64, verbose=false::Bool, kwargs...)
     matrix = _finite_element_matrix(args...;
                 adaptive_quadrature_points=0,
                 kwargs...)
@@ -90,6 +90,10 @@ function finite_element_matrix(args...;
         # update the matrix
         matrix .= matrixnew
     end
+    if verbose && !(maxerr > atol + rtol*maxvalue)
+        println(
+            "finite_element_matrix() converged after $iteration iterations \n with maxerr = $maxerr, adaptive_quadrature_points=$(quadrature_increment*iteration) \n")
+    end
     if maxerr > atol + rtol*maxvalue && max_iterations > 0
         error("Unable to achieve specified error bounds: \n
               maxerr < atol + rtol*maxvalue \n

From 0456bbdfc2ee98fbbd6feccff8ad7736eccb8467 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Thu, 29 Jan 2026 12:01:37 +0000
Subject: [PATCH 11/12] Tidy up while loop and message conditions so that
 `tolerance_exceeded` is defined only on one line.

---
 src/FiniteElementMatrices.jl | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index 500fb1e..4059c02 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -79,7 +79,8 @@ function finite_element_matrix(args...;
     iteration = 0
     maxvalue = maximum(abs.(matrix))
     maxerr = maxvalue
-    while iteration < max_iterations && maxerr > atol + rtol*maxvalue
+    tolerance_exceeded = max_iterations > 0
+    while iteration < max_iterations && tolerance_exceeded
         # increment the iteration
         iteration += 1
         matrixnew = _finite_element_matrix(args...;
@@ -87,14 +88,18 @@ function finite_element_matrix(args...;
                         kwargs...)
         maxvalue = maximum(abs.(matrixnew))
         maxerr = maximum(abs.(matrixnew .- matrix))
+        tolerance_exceeded = maxerr > atol + rtol*maxvalue
         # update the matrix
         matrix .= matrixnew
     end
-    if verbose && !(maxerr > atol + rtol*maxvalue)
-        println(
-            "finite_element_matrix() converged after $iteration iterations \n with maxerr = $maxerr, adaptive_quadrature_points=$(quadrature_increment*iteration) \n")
+    if verbose && !(tolerance_exceeded)
+        if iteration > 0
+            println("finite_element_matrix() converged after $iteration iterations \n with maxerr = $maxerr, adaptive_quadrature_points=$(quadrature_increment*iteration) \n")
+        else
+            println("finite_element_matrix() returned a matrix without p-adaptive quadrature refinement.")
+        end
     end
-    if maxerr > atol + rtol*maxvalue && max_iterations > 0
+    if tolerance_exceeded
         error("Unable to achieve specified error bounds: \n
               maxerr < atol + rtol*maxvalue \n
               after max_iterations = $max_iterations with atol=$atol and rtol=$rtol we obtain \n

From c51e83b325cdca22d757ad31df8b428d26be3b77 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Tue, 24 Feb 2026 12:15:49 +0000
Subject: [PATCH 12/12] Update README.md.

---
 README.md | 49 +++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 43 insertions(+), 6 deletions(-)

diff --git a/README.md b/README.md
index c142b8e..c8d3088 100644
--- a/README.md
+++ b/README.md
@@ -3,7 +3,7 @@
 A package for computing unassembled matrices appearing in $`C^0`$ continuous Galerkin
 finite element models.
 
-# Methods
+# Matrices in one coordinate (1D)
 
 This package computes matrices of the following form.
 ```math
@@ -45,6 +45,10 @@ N_{ijk} = \int^{x_u}_{x_l} P_i(x) Q_j(x) R_k(x) \rho(x) d x = \int^{1}_{-1} S_i(
 where $`R_i(x)`$ may be either of $`\Phi(x)`$ or $`d\Phi/dx`$ and
 $`U_i(z)`$ may be either of $`l_i(z)`$ or $`(1/s)d l_i(z)/d z`$.
 
+# p-adaptive Gaussian quadrature
+
+In order to guarantee that the matrices $`A_{ij}`$ with non-polynomial kernels $`\rho(x)`$ are calculated with sufficient accuracy, we use p-adaptive quadrature with Gauss-Legendre quadrature rules from `FastGaussQuadrature` at each iteration. We iteratively calculate the matrices with an increasingly large number of quadrature points and iterate until the result is constant up to absolute and relative tolerances. The inputs used to control this process are described in the section on usage below. For polynomial kernels we can deduce the number of quadrature points for exact results, avoiding the need for the iterative calculation in the pure polynomial case.
+
 # Usage
 
 To create these matrices, the user must supply the key grid
@@ -75,6 +79,16 @@ using FiniteElementMatrices
 coordinate = ElementCoordinates(x,scale,shift)
 M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
 ```
+Note that the order in which the `@enum` Lagrange polynomials are given corresponds to the indices of `M`. E.g., for a matrix describing a first derivative
+```
+l_i = lagrange_x
+l_j = d_lagrange_dx
+P = finite_element_matrix(l_i,l_j,0,coordinate)
+```
+the number `P[i,j]` corresponds to the matrix element
+```math
+P_{ij} = \int^{x_u}_{x_l} \Phi_i(x) \frac{d \Phi_j}{d x} d x.
+```
 For a stiffness matrix describing a second derivative, we could
 use the following.
 ```
@@ -108,13 +122,36 @@ To specify an arbitrary kernel function $`\rho(x)`$, the following function call
 ```
 coordinate = ElementCoordinates(x,scale,shift)
 M = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
-            kernel_function=(v -> rho(v)),additional_quadrature_points=n)
+            kernel_function=(v -> rho(v)),additional_quadrature_points=n,
+            quadrature_increment=m,
+            max_iterations=q, atol=atol, rtol=rtol)
+```
+where we specify some function `rho(v)` to be the kernel, and we specify that we `n` additional quadrature
+points should be used in addition to the number assumed from the size of the reference grid `x`. If `q` is greater than zero, for each iteration we proceed to calculate the matrix again with a further additional `m` quadrature points, until the matrix is calculated to the specified absolute tolerance `atol` and relative tolerance `rtol`, or we call `error()` with an appropriate message for the user. Note that `rho(v)` should be specified in terms of the physical coordinate including the scale and shift factors `v = scale * x + shift`, not in terms of the reference coordinate `x`.
+
+# Matrices in two coordinates (2D)
+
+It is possible to write down weak forms where the multi-dimensional integrals do not separate into a product of one-dimensional integrals, even for a tensor-product basis. This occurs where the kernel is not a separable function, e.g., for the matrix
+```math
+K_{ikjn} = \int^{x_{2u}}_{x_{2l}}\int^{x_{1u}}_{x_{1l}} P_i(x_1) R_k(x_2) Q_j(x_1) S_n(x_2)\rho(x_1,x_2) d x_1 d x_2,
 ```
-where we specify some function `rho(v)` to be the kernel, and we specify than `n` additional quadrature
-points should be used in addition to the number assumed from the size of the reference grid `x`.
-Note that `rho(v)` should be specified in terms of the physical coordinate including the scale and shift
-factors `v = scale * x + shift`, not in terms of the reference coordinate `x`.
+where $`P_i(x_1)`$, $`R_i(x_2)`$, $`Q_i(x_1)`$, $`S_i(x_2)`$ are from the set of Lagrange polynomials basis
+functions (and their first derivatives) in $`x_1`$ and $`x_2`$, respectively.
 
+The syntax that we use to create matrices in two coordinates is as follows
+```
+l_i = lagrange_x
+l_j = lagrange_x
+l_k = d_lagrange_dx
+l_n = d_lagrange_dx
+K_2D = finite_element_matrix(l_i,l_j,coordinate_x1,
+                            l_k,l_n,coordinate_x2;
+                            kernel_function=((x1,x2)->rho(x1,x2)),
+                            max_iterations=10, quadrature_increment=10,
+                            atol=1.0e-13, rtol=1.0e-13)
+```
+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.
 # Examples
 
 There are several examples of taking first and second derivatives