--- a/src/BoundaryConditions/BoundaryConditions.jl	Sat Jun 01 17:39:54 2024 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,22 +0,0 @@
-module BoundaryConditions
-
-# REVIEW: Does this need to be in a separate module? I feel like it fits quite well into SbpOperators.
-
-export BoundaryCondition
-export discretize_data
-export boundary_data
-export boundary
-
-export NeumannCondition
-export DirichletCondition
-
-export sat
-export sat_tensors
-
-using Sbplib.Grids
-using Sbplib.LazyTensors
-
-include("boundary_condition.jl")
-include("sat.jl")
-
-end # module

--- a/src/BoundaryConditions/boundary_condition.jl	Sat Jun 01 17:39:54 2024 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,48 +0,0 @@
-"""
-    BoundaryCondition{BID}
-
-A type for implementing boundary_data needed in order to impose a boundary condition.
-Subtypes refer to perticular types of boundary conditions, e.g. Neumann conditions.
-"""
-abstract type BoundaryCondition{BID} end
-
-"""
-    boundary(::BoundaryCondition)
-
-The boundary identifier of the BoundaryCondition.
-"""
-boundary(::BoundaryCondition{BID}) where {BID} = BID()
-
-"""
-    boundary_data(::BoundaryCondition)
-
-If implemented, the data associated with the BoundaryCondition
-"""
-function boundary_data end
-
-"""
-    discretize(grid, bc::BoundaryCondition)
-
-The grid function obtained from discretizing the `bc` boundary_data on the boundary grid
-    specified the by bc `id`.
-"""
-function discretize_data(grid, bc::BoundaryCondition)
-    return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc))
-end
-
-struct DirichletCondition{T1,T2} <: BoundaryCondition{T2}
-    data::T1
-    function DirichletCondition(data, id)
-        return new{typeof(data),typeof(id)}(data)
-    end
-end
-boundary_data(bc::DirichletCondition) = bc.data
-
-struct NeumannCondition{T1,T2} <: BoundaryCondition{T2}
-    data::T1
-    function NeumannCondition(data, id)
-        return new{typeof(data),typeof(id)}(data)
-    end
-end
-boundary_data(bc::NeumannCondition) = bc.data
-

--- a/src/BoundaryConditions/sat.jl	Sat Jun 01 17:39:54 2024 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-"""
-    sat_tensors(op, grid, bc::BoundaryCondition, params...)
-
-The tensor and boundary operator used to construct a simultaneous-approximation-term
-for imposing `bc` related to `op`.
-
-For `penalty_tensor, L  = sat_tensors(...)` then `SAT = penalty_tensor*(L*u - g)`  where `g` 
-is the boundary data.
-"""
-function sat_tensors end
-
-
-"""
-    sat(op, grid, bc::BoundaryCondition; kwargs...)
-
-Simultaneous-Approximation-Term for a general `bc` to `op`. 
-Returns a function `SAT(u,g)` weakly imposing `bc` when added to `op*u`.
-
-`op` must implement the function `sat_tensors`.
-"""
-function sat(op, grid, bc::BoundaryCondition; kwargs...)
-    penalty_tensor, L = sat_tensors(op, grid, bc; kwargs...)
-    return SAT(u, g) = penalty_tensor*(L*u - g)
-end

--- a/src/SbpOperators/SbpOperators.jl	Sat Jun 01 17:39:54 2024 -0700
+++ b/src/SbpOperators/SbpOperators.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -23,21 +23,34 @@
 export second_derivative
 export second_derivative_variable
 export undivided_skewed04
-
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
-using Sbplib.Grids
-using Sbplib.BoundaryConditions
+export closure_size
 
 @enum Parity begin
     odd = -1
     even = 1
 end
 
-export closure_size
 
+# Boundary conditions
+export BoundaryCondition
+export NeumannCondition
+export DirichletCondition
+export discretize_data
+export boundary_data
+export boundary
+export sat
+export sat_tensors
+
+# Using
+using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+using Sbplib.Grids
+
+# Includes
 include("stencil.jl")
 include("stencil_set.jl")
+include("boundary_conditions/boundary_condition.jl")
+include("boundary_conditions/sat.jl")
 include("volumeops/volume_operator.jl")
 include("volumeops/stencil_operator_distinct_closures.jl")
 include("volumeops/constant_interior_scaling_operator.jl")

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundary_conditions/boundary_condition.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -0,0 +1,48 @@
+"""
+    BoundaryCondition{BID}
+
+A type for implementing boundary_data needed in order to impose a boundary condition.
+Subtypes refer to perticular types of boundary conditions, e.g. Neumann conditions.
+"""
+abstract type BoundaryCondition{BID} end
+
+"""
+    boundary(::BoundaryCondition)
+
+The boundary identifier of the BoundaryCondition.
+"""
+boundary(::BoundaryCondition{BID}) where {BID} = BID()
+
+"""
+    boundary_data(::BoundaryCondition)
+
+If implemented, the data associated with the BoundaryCondition
+"""
+function boundary_data end
+
+"""
+    discretize(grid, bc::BoundaryCondition)
+
+The grid function obtained from discretizing the `bc` boundary_data on the boundary grid
+    specified the by bc `id`.
+"""
+function discretize_data(grid, bc::BoundaryCondition)
+    return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc))
+end
+
+struct DirichletCondition{T1,T2} <: BoundaryCondition{T2}
+    data::T1
+    function DirichletCondition(data, id)
+        return new{typeof(data),typeof(id)}(data)
+    end
+end
+boundary_data(bc::DirichletCondition) = bc.data
+
+struct NeumannCondition{T1,T2} <: BoundaryCondition{T2}
+    data::T1
+    function NeumannCondition(data, id)
+        return new{typeof(data),typeof(id)}(data)
+    end
+end
+boundary_data(bc::NeumannCondition) = bc.data
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundary_conditions/sat.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -0,0 +1,24 @@
+"""
+    sat_tensors(op, grid, bc::BoundaryCondition, params...)
+
+The tensor and boundary operator used to construct a simultaneous-approximation-term
+for imposing `bc` related to `op`.
+
+For `penalty_tensor, L  = sat_tensors(...)` then `SAT = penalty_tensor*(L*u - g)`  where `g` 
+is the boundary data.
+"""
+function sat_tensors end
+
+
+"""
+    sat(op, grid, bc::BoundaryCondition; kwargs...)
+
+Simultaneous-Approximation-Term for a general `bc` to `op`. 
+Returns a function `SAT(u,g)` weakly imposing `bc` when added to `op*u`.
+
+`op` must implement the function `sat_tensors`.
+"""
+function sat(op, grid, bc::BoundaryCondition; kwargs...)
+    penalty_tensor, L = sat_tensors(op, grid, bc; kwargs...)
+    return SAT(u, g) = penalty_tensor*(L*u - g)
+end

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sat Jun 01 17:39:54 2024 -0700
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -62,7 +62,7 @@
 
 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
 """
-function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning = (1., 1.))
+function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning = (1., 1.))
     id = boundary(bc)
     set  = Δ.stencil_set
     H⁻¹ = inverse_inner_product(g,set)
@@ -82,7 +82,7 @@
 
 See also: [`sat`,`NeumannCondition`](@ref).
 """
-function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
     id = boundary(bc)
     set  = Δ.stencil_set
     H⁻¹ = inverse_inner_product(g,set)

--- a/src/Sbplib.jl	Sat Jun 01 17:39:54 2024 -0700
+++ b/src/Sbplib.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -4,7 +4,6 @@
 include("RegionIndices/RegionIndices.jl")
 include("LazyTensors/LazyTensors.jl")
 include("Grids/Grids.jl")
-include("BoundaryConditions/BoundaryConditions.jl")
 include("SbpOperators/SbpOperators.jl")
 include("DiffOps/DiffOps.jl")
 
@@ -13,6 +12,5 @@
 export Grids
 export SbpOperators
 export DiffOps
-export BoundaryConditions
 
 end

--- a/test/BoundaryConditions/boundary_condition_test.jl	Sat Jun 01 17:39:54 2024 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-using Test
-
-using Sbplib.BoundaryConditions
-using Sbplib.Grids
-using Sbplib.RegionIndices
-
-@testset "BoundaryCondition" begin
-    grid_1d = equidistant_grid(0.0, 1.0, 11)
-    grid_2d = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 15)
-    grid_3d = equidistant_grid((0.0, 0.0, 0.0), (1.0,1.0, 1.0), 11, 15, 13)
-    (id_l,_) = boundary_identifiers(grid_1d)
-    (_,_,_,id_n) = boundary_identifiers(grid_2d)
-    (_,_,_,_,id_b,_) = boundary_identifiers(grid_3d)
-
-    g = 3.14
-    f(x,y,z) = x^2+y^2+z^2
-    @testset "Constructors" begin
-        @test DirichletCondition(g,id_l) isa BoundaryCondition{Lower}
-        @test DirichletCondition(g,id_n) isa BoundaryCondition{CartesianBoundary{2,Upper}}
-        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower}
-        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function}
-    end
-
-    @testset "boundary" begin
-        @test boundary(DirichletCondition(g,id_l)) == id_l
-        @test boundary(NeumannCondition(f,id_b)) == id_b
-    end
-
-    @testset "boundary_data" begin
-        @test boundary_data(DirichletCondition(g,id_l)) == g
-        @test boundary_data(NeumannCondition(f,id_b)) == f
-    end
-
-    @testset "discretize_data" begin
-        @test fill(g) ≈ discretize_data(grid_1d,DirichletCondition(g,id_l))
-        @test g*ones(11,1) ≈ discretize_data(grid_2d,DirichletCondition(g,id_n))
-        X = repeat(0:1/10:1, inner = (1,15))
-        Y = repeat(0:1/14:1, outer = (1,11))
-        @test map((x,y)->f(x,y,0), X,Y') ≈ discretize_data(grid_3d,NeumannCondition(f,id_b))
-    end
-end

--- a/test/BoundaryConditions/sat_test.jl	Sat Jun 01 17:39:54 2024 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,73 +0,0 @@
-using Test
-
-
-using Sbplib.BoundaryConditions
-using Sbplib.Grids
-using Sbplib.LazyTensors
-using Sbplib.SbpOperators
-
-stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
-
-struct MockOp end
-
-function BoundaryConditions.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition; a = 1.)
-    e = boundary_restriction(g, stencil_set, boundary(bc))
-    L = a*e
-    sat_op = e'
-    return sat_op, L
-end
-
-function BoundaryConditions.sat_tensors(op::MockOp, g::Grid, bc::NeumannCondition)
-    e = boundary_restriction(g, stencil_set, boundary(bc))
-    d = normal_derivative(g, stencil_set, boundary(bc))
-    L = d
-    sat_op = e'
-    return sat_op, L
-end
-
-@testset "sat" begin
-    op = MockOp()
-    @testset "1D" begin
-        grid  = equidistant_grid(0., 1., 11)
-        l, r = boundary_identifiers(grid)
-        u = eval_on(grid, x-> 1. + 2x^2)
-        dc = DirichletCondition(1.0, l)
-        g_l = discretize_data(grid, dc)
-        SAT_l = sat(op, grid, dc)
-        @test SAT_l(u, g_l) ≈ zeros((size(grid))) atol = 1e-13
-        
-        nc = NeumannCondition(4.0, r)
-        g_r = discretize_data(grid, nc)
-        SAT_r = sat(op, grid, nc)
-        @test SAT_r(u, g_r) ≈ zeros((size(grid))) atol = 1e-13
-    end
-    @testset "2D" begin
-        grid  = equidistant_grid((0.,0.), (1.,1.), 11, 13)
-        W, E, S, N = boundary_identifiers(grid)
-        u = eval_on(grid, (x,y) -> x+y^2)
-
-        dc_W = DirichletCondition(1.0, W)
-        SAT_W = sat(op, grid, dc_W)
-        g_W = discretize_data(grid, dc_W)
-        r_W = zeros(size(grid))
-        r_W[1,:] .= map(y -> (y^2-1.), range(0., 1., length=13))
-        @test SAT_W(u, g_W) ≈ r_W atol = 1e-13
-
-        dc_E = DirichletCondition(2, E)
-        SAT_E = sat(op, grid, dc_E; a = 2.)
-        g_E = discretize_data(grid, dc_E)
-        r_E = zeros(size(grid))
-        r_E[end,:] .= map(y -> (2*(1. + y^2)-2.), range(0., 1., length=13))
-        @test SAT_E(u, g_E) ≈ r_E atol = 1e-13
-
-        nc_S = NeumannCondition(.0, S)
-        SAT_S = sat(op, grid, nc_S)
-        g_S = discretize_data(grid, nc_S)
-        @test SAT_S(u, g_S) ≈ zeros(size(grid)) atol = 1e-13
-
-        nc_N = NeumannCondition(2.0, N)
-        SAT_N = sat(op, grid, nc_N)
-        g_N = discretize_data(grid, nc_N)
-        @test SAT_N(u, g_N) ≈ zeros(size(grid)) atol = 1e-13
-    end
-end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -0,0 +1,41 @@
+using Test
+
+using Sbplib.Grids
+using Sbplib.RegionIndices
+using Sbplib.SbpOperators
+
+@testset "BoundaryCondition" begin
+    grid_1d = equidistant_grid(0.0, 1.0, 11)
+    grid_2d = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 15)
+    grid_3d = equidistant_grid((0.0, 0.0, 0.0), (1.0,1.0, 1.0), 11, 15, 13)
+    (id_l,_) = boundary_identifiers(grid_1d)
+    (_,_,_,id_n) = boundary_identifiers(grid_2d)
+    (_,_,_,_,id_b,_) = boundary_identifiers(grid_3d)
+
+    g = 3.14
+    f(x,y,z) = x^2+y^2+z^2
+    @testset "Constructors" begin
+        @test DirichletCondition(g,id_l) isa BoundaryCondition{Lower}
+        @test DirichletCondition(g,id_n) isa BoundaryCondition{CartesianBoundary{2,Upper}}
+        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower}
+        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function}
+    end
+
+    @testset "boundary" begin
+        @test boundary(DirichletCondition(g,id_l)) == id_l
+        @test boundary(NeumannCondition(f,id_b)) == id_b
+    end
+
+    @testset "boundary_data" begin
+        @test boundary_data(DirichletCondition(g,id_l)) == g
+        @test boundary_data(NeumannCondition(f,id_b)) == f
+    end
+
+    @testset "discretize_data" begin
+        @test fill(g) ≈ discretize_data(grid_1d,DirichletCondition(g,id_l))
+        @test g*ones(11,1) ≈ discretize_data(grid_2d,DirichletCondition(g,id_n))
+        X = repeat(0:1/10:1, inner = (1,15))
+        Y = repeat(0:1/14:1, outer = (1,11))
+        @test map((x,y)->f(x,y,0), X,Y') ≈ discretize_data(grid_3d,NeumannCondition(f,id_b))
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/boundary_conditions/sat_test.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -0,0 +1,71 @@
+using Test
+
+using Sbplib.Grids
+using Sbplib.LazyTensors
+using Sbplib.SbpOperators
+
+stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
+
+struct MockOp end
+
+function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition; a = 1.)
+    e = boundary_restriction(g, stencil_set, boundary(bc))
+    L = a*e
+    sat_op = e'
+    return sat_op, L
+end
+
+function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::NeumannCondition)
+    e = boundary_restriction(g, stencil_set, boundary(bc))
+    d = normal_derivative(g, stencil_set, boundary(bc))
+    L = d
+    sat_op = e'
+    return sat_op, L
+end
+
+@testset "sat" begin
+    op = MockOp()
+    @testset "1D" begin
+        grid  = equidistant_grid(0., 1., 11)
+        l, r = boundary_identifiers(grid)
+        u = eval_on(grid, x-> 1. + 2x^2)
+        dc = DirichletCondition(1.0, l)
+        g_l = discretize_data(grid, dc)
+        SAT_l = sat(op, grid, dc)
+        @test SAT_l(u, g_l) ≈ zeros((size(grid))) atol = 1e-13
+        
+        nc = NeumannCondition(4.0, r)
+        g_r = discretize_data(grid, nc)
+        SAT_r = sat(op, grid, nc)
+        @test SAT_r(u, g_r) ≈ zeros((size(grid))) atol = 1e-13
+    end
+    @testset "2D" begin
+        grid  = equidistant_grid((0.,0.), (1.,1.), 11, 13)
+        W, E, S, N = boundary_identifiers(grid)
+        u = eval_on(grid, (x,y) -> x+y^2)
+
+        dc_W = DirichletCondition(1.0, W)
+        SAT_W = sat(op, grid, dc_W)
+        g_W = discretize_data(grid, dc_W)
+        r_W = zeros(size(grid))
+        r_W[1,:] .= map(y -> (y^2-1.), range(0., 1., length=13))
+        @test SAT_W(u, g_W) ≈ r_W atol = 1e-13
+
+        dc_E = DirichletCondition(2, E)
+        SAT_E = sat(op, grid, dc_E; a = 2.)
+        g_E = discretize_data(grid, dc_E)
+        r_E = zeros(size(grid))
+        r_E[end,:] .= map(y -> (2*(1. + y^2)-2.), range(0., 1., length=13))
+        @test SAT_E(u, g_E) ≈ r_E atol = 1e-13
+
+        nc_S = NeumannCondition(.0, S)
+        SAT_S = sat(op, grid, nc_S)
+        g_S = discretize_data(grid, nc_S)
+        @test SAT_S(u, g_S) ≈ zeros(size(grid)) atol = 1e-13
+
+        nc_N = NeumannCondition(2.0, N)
+        SAT_N = sat(op, grid, nc_N)
+        g_N = discretize_data(grid, nc_N)
+        @test SAT_N(u, g_N) ≈ zeros(size(grid)) atol = 1e-13
+    end
+end

--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sat Jun 01 17:39:54 2024 -0700
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Jun 04 16:46:14 2024 -0700
@@ -3,7 +3,6 @@
 using Sbplib.SbpOperators
 using Sbplib.Grids
 using Sbplib.LazyTensors
-using Sbplib.BoundaryConditions
 
 @testset "Laplace" begin
     # Default stencils (4th order)