From 6a272dcf9432749f370c27dafda731d9cb82e385 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 5 Jan 2026 09:16:44 +0000
Subject: [PATCH 1/5] Add the option to pass coordinate functions into
 `finite_element_matrix()`. The functions passed in must be expressed in the
 physical coordinate `v = s z + c` rather than the reference coordinate `z` on
 the element which has a range [-1,1].

---
 src/FiniteElementMatrices.jl | 53 +++++++++++++++++++++++-------------
 1 file changed, 34 insertions(+), 19 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index e52fb5e..b0b100c 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -68,25 +68,26 @@ function select_lagrange_prefactor(fn_type::LagrangeFunctionType,
     return prefactor
 end
 
-function get_polyz(power::Int64,
-                    scale::Float64,
-                    shift::Float64,
-                    zz::Float64)
-    polyz = 0.0
-    for q in 0:power
-        polyz += (binomial(power,q)*
-                    ((scale*zz)^q)*
-                    (shift^(power-q)))
-    end
-    return polyz
-end
-
 function finite_element_matrix(
     fn1_type::LagrangeFunctionType,
     fn2_type::LagrangeFunctionType,
     power::Int64,
     coordinate::ElementCoordinates
     )
+    return finite_element_matrix(fn1_type, fn2_type, coordinate;
+                function_of_v=((v -> v^power)),
+                additional_quadrature_points=power)
+end
+
+function finite_element_matrix(
+    fn1_type::LagrangeFunctionType,
+    fn2_type::LagrangeFunctionType,
+    coordinate::ElementCoordinates;
+    # function of the "physical" coord v = s z + c
+    # rather than of the reference coordinate z on [-1,1]
+    function_of_v=((v -> 1.0))::Function,
+    additional_quadrature_points=0::Int64,
+    )
     lpoly_data = coordinate.lpoly_data
     ngrid = length(coordinate.lpoly_data.x_nodes)
     scale = coordinate.scale
@@ -100,7 +101,7 @@ function finite_element_matrix(
     prefactor1 = select_lagrange_prefactor(fn1_type,scale)
     prefactor2 = select_lagrange_prefactor(fn2_type,scale)
     # nquad chosen for exact results for all power, ngrid
-    nquad = ngrid + power
+    nquad = ngrid + additional_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
     # int P_i(z) Q_j(z) poly(z) s d z
@@ -113,8 +114,8 @@ function finite_element_matrix(
             ith_lpoly_data = lpoly_data.lpoly_data[i]
             for l in 1:nquad
                 zzl = zz[l]
-                polyz = get_polyz(power,scale,shift,zzl)
-                matrix[i,j] += (scale*wz[l]*polyz*
+                funcz = function_of_v(scale*zzl+shift)::Float64
+                matrix[i,j] += (scale*wz[l]*funcz*
                            prefactor1*lagrange1(ith_lpoly_data,zzl)*
                            prefactor2*lagrange2(jth_lpoly_data,zzl))
             end
@@ -130,6 +131,20 @@ function finite_element_matrix(
     power::Int64,
     coordinate::ElementCoordinates
     )
+    return finite_element_matrix(fn1_type, fn2_type, fn3_type,
+                coordinate; function_of_v=((v -> v^power)),
+                additional_quadrature_points=power)
+end
+function finite_element_matrix(
+    fn1_type::LagrangeFunctionType,
+    fn2_type::LagrangeFunctionType,
+    fn3_type::LagrangeFunctionType,
+    coordinate::ElementCoordinates;
+    # function of the "physical" coord v = s z + c
+    # rather than of the reference coordinate z on [-1,1]
+    function_of_v=((v -> 1.0))::Function,
+    additional_quadrature_points=0::Int64,
+    )
     lpoly_data = coordinate.lpoly_data
     ngrid = length(coordinate.lpoly_data.x_nodes)
     scale = coordinate.scale
@@ -145,7 +160,7 @@ function finite_element_matrix(
     prefactor2 = select_lagrange_prefactor(fn2_type,scale)
     prefactor3 = select_lagrange_prefactor(fn3_type,scale)
     # nquad chosen for exact results for all power, ngrid
-    nquad = 2*ngrid + power
+    nquad = 2*ngrid + additional_quadrature_points
     zz, wz = gausslegendre(nquad)
     # compute integral
     # int P_i(z) Q_j(z) S_k(z) poly(z) d z
@@ -160,8 +175,8 @@ function finite_element_matrix(
                 ith_lpoly_data = lpoly_data.lpoly_data[i]
                 for l in 1:nquad
                     zzl = zz[l]
-                    polyz = get_polyz(power,scale,shift,zzl)
-                    matrix[i,j,k] += (scale*wz[l]*polyz*
+                    funcz = function_of_v(scale*zzl+shift)
+                    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))

From b83cc739656832732d59051e7b564f2df9d6ce8c Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 5 Jan 2026 10:31:16 +0000
Subject: [PATCH 2/5] Test nonpolynomial kernels in `finite_element_matrix()`
 for both the linear and nonlinear operators.

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

diff --git a/test/runtests.jl b/test/runtests.jl
index 894d219..ea7fe76 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -290,6 +290,95 @@ function test_nonlinear_operators(;nodes::node_type=GLL,
     @test maxerr < atol
 end
 
+function test_nonpolynomial_kernel(;nodes::node_type=GLL,
+                            ngrid::Int64=20,
+                            func::Function=(x -> sin(x)),
+                            jacobian::Function=(x -> exp(-x)),
+                            additional_quadrature_points::Int64 = 5,
+                            y0::Float64=-1.0,
+                            y1::Float64=1.0,
+                            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)
+    M_with_kernel = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
+            function_of_v=jacobian,additional_quadrature_points=5)
+    M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
+    # function to test
+    f = Array{Float64,1}(undef,ngrid)
+    h = Array{Float64,1}(undef,ngrid)
+    dummy = Array{Float64,1}(undef,ngrid)
+    result = Array{Float64,1}(undef,ngrid)
+    err = Array{Float64,1}(undef,ngrid)
+    for i in 1:ngrid
+        # change to physical coordinate y
+        # from reference coordinate x
+        y = scale*x[i] + shift
+        f[i] = func(y)
+        h[i] = jacobian(y)*func(y)
+    end
+    luM = lu(M)
+    mul!(dummy,M_with_kernel,f)
+    ldiv!(result,luM,dummy)
+    @. err = result - h
+    # println(f)
+    # println(h)
+    # println(result)
+    maxerr = maximum(abs.(err))
+    @test maxerr < atol
+    #println(maxerr)
+end
+
+function test_nonlinear_operator_nonpolynomial_kernel(;nodes::node_type=GLL,
+                            ngrid::Int64=20,
+                            func1::Function=(x -> sin(x)),
+                            func2::Function=(x -> x^2),
+                            jacobian::Function=( x -> exp(-x)),
+                            additional_quadrature_points::Int64 = 5,
+                            y0::Float64=-1.0,
+                            y1::Float64=1.0,
+                            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)
+    Y000_with_kernel = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,coordinate,
+            function_of_v=jacobian,additional_quadrature_points=5)
+    M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
+    # function to test
+    u = Array{Float64,1}(undef,ngrid)
+    v = Array{Float64,1}(undef,ngrid)
+    h = Array{Float64,1}(undef,ngrid)
+    dummy = Array{Float64,1}(undef,ngrid)
+    result = Array{Float64,1}(undef,ngrid)
+    err = Array{Float64,1}(undef,ngrid)
+    for i in 1:ngrid
+        # change to physical coordinate y
+        # from reference coordinate x
+        y = scale*x[i] + shift
+        u[i] = func1(y)
+        v[i] = func2(y)
+        h[i] = jacobian(y)*func1(y)*func2(y)
+    end
+    luM = lu(M)
+    @. dummy = 0.0
+    for k in 1:ngrid
+        for j in 1:ngrid
+            for i in 1:ngrid
+                dummy[k] += Y000_with_kernel[i,j,k]*u[i]*v[j]
+            end
+        end
+    end
+    ldiv!(result,luM,dummy)
+    @. err = result - h
+    # println(h)
+    # println(result)
+    maxerr = maximum(abs.(err))
+    @test maxerr < atol
+    #println(maxerr)
+end
+
 function runtests()
     @testset "FiniteElementMatrices" begin
         println("test FiniteElementMatrices")
@@ -343,6 +432,15 @@ function runtests()
                                     y0 = -0.6, y1=1.0,
                                     ngrid=n)
         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
     end
 end
 

From e48af0fcee912d6ee305af88992a2a6cf38a8635 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 5 Jan 2026 10:43:02 +0000
Subject: [PATCH 3/5] Rename `function_of_v` to `kernel_function`.

---
 src/FiniteElementMatrices.jl | 12 ++++++------
 test/runtests.jl             |  4 ++--
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/FiniteElementMatrices.jl b/src/FiniteElementMatrices.jl
index b0b100c..39234ac 100644
--- a/src/FiniteElementMatrices.jl
+++ b/src/FiniteElementMatrices.jl
@@ -75,7 +75,7 @@ function finite_element_matrix(
     coordinate::ElementCoordinates
     )
     return finite_element_matrix(fn1_type, fn2_type, coordinate;
-                function_of_v=((v -> v^power)),
+                kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
 end
 
@@ -85,7 +85,7 @@ function finite_element_matrix(
     coordinate::ElementCoordinates;
     # function of the "physical" coord v = s z + c
     # rather than of the reference coordinate z on [-1,1]
-    function_of_v=((v -> 1.0))::Function,
+    kernel_function=((v -> 1.0))::Function,
     additional_quadrature_points=0::Int64,
     )
     lpoly_data = coordinate.lpoly_data
@@ -114,7 +114,7 @@ function finite_element_matrix(
             ith_lpoly_data = lpoly_data.lpoly_data[i]
             for l in 1:nquad
                 zzl = zz[l]
-                funcz = function_of_v(scale*zzl+shift)::Float64
+                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))
@@ -132,7 +132,7 @@ function finite_element_matrix(
     coordinate::ElementCoordinates
     )
     return finite_element_matrix(fn1_type, fn2_type, fn3_type,
-                coordinate; function_of_v=((v -> v^power)),
+                coordinate; kernel_function=((v -> v^power)),
                 additional_quadrature_points=power)
 end
 function finite_element_matrix(
@@ -142,7 +142,7 @@ function finite_element_matrix(
     coordinate::ElementCoordinates;
     # function of the "physical" coord v = s z + c
     # rather than of the reference coordinate z on [-1,1]
-    function_of_v=((v -> 1.0))::Function,
+    kernel_function=((v -> 1.0))::Function,
     additional_quadrature_points=0::Int64,
     )
     lpoly_data = coordinate.lpoly_data
@@ -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 = function_of_v(scale*zzl+shift)
+                    funcz = kernel_function(scale*zzl+shift)
                     matrix[i,j,k] += (scale*wz[l]*funcz*
                             prefactor1*lagrange1(ith_lpoly_data,zzl)*
                             prefactor2*lagrange2(jth_lpoly_data,zzl)*
diff --git a/test/runtests.jl b/test/runtests.jl
index ea7fe76..1706ed3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -303,7 +303,7 @@ function test_nonpolynomial_kernel(;nodes::node_type=GLL,
     shift = 0.5*(y0+y1)
     coordinate = ElementCoordinates(x,scale,shift)
     M_with_kernel = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
-            function_of_v=jacobian,additional_quadrature_points=5)
+            kernel_function=jacobian,additional_quadrature_points=5)
     M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
     # function to test
     f = Array{Float64,1}(undef,ngrid)
@@ -344,7 +344,7 @@ function test_nonlinear_operator_nonpolynomial_kernel(;nodes::node_type=GLL,
     shift = 0.5*(y0+y1)
     coordinate = ElementCoordinates(x,scale,shift)
     Y000_with_kernel = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,coordinate,
-            function_of_v=jacobian,additional_quadrature_points=5)
+            kernel_function=jacobian,additional_quadrature_points=5)
     M = finite_element_matrix(lagrange_x,lagrange_x,0,coordinate)
     # function to test
     u = Array{Float64,1}(undef,ngrid)

From df958e2a48bc3e83cb8cc3813d5de68233c3f674 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Mon, 5 Jan 2026 10:46:44 +0000
Subject: [PATCH 4/5] Update README.md

---
 README.md | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index db3fb2b..12382c2 100644
--- a/README.md
+++ b/README.md
@@ -7,14 +7,16 @@ finite element models.
 
 This package computes matrices of the following form.
 ```math
-A_{ij} = \int^{x_u}_{x_l} P_i(x) Q_j(x) x^p d x,
+A_{ij} = \int^{x_u}_{x_l} P_i(x) Q_j(x) \rho(x) d x,
 ```
 where $`P_i`$ and $`Q_i`$ are from the set of Lagrange polynomials basis
 functions (and their first derivatives) used for
 $`C^0`$ continuous Galerkin
 finite element models, $`x_u`$ and $`x_l`$
 are the upper and lower limits of the element in the
-coordinate $`x`$, and $`p`$ is an integer.
+coordinate $`x`$, and $`\rho(x)`$ is a function. We provide an interface
+where the user can specify the function $`\rho(x)`$ as an argument, or
+the user can pass an integer $`p`$, in which case $`\rho(x) = x^p`$.
 
 The basis functions we use here are
 ```math
@@ -33,11 +35,11 @@ is required for the problem of interest.
 To compute the matrix $`A_{ij}`$, we use the reference coordinate
 $`z`$ to write
 ```math
-A_{ij} = \int^{1}_{-1} P_i(z) Q_j(z) (s z + c)^p d z.
+A_{ij} = \int^{1}_{-1} P_i(z) Q_j(z) \rho(s z + c) d z.
 ```
 We can also compute nonlinear matrices which have the form
 ```math
-N_{ijk} = \int^{x_u}_{x_l} P_i(x) Q_j(x) R_k(x) x^p d x = \int^{1}_{-1} P_i(z) Q_j(z) Q_k(z) (s z + c)^p d z.
+N_{ijk} = \int^{x_u}_{x_l} P_i(x) Q_j(x) R_k(x) x^p d x = \int^{1}_{-1} P_i(z) Q_j(z) Q_k(z) \rho(s z + c) d z.
 ```
 
 # Usage
@@ -99,6 +101,17 @@ Finally, we can form a nonlinear stiffness matrix by inserting an extra enum arg
 N = finite_element_matrix(lagrange_x,lagrange_x,lagrange_x,1,coordinate)
 ```
 
+To specify an arbitrary kernel function $`\rho(x)`$, the following function call should be used
+```
+coordinate = ElementCoordinates(x,scale,shift)
+M = finite_element_matrix(lagrange_x,lagrange_x,coordinate,
+            kernel_function=(v -> rho(v)),additional_quadrature_points=n)
+```
+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`.
+
 # Examples
 
 There are several examples of taking first and second derivatives

From ceb41c5c59ebbfbaf1b90f6a11bc989bb8c54f99 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Tue, 6 Jan 2026 16:30:07 +0000
Subject: [PATCH 5/5] Correct README.md.

---
 README.md | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 12382c2..6471951 100644
--- a/README.md
+++ b/README.md
@@ -35,12 +35,15 @@ is required for the problem of interest.
 To compute the matrix $`A_{ij}`$, we use the reference coordinate
 $`z`$ to write
 ```math
-A_{ij} = \int^{1}_{-1} P_i(z) Q_j(z) \rho(s z + c) d z.
+A_{ij} = \int^{1}_{-1} S_i(z) T_j(z) \rho(s z + c) d z,
 ```
+where $`S_i(z)`$, $`T_i(z)`$ may be either of $`l_i(z)`$ or $`(1/s)d l_i(z)/d z`$.
 We can also compute nonlinear matrices which have the form
 ```math
-N_{ijk} = \int^{x_u}_{x_l} P_i(x) Q_j(x) R_k(x) x^p d x = \int^{1}_{-1} P_i(z) Q_j(z) Q_k(z) \rho(s z + c) d z.
+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(z) T_j(z) U_k(z) \rho(s z + c) d z.
 ```
+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`$.
 
 # Usage