Skip to content

Type specialization #1

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 6 commits into from
Mar 4, 2026
Merged

Type specialization #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
339b902
Attempt to specialize on types in structs.
mrhardman Feb 27, 2026
890e53c
Use where statement for function arguments.
mrhardman Mar 4, 2026
77a2af7
Use where with AbstractBoundaryConditionType.
mrhardman Mar 4, 2026
865ce48
Add type specialization to function interfaces for boundary conditions.
mrhardman Mar 4, 2026
72a38b9
Rename `AbstractBoundaryConditionType` to `AbstractBoundaryCondition`.
mrhardman Mar 4, 2026
8e427ef
Update README.md
mrhardman Mar 4, 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
26 changes: 15 additions & 11 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,14 +69,18 @@ Once we have created a set of 1D coordinates that represent a N-D physical syste
```
# method for assembling operators in one coordinate
assemble_operator(weak_form::Array{Float64,3},
x::FiniteElementCoordinate,
boundary_condition_x::BoundaryConditionType)
x::FiniteElementCoordinate,
boundary_condition_x::Tbc
) where Tbc <: AbstractBoundaryCondition
# method for assembling operators in two coordinates
assemble_operator(weak_form::Function,
x::FiniteElementCoordinate,
y::FiniteElementCoordinate,
boundary_condition_x::BoundaryConditionType,
boundary_condition_y::BoundaryConditionType)
assemble_operator(weak_form::TF,
x::FiniteElementCoordinate,
y::FiniteElementCoordinate,
boundary_condition_x::Tbcx,
boundary_condition_y::Tbcy
) where {TF <: Function,
Tbcx <: AbstractBoundaryCondition,
Tbcy <: AbstractBoundaryCondition}
```
The returned value is of type `AbstractSparseArray{Float64,Int64,2}`.

Expand Down Expand Up @@ -114,10 +118,10 @@ and use `lu_Laplacian_Dirichlet_bc` with `ldiv!` as usual. See `test/PoissonSolv

The internally supported boundary conditions are
```
abstract type BoundaryConditionType end
struct NaturalBC <: BoundaryConditionType end
struct DirichletBC <: BoundaryConditionType end
struct PeriodicBC <: BoundaryConditionType end
abstract type AbstractBoundaryCondition end
struct NaturalBC <: AbstractBoundaryCondition end
struct DirichletBC <: AbstractBoundaryCondition end
struct PeriodicBC <: AbstractBoundaryCondition end
```
The user can define their own boundary condition types and extend the boundary condition function `impose_boundary_condition()` as required. Note that we explicitly supply the boundary condition choice for each call to `assemble_operator()`. The option to pass a boundary condition instance to the coordinate struct is given, e.g., from `test/PeriodicBcTests.jl`,
```
Expand Down
10 changes: 6 additions & 4 deletions src/Calculus.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,10 @@ points in vx. As a result, this means that this integral
is only exact for polynomials with order up to 2N-3 where
N=`ngrid` is the number of collocation points per element.
"""
function integral(prefactor::Function,
function integral(prefactor::TF,
integrand::AbstractArray{Float64,1},
vx::FiniteElementCoordinate)
vx::FiniteElementCoordinate
) where TF <: Function
@boundscheck (vx.n,) == size(integrand) || throw(BoundsError(integrand))
integral = 0.0
for ivx ∈ eachindex(vx.grid)
Expand All @@ -29,10 +30,11 @@ points in vx and vy. As a result, this means that this integral
is only exact for polynomials with order up to 2N-3 where
N=`ngrid` is the number of collocation points per element.
"""
function integral(prefactor::Function,
function integral(prefactor::TF,
integrand::AbstractArray{Float64,2},
vx::FiniteElementCoordinate,
vy::FiniteElementCoordinate)
vy::FiniteElementCoordinate
) where TF <: Function
@boundscheck (vx.n, vy.n) == size(integrand) || throw(BoundsError(integrand))
integral = 0.0
for ivy ∈ eachindex(vy.grid), ivx ∈ eachindex(vx.grid)
Expand Down
105 changes: 61 additions & 44 deletions src/FiniteElementAssembly.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ export FiniteElementCoordinate,
exclude_lower_boundary_point, exclude_upper_boundary_point,
BoundaryPointType,
impose_boundary_condition,
BoundaryConditionType, NaturalBC, DirichletBC, PeriodicBC,
AbstractBoundaryCondition, NaturalBC, DirichletBC, PeriodicBC,
assemble_operator,
# testing
set_element_boundaries,
Expand All @@ -23,15 +23,15 @@ using FiniteElementMatrices: ElementCoordinates,
finite_element_matrix
using FastGaussQuadrature: gausslegendre, gausslobatto, gaussradau

abstract type BoundaryConditionType end
struct NaturalBC <: BoundaryConditionType end
struct DirichletBC <: BoundaryConditionType end
struct PeriodicBC <: BoundaryConditionType end
abstract type AbstractBoundaryCondition end
struct NaturalBC <: AbstractBoundaryCondition end
struct DirichletBC <: AbstractBoundaryCondition end
struct PeriodicBC <: AbstractBoundaryCondition end

"""
struct containing information for first derivatives
"""
struct FirstDerivativeData
struct FirstDerivativeData{TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# Local mass matrix = \int phi_i(x) phi_j(x) dx
MM::Array{Float64,3}
# Local first derivative matrix \int phi_i(x) phi'_j(x) dx
Expand All @@ -41,9 +41,9 @@ struct FirstDerivativeData
# dummy array for storing solution
dummy_df::Array{Float64,1}
# Assembled 1D mass matrix
MM1D::AbstractSparseArray{Float64,Int64,2}
MM1D::TMatrix
# Assembled 1D first derivative matrix
PP1D::AbstractSparseArray{Float64,Int64,2}
PP1D::TMatrix
# LU object for mass matrix solve
lu_MM1D::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}
"""
Expand All @@ -55,7 +55,7 @@ struct FirstDerivativeData
n::Int64,
igrid_full::Array{Int64,2},
element_data::Array{ElementCoordinates,1},
boundary_condition::BoundaryConditionType)
boundary_condition::Tbc) where Tbc <: AbstractBoundaryCondition
# Local mass matrix MM[j,i] = \int phi_i(x) phi_j(x) dx
MM = Array{Float64,3}(undef,ngrid,ngrid,nelement)
# PP[j,i] = \int phi_i(x) phi'_j(x) dx
Expand All @@ -78,7 +78,7 @@ struct FirstDerivativeData
MM1D_sparse = sparse(MM1D)
PP1D_sparse = sparse(PP1D)
lu_MM1D = lu(MM1D_sparse)
return new(MM,PP,dummy_rhs,dummy_df,
return new{Float64,typeof(MM1D)}(MM,PP,dummy_rhs,dummy_df,
MM1D_sparse,PP1D_sparse,lu_MM1D)
end
end
Expand Down Expand Up @@ -107,7 +107,7 @@ function assemble_1D_operator(QQ1D_local::Array{Float64,3},
nelement::Int64,
n::Int64,
igrid_full::Array{Int64,2},
boundary_condition::BoundaryConditionType)
boundary_condition::Tbc) where Tbc <: AbstractBoundaryCondition
# create the 1D constructor arrays
nsparse = (nelement - 1)*(ngrid^2 - 1) + ngrid^2
II = zeros(Int64,nsparse)
Expand Down Expand Up @@ -218,7 +218,7 @@ end
"""
structure containing basic information related to coordinates
"""
struct FiniteElementCoordinate
struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition}
# name is the name of the variable associated with this coordiante
name::String
# n is the total number of grid points associated with this coordinate
Expand Down Expand Up @@ -246,7 +246,7 @@ struct FiniteElementCoordinate
# Note that the boundary condition must be separately supplied in the construction
# of other operators, since not all operators have the same boundary conditions on
# the same coordinates.
bc::BoundaryConditionType
bc::Tbc
# wgts contains the integration weights associated with each grid point
wgts::Array{Float64,1}
# scale for each element
Expand Down Expand Up @@ -277,9 +277,10 @@ struct FiniteElementCoordinate
scalar_input::ScalarCoordinateInputs;
# the 1D Jacobian or weight function rho(x) that
# appears in 1D integrals as \int (.) rho(x) d x
weight_function=((x)-> 1.0 )::Function,
weight_function::TF=((x)-> 1.0 ),
# which boundary condition to store and use for first derivatives
bc=NaturalBC()::BoundaryConditionType)
bc::Tbc=NaturalBC()
) where {TF <: Function, Tbc <: AbstractBoundaryCondition}
ngrid = scalar_input.ngrid
nelement = scalar_input.nelement
# initialise the data used to construct the grid
Expand Down Expand Up @@ -324,9 +325,10 @@ struct FiniteElementCoordinate
element_data::Union{Array{ElementCoordinates,1},Nothing};
# the 1D Jacobian or weight function rho(x) that
# appears in 1D integrals as \int (.) rho(x) d x
weight_function=((x)-> 1.0 )::Function,
weight_function::TF=((x)-> 1.0 ),
# which boundary condition to store and use in first derivatives
bc=NaturalBC()::BoundaryConditionType)
bc::Tbc=NaturalBC()
) where {TF <: Function, Tbc <: AbstractBoundaryCondition}
if typeof(element_data) == Nothing
# this is a trivial coordinate of length 1
nelement = 1
Expand Down Expand Up @@ -430,7 +432,7 @@ struct FiniteElementCoordinate
if typeof(bc) == PeriodicBC && (!include_lower_boundary || !include_upper_boundary)
error("Must include upper and lower boundaries on the grid for typeof(bc) = $(typeof(bc))")
end
return new(name, n_global, ngrid,
return new{typeof(bc)}(name, n_global, ngrid,
nelement, domainLength, grid, igrid, ielement, imin, imax,
igrid_full, bc, wgts,
element_scale, element_shift, element_boundaries,
Expand Down Expand Up @@ -607,14 +609,16 @@ end
Methods for imposing boundary conditions on 1D operators
"""
function impose_boundary_condition_x(boundary_condition::Union{NaturalBC,DirichletBC},
operator_sparse::AbstractSparseArray{Float64,Int64,2},
n::Int64, ngrid::Int64, QQ::AbstractArray{Float64,3})
operator_sparse::TMatrix,
n::Int64, ngrid::Int64, QQ::Array{Float64,3}
) where {TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# do nothing
return nothing
end
function impose_boundary_condition_x(boundary_condition::PeriodicBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
n::Int64, ngrid::Int64, QQ::AbstractArray{Float64,3})
operator_sparse::TMatrix,
n::Int64, ngrid::Int64, QQ::Array{Float64,3}
) where {TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# periodic BC contributions
# set lower row to (1, 0, ..., 0, -1)
operator_sparse[1,:] .= 0.0
Expand All @@ -630,16 +634,18 @@ end
Methods for imposing boundary conditions on 2D operators
"""
function impose_boundary_condition_x(boundary_condition::NaturalBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# do nothing
return nothing
end
function impose_boundary_condition_x(boundary_condition::DirichletBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
if x.include_lower_boundary
# set row to (1, 0, 0, ..., 0)
for iy in 1:y.n
Expand All @@ -659,9 +665,10 @@ function impose_boundary_condition_x(boundary_condition::DirichletBC,
return nothing
end
function impose_boundary_condition_x(boundary_condition::PeriodicBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
if x.include_lower_boundary
# set row to (1, 0, 0, ..., 0, -1)
for iy in 1:y.n
Expand Down Expand Up @@ -707,16 +714,18 @@ function impose_boundary_condition_x(boundary_condition::PeriodicBC,
end

function impose_boundary_condition_y(boundary_condition::NaturalBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# do nothing
return nothing
end
function impose_boundary_condition_y(boundary_condition::DirichletBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
if y.include_lower_boundary
# set row to (1, 0, 0, ..., 0)
for ix in 1:x.n
Expand All @@ -736,9 +745,10 @@ function impose_boundary_condition_y(boundary_condition::DirichletBC,
return nothing
end
function impose_boundary_condition_y(boundary_condition::PeriodicBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
if y.include_lower_boundary
# set row to (1, 0, 0, ..., 0, -1)
for ix in 1:x.n
Expand Down Expand Up @@ -785,9 +795,10 @@ end

function impose_boundary_condition(
boundary_condition_x::PeriodicBC,boundary_condition_y::Union{NaturalBC,DirichletBC},
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# impose periodic BC first to ensure Dirichlet BC is imposed last
impose_boundary_condition_x(boundary_condition_x,operator_sparse,x,y,weak_form)
impose_boundary_condition_y(boundary_condition_y,operator_sparse,x,y,weak_form)
Expand All @@ -796,19 +807,21 @@ function impose_boundary_condition(
end
function impose_boundary_condition(
boundary_condition_x::Union{NaturalBC,DirichletBC},boundary_condition_y::PeriodicBC,
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# impose periodic BC first to ensure Dirichlet BC is imposed last
impose_boundary_condition_y(boundary_condition_y,operator_sparse,x,y,weak_form)
impose_boundary_condition_x(boundary_condition_x,operator_sparse,x,y,weak_form)
return nothing
end
function impose_boundary_condition(
boundary_condition_x::Union{NaturalBC,DirichletBC},boundary_condition_y::Union{NaturalBC,DirichletBC},
operator_sparse::AbstractSparseArray{Float64,Int64,2},
operator_sparse::TMatrix,
x::FiniteElementCoordinate,y::FiniteElementCoordinate,
weak_form::Function)
weak_form::TF) where {TF <: Function, TFloat <: Real,
TMatrix <: AbstractSparseArray{TFloat,Int64,2}}
# order of imposition is unimportant
impose_boundary_condition_x(boundary_condition_x,operator_sparse,x,y,weak_form)
impose_boundary_condition_y(boundary_condition_y,operator_sparse,x,y,weak_form)
Expand All @@ -819,7 +832,8 @@ External method for assembling operators in one coordinate
"""
function assemble_operator(weak_form::Array{Float64,3},
x::FiniteElementCoordinate,
boundary_condition_x::BoundaryConditionType)
boundary_condition_x::Tbc
) where Tbc <: AbstractBoundaryCondition
return assemble_1D_operator(weak_form,
x.ngrid, x.nelement, x.n, x.igrid_full,
boundary_condition_x)
Expand All @@ -828,11 +842,14 @@ end
"""
Method for assembling operators in two coordinates
"""
function assemble_operator(weak_form::Function,
function assemble_operator(weak_form::TF,
x::FiniteElementCoordinate,
y::FiniteElementCoordinate,
boundary_condition_x::BoundaryConditionType,
boundary_condition_y::BoundaryConditionType)
boundary_condition_x::Tbcx,
boundary_condition_y::Tbcy
) where {TF <: Function,
Tbcx <: AbstractBoundaryCondition,
Tbcy <: AbstractBoundaryCondition}
# Assemble a 2D mass matrix in the global compound coordinate
# total number of non-zero element is the maximum index of
# the sparse matrix index, by construction
Expand Down