Ensure that first derivative operators within the FiniteElementCoordinate struct can only see either a NaturalBC or a PeriodicBC. Create a function interface that returns a NaturalBC if a DirichletBC is stored in the FiniteElementCoordinate struct.

Author: mrhardman

src/FiniteElementAssembly.jl

@@ -28,6 +28,16 @@ struct NaturalBC <: AbstractBoundaryCondition end
 struct DirichletBC <: AbstractBoundaryCondition end
 struct PeriodicBC <: AbstractBoundaryCondition end
 
+# logic to ensure that only NaturalBC or PeriodicBC are used in first derivatives
+# - periodic boundary conditions are passed on
+function first_derivative_boundary_condition(bc::PeriodicBC)
+    return PeriodicBC()
+end
+# - Dirichlet boundary conditions are converted to natural boundary conditions
+function first_derivative_boundary_condition(bc::Union{NaturalBC,DirichletBC})
+    return NaturalBC()
+end
+
 """
 struct containing information for first derivatives
 """
@@ -69,12 +79,13 @@ struct FirstDerivativeData{TFloat <: Real, TMatrix <: AbstractSparseArray{TFloat
         # dummy arrays
         dummy_rhs = Array{Float64,1,}(undef,n)
         dummy_df = Array{Float64,1,}(undef,n)
-
+        # boundary condition for first derivative matrices, only NaturalBC or PeriodicBC allowed
+        first_derivative_bc = first_derivative_boundary_condition(boundary_condition)
         # assemble the global matrices with a sparse construction
         MM1D = assemble_1D_operator(MM, ngrid, nelement,
-                            n, igrid_full, boundary_condition)
+                            n, igrid_full, first_derivative_bc)
         PP1D = assemble_1D_operator(PP, ngrid, nelement,
-                            n, igrid_full, boundary_condition)
+                            n, igrid_full, first_derivative_bc)
         MM1D_sparse = sparse(MM1D)
         PP1D_sparse = sparse(PP1D)
         lu_MM1D = lu(MM1D_sparse)

test/CalculusTests.jl

@@ -9,7 +9,8 @@ using FiniteElementAssembly: FiniteElementCoordinate,
                             first_derivative!,
                             ScalarCoordinateInputs,
                             include_boundary_points, exclude_boundary_points,
-                            exclude_lower_boundary_point, exclude_upper_boundary_point
+                            exclude_lower_boundary_point, exclude_upper_boundary_point,
+                            NaturalBC, DirichletBC
 using FiniteElementAssembly: integral
 
 """
@@ -24,16 +25,18 @@ function runtests()
         rng = StableRNG(42)
 
         @testset "GaussLegendre derivatives and fundamental theorem of calculus, testing exact polynomials" begin
-            @testset "$nelement $ngrid $boundary_points" for nelement ∈ (1:5), ngrid ∈ (3:17),
+            @testset "$nelement $ngrid $boundary_points $bc" for nelement ∈ (1:5), ngrid ∈ (3:17),
                     boundary_points in (include_boundary_points, exclude_boundary_points,
-                            exclude_lower_boundary_point, exclude_upper_boundary_point)
+                            exclude_lower_boundary_point, exclude_upper_boundary_point),
+                    bc in (NaturalBC(), DirichletBC())
 
                 # define inputs needed for the test
                 coord_min = -1.0
                 coord_max = 1.0
                 # create the coordinate struct 'x'
                 x = FiniteElementCoordinate("coord", ScalarCoordinateInputs(ngrid,
-                                      nelement, coord_min, coord_max, boundary_points))
+                                      nelement, coord_min, coord_max, boundary_points),
+                                      bc=bc)
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
                     # create array for the function f(x) to be differentiated/integrated