Skip to content

Integrand function inputs #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mrhardman merged 5 commits into from
Jan 7, 2026
Merged

Integrand function inputs #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
6a272dc
Add the option to pass coordinate functions into `finite_element_matr…
mrhardman Jan 5, 2026
b83cc73
Test nonpolynomial kernels in `finite_element_matrix()` for both the …
mrhardman Jan 5, 2026
e48af0f
Rename `function_of_v` to `kernel_function`.
mrhardman Jan 5, 2026
df958e2
Update README.md
mrhardman Jan 5, 2026
ceb41c5
Correct README.md.
mrhardman Jan 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 20 additions & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -33,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) (s z + c)^p 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) (s z + c)^p 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

Expand Down Expand Up @@ -99,6 +104,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
Expand Down
53 changes: 34 additions & 19 deletions src/FiniteElementMatrices.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
kernel_function=((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]
kernel_function=((v -> 1.0))::Function,
additional_quadrature_points=0::Int64,
)
lpoly_data = coordinate.lpoly_data
ngrid = length(coordinate.lpoly_data.x_nodes)
scale = coordinate.scale
Expand All @@ -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
Expand All @@ -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 = 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))
end
Expand All @@ -130,6 +131,20 @@ function finite_element_matrix(
power::Int64,
coordinate::ElementCoordinates
)
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(
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]
kernel_function=((v -> 1.0))::Function,
additional_quadrature_points=0::Int64,
)
lpoly_data = coordinate.lpoly_data
ngrid = length(coordinate.lpoly_data.x_nodes)
scale = coordinate.scale
Expand All @@ -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
Expand All @@ -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 = 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)*
prefactor3*lagrange3(kth_lpoly_data,zzl))
Expand Down
98 changes: 98 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
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)
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,
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)
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")
Expand Down Expand Up @@ -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

Expand Down