From 339b902071bb118db6eb325eee45ccfb87babd7f Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Fri, 27 Feb 2026 16:46:30 +0000
Subject: [PATCH 1/6] Attempt to specialize on types in structs.

---
 src/FiniteElementAssembly.jl | 38 ++++++++++++++++++------------------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl
index c48e91d..8f44473 100644
--- a/src/FiniteElementAssembly.jl
+++ b/src/FiniteElementAssembly.jl
@@ -7,7 +7,7 @@ export FiniteElementCoordinate,
     exclude_lower_boundary_point, exclude_upper_boundary_point,
     BoundaryPointType,
     impose_boundary_condition,
-    BoundaryConditionType, NaturalBC, DirichletBC, PeriodicBC,
+    AbstractBoundaryConditionType, NaturalBC, DirichletBC, PeriodicBC,
     assemble_operator,
     # testing
     set_element_boundaries,
@@ -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 AbstractBoundaryConditionType end
+struct NaturalBC <: AbstractBoundaryConditionType end
+struct DirichletBC <: AbstractBoundaryConditionType end
+struct PeriodicBC <: AbstractBoundaryConditionType 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
@@ -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}
     """
@@ -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 <: AbstractBoundaryConditionType
         # 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
@@ -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
@@ -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 <: AbstractBoundaryConditionType
     # create the 1D constructor arrays
     nsparse = (nelement - 1)*(ngrid^2 - 1) + ngrid^2
     II = zeros(Int64,nsparse)
@@ -218,7 +218,7 @@ end
 """
 structure containing basic information related to coordinates
 """
-struct FiniteElementCoordinate
+struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
     # 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
@@ -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
@@ -279,7 +279,7 @@ struct FiniteElementCoordinate
         # appears in 1D integrals as \int (.) rho(x) d x
         weight_function=((x)-> 1.0 )::Function,
         # which boundary condition to store and use for first derivatives
-        bc=NaturalBC()::BoundaryConditionType)
+        bc=NaturalBC()::AbstractBoundaryConditionType)
         ngrid = scalar_input.ngrid
         nelement = scalar_input.nelement
         # initialise the data used to construct the grid
@@ -326,7 +326,7 @@ struct FiniteElementCoordinate
         # appears in 1D integrals as \int (.) rho(x) d x
         weight_function=((x)-> 1.0 )::Function,
         # which boundary condition to store and use in first derivatives
-        bc=NaturalBC()::BoundaryConditionType)
+        bc=NaturalBC()::AbstractBoundaryConditionType)
         if typeof(element_data) == Nothing
             # this is a trivial coordinate of length 1
             nelement = 1
@@ -430,7 +430,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,
@@ -819,7 +819,7 @@ 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::AbstractBoundaryConditionType)
     return assemble_1D_operator(weak_form,
             x.ngrid, x.nelement, x.n, x.igrid_full,
             boundary_condition_x)
@@ -831,8 +831,8 @@ Method for assembling operators in two coordinates
 function assemble_operator(weak_form::Function,
                     x::FiniteElementCoordinate,
                     y::FiniteElementCoordinate,
-                    boundary_condition_x::BoundaryConditionType,
-                    boundary_condition_y::BoundaryConditionType)
+                    boundary_condition_x::AbstractBoundaryConditionType,
+                    boundary_condition_y::AbstractBoundaryConditionType)
     # 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

From 890e53cbb673ead1a0031f858717c86af67c70b2 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 4 Mar 2026 10:50:57 +0000
Subject: [PATCH 2/6] Use where statement for function arguments.

---
 src/Calculus.jl              | 10 ++++++----
 src/FiniteElementAssembly.jl | 33 ++++++++++++++++++---------------
 2 files changed, 24 insertions(+), 19 deletions(-)

diff --git a/src/Calculus.jl b/src/Calculus.jl
index 8704ffc..5254f68 100644
--- a/src/Calculus.jl
+++ b/src/Calculus.jl
@@ -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)
@@ -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)
diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl
index 8f44473..24ed6f4 100644
--- a/src/FiniteElementAssembly.jl
+++ b/src/FiniteElementAssembly.jl
@@ -277,9 +277,10 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
         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()::AbstractBoundaryConditionType)
+        bc::AbstractBoundaryConditionType=NaturalBC()
+        ) where TF <: Function
         ngrid = scalar_input.ngrid
         nelement = scalar_input.nelement
         # initialise the data used to construct the grid
@@ -324,9 +325,10 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
         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()::AbstractBoundaryConditionType)
+        bc::AbstractBoundaryConditionType=NaturalBC()
+        ) where TF <: Function
         if typeof(element_data) == Nothing
             # this is a trivial coordinate of length 1
             nelement = 1
@@ -632,14 +634,14 @@ Methods for imposing boundary conditions on 2D operators
 function impose_boundary_condition_x(boundary_condition::NaturalBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     # do nothing
     return nothing
 end
 function impose_boundary_condition_x(boundary_condition::DirichletBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     if x.include_lower_boundary
         # set row to (1, 0, 0, ..., 0)
         for iy in 1:y.n
@@ -661,7 +663,7 @@ end
 function impose_boundary_condition_x(boundary_condition::PeriodicBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     if x.include_lower_boundary
         # set row to (1, 0, 0, ..., 0, -1)
         for iy in 1:y.n
@@ -709,14 +711,14 @@ end
 function impose_boundary_condition_y(boundary_condition::NaturalBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     # do nothing
     return nothing
 end
 function impose_boundary_condition_y(boundary_condition::DirichletBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     if y.include_lower_boundary
         # set row to (1, 0, 0, ..., 0)
         for ix in 1:x.n
@@ -738,7 +740,7 @@ end
 function impose_boundary_condition_y(boundary_condition::PeriodicBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     if y.include_lower_boundary
         # set row to (1, 0, 0, ..., 0, -1)
         for ix in 1:x.n
@@ -787,7 +789,7 @@ function impose_boundary_condition(
             boundary_condition_x::PeriodicBC,boundary_condition_y::Union{NaturalBC,DirichletBC},
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     # 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)
@@ -798,7 +800,7 @@ function impose_boundary_condition(
             boundary_condition_x::Union{NaturalBC,DirichletBC},boundary_condition_y::PeriodicBC,
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     # 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)
@@ -808,7 +810,7 @@ function impose_boundary_condition(
             boundary_condition_x::Union{NaturalBC,DirichletBC},boundary_condition_y::Union{NaturalBC,DirichletBC},
             operator_sparse::AbstractSparseArray{Float64,Int64,2},
             x::FiniteElementCoordinate,y::FiniteElementCoordinate,
-            weak_form::Function)
+            weak_form::TF) where TF <: Function
     # 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)
@@ -828,11 +830,12 @@ 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::AbstractBoundaryConditionType,
-                    boundary_condition_y::AbstractBoundaryConditionType)
+                    boundary_condition_y::AbstractBoundaryConditionType
+                    ) where TF <: Function
     # 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

From 77a2af705615ad0f5d6b405f147cf377de530633 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 4 Mar 2026 11:02:00 +0000
Subject: [PATCH 3/6] Use where with AbstractBoundaryConditionType.

---
 src/FiniteElementAssembly.jl | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl
index 24ed6f4..2966478 100644
--- a/src/FiniteElementAssembly.jl
+++ b/src/FiniteElementAssembly.jl
@@ -279,8 +279,8 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
         # appears in 1D integrals as \int (.) rho(x) d x
         weight_function::TF=((x)-> 1.0 ),
         # which boundary condition to store and use for first derivatives
-        bc::AbstractBoundaryConditionType=NaturalBC()
-        ) where TF <: Function
+        bc::Tbc=NaturalBC()
+        ) where {TF <: Function, Tbc <: AbstractBoundaryConditionType}
         ngrid = scalar_input.ngrid
         nelement = scalar_input.nelement
         # initialise the data used to construct the grid
@@ -327,8 +327,8 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
         # appears in 1D integrals as \int (.) rho(x) d x
         weight_function::TF=((x)-> 1.0 ),
         # which boundary condition to store and use in first derivatives
-        bc::AbstractBoundaryConditionType=NaturalBC()
-        ) where TF <: Function
+        bc::Tbc=NaturalBC()
+        ) where {TF <: Function, Tbc <: AbstractBoundaryConditionType}
         if typeof(element_data) == Nothing
             # this is a trivial coordinate of length 1
             nelement = 1
@@ -821,7 +821,8 @@ External method for assembling operators in one coordinate
 """
 function assemble_operator(weak_form::Array{Float64,3},
                             x::FiniteElementCoordinate,
-                            boundary_condition_x::AbstractBoundaryConditionType)
+                            boundary_condition_x::Tbc
+    ) where Tbc <: AbstractBoundaryConditionType
     return assemble_1D_operator(weak_form,
             x.ngrid, x.nelement, x.n, x.igrid_full,
             boundary_condition_x)
@@ -833,9 +834,11 @@ Method for assembling operators in two coordinates
 function assemble_operator(weak_form::TF,
                     x::FiniteElementCoordinate,
                     y::FiniteElementCoordinate,
-                    boundary_condition_x::AbstractBoundaryConditionType,
-                    boundary_condition_y::AbstractBoundaryConditionType
-                    ) where TF <: Function
+                    boundary_condition_x::Tbcx,
+                    boundary_condition_y::Tbcy
+                    ) where {TF <: Function,
+                        Tbcx <: AbstractBoundaryConditionType,
+                        Tbcy <: AbstractBoundaryConditionType}
     # 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

From 865ce4818b7cb017f143bb4fb9869e10e9a78ae0 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 4 Mar 2026 11:16:32 +0000
Subject: [PATCH 4/6] Add type specialization to function interfaces for
 boundary conditions.

---
 src/FiniteElementAssembly.jl | 55 +++++++++++++++++++++---------------
 1 file changed, 33 insertions(+), 22 deletions(-)

diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl
index 2966478..4d48255 100644
--- a/src/FiniteElementAssembly.jl
+++ b/src/FiniteElementAssembly.jl
@@ -609,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
@@ -632,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::TF) where TF <: 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::TF) where TF <: 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
@@ -661,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::TF) where TF <: 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
@@ -709,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::TF) where TF <: 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::TF) where TF <: 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
@@ -738,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::TF) where TF <: 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
@@ -787,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::TF) where TF <: 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)
@@ -798,9 +807,10 @@ 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::TF) where TF <: 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)
@@ -808,9 +818,10 @@ function impose_boundary_condition(
 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::TF) where TF <: 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)

From 72a38b937bb431e9854002baf87d52fd2e080acf Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 4 Mar 2026 11:41:39 +0000
Subject: [PATCH 5/6] Rename `AbstractBoundaryConditionType` to
 `AbstractBoundaryCondition`.

---
 src/FiniteElementAssembly.jl | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl
index 4d48255..0f215ae 100644
--- a/src/FiniteElementAssembly.jl
+++ b/src/FiniteElementAssembly.jl
@@ -7,7 +7,7 @@ export FiniteElementCoordinate,
     exclude_lower_boundary_point, exclude_upper_boundary_point,
     BoundaryPointType,
     impose_boundary_condition,
-    AbstractBoundaryConditionType, NaturalBC, DirichletBC, PeriodicBC,
+    AbstractBoundaryCondition, NaturalBC, DirichletBC, PeriodicBC,
     assemble_operator,
     # testing
     set_element_boundaries,
@@ -23,10 +23,10 @@ using FiniteElementMatrices: ElementCoordinates,
                              finite_element_matrix
 using FastGaussQuadrature: gausslegendre, gausslobatto, gaussradau
 
-abstract type AbstractBoundaryConditionType end
-struct NaturalBC <: AbstractBoundaryConditionType end
-struct DirichletBC <: AbstractBoundaryConditionType end
-struct PeriodicBC <: AbstractBoundaryConditionType end
+abstract type AbstractBoundaryCondition end
+struct NaturalBC <: AbstractBoundaryCondition end
+struct DirichletBC <: AbstractBoundaryCondition end
+struct PeriodicBC <: AbstractBoundaryCondition end
 
 """
 struct containing information for first derivatives
@@ -55,7 +55,7 @@ struct FirstDerivativeData{TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat
                                         n::Int64,
                                         igrid_full::Array{Int64,2},
                                         element_data::Array{ElementCoordinates,1},
-                                        boundary_condition::Tbc) where Tbc <: AbstractBoundaryConditionType
+                                        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
@@ -107,7 +107,7 @@ function assemble_1D_operator(QQ1D_local::Array{Float64,3},
                             nelement::Int64,
                             n::Int64,
                             igrid_full::Array{Int64,2},
-                            boundary_condition::Tbc) where Tbc <: AbstractBoundaryConditionType
+                            boundary_condition::Tbc) where Tbc <: AbstractBoundaryCondition
     # create the 1D constructor arrays
     nsparse = (nelement - 1)*(ngrid^2 - 1) + ngrid^2
     II = zeros(Int64,nsparse)
@@ -218,7 +218,7 @@ end
 """
 structure containing basic information related to coordinates
 """
-struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
+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
@@ -280,7 +280,7 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
         weight_function::TF=((x)-> 1.0 ),
         # which boundary condition to store and use for first derivatives
         bc::Tbc=NaturalBC()
-        ) where {TF <: Function, Tbc <: AbstractBoundaryConditionType}
+        ) where {TF <: Function, Tbc <: AbstractBoundaryCondition}
         ngrid = scalar_input.ngrid
         nelement = scalar_input.nelement
         # initialise the data used to construct the grid
@@ -328,7 +328,7 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryConditionType}
         weight_function::TF=((x)-> 1.0 ),
         # which boundary condition to store and use in first derivatives
         bc::Tbc=NaturalBC()
-        ) where {TF <: Function, Tbc <: AbstractBoundaryConditionType}
+        ) where {TF <: Function, Tbc <: AbstractBoundaryCondition}
         if typeof(element_data) == Nothing
             # this is a trivial coordinate of length 1
             nelement = 1
@@ -833,7 +833,7 @@ External method for assembling operators in one coordinate
 function assemble_operator(weak_form::Array{Float64,3},
                             x::FiniteElementCoordinate,
                             boundary_condition_x::Tbc
-    ) where Tbc <: AbstractBoundaryConditionType
+    ) where Tbc <: AbstractBoundaryCondition
     return assemble_1D_operator(weak_form,
             x.ngrid, x.nelement, x.n, x.igrid_full,
             boundary_condition_x)
@@ -848,8 +848,8 @@ function assemble_operator(weak_form::TF,
                     boundary_condition_x::Tbcx,
                     boundary_condition_y::Tbcy
                     ) where {TF <: Function,
-                        Tbcx <: AbstractBoundaryConditionType,
-                        Tbcy <: AbstractBoundaryConditionType}
+                        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

From 8e427eff491e8c5d36e8681bdc36bbce32b547b7 Mon Sep 17 00:00:00 2001
From: Michael Hardman <29800382+mrhardman@users.noreply.github.com>
Date: Wed, 4 Mar 2026 11:46:09 +0000
Subject: [PATCH 6/6] Update README.md

---
 README.md | 26 +++++++++++++++-----------
 1 file changed, 15 insertions(+), 11 deletions(-)

diff --git a/README.md b/README.md
index 6185960..aaf9f9c 100644
--- a/README.md
+++ b/README.md
@@ -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}`.
 
@@ -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`,
 ```