--- a/src/BoundaryConditions/BoundaryConditions.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/BoundaryConditions/BoundaryConditions.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -4,8 +4,8 @@
 
 export BoundaryCondition
 export discretize_data
-export data # REVIEW: Name too generic to export.
-export id # REVIEW: Also to generic. Change to `boundary`?
+export boundary_data
+export boundary
 
 export NeumannCondition
 export DirichletCondition

--- a/src/BoundaryConditions/boundary_condition.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/BoundaryConditions/boundary_condition.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -1,47 +1,48 @@
 """
-    BoundaryCondition
+    BoundaryCondition{BID}
 
-A type for implementing data needed in order to impose a boundary condition.
+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{T1,T2} end # REVIEW: No type parameters needed here.
+abstract type BoundaryCondition{BID} end
 
 """
-    id(::BoundaryCondition)
+    boundary(::BoundaryCondition)
 
 The boundary identifier of the BoundaryCondition.
-Must be implemented by subtypes.
 """
-function id end
+boundary(::BoundaryCondition{BID}) where {BID} = BID()
 
 """
-    data(::BoundaryCondition)
+    boundary_data(::BoundaryCondition)
 
 If implemented, the data associated with the BoundaryCondition
 """
-function data end
+function boundary_data end
 
 """
     discretize(grid, bc::BoundaryCondition)
 
-The grid function obtained from discretizing the `bc` data on the boundary grid
+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, id(bc)), data(bc))
+    return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc))
 end
 
-struct DirichletCondition{T1,T2} <: BoundaryCondition{T1,T2}
+struct DirichletCondition{T1,T2} <: BoundaryCondition{T2}
     data::T1
-    id::T2 # REVIEW: This field not needed since BoundaryId are usually type parameters?
+    function DirichletCondition(data, id)
+        return new{typeof(data),typeof(id)}(data)
+    end
 end
-id(bc::DirichletCondition) = bc.id
-data(bc::DirichletCondition) = bc.data
+boundary_data(bc::DirichletCondition) = bc.data
 
-struct NeumannCondition{T1,T2} <: BoundaryCondition{T1,T2}
+struct NeumannCondition{T1,T2} <: BoundaryCondition{T2}
     data::T1
-    id::T2 # REVIEW: This field not needed since BoundaryId are usually type parameters?
+    function NeumannCondition(data, id)
+        return new{typeof(data),typeof(id)}(data)
+    end
 end
-id(bc::NeumannCondition) = bc.id
-data(bc::NeumannCondition) = bc.data
+boundary_data(bc::NeumannCondition) = bc.data
 

--- a/src/BoundaryConditions/sat.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/BoundaryConditions/sat.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -4,22 +4,21 @@
 The tensor and boundary operator used to construct a simultaneous-approximation-term
 for imposing `bc` related to `op`.
 
-For `sat_op, L  = sat_tensors(...)` then `SAT = sat_op*(L*u - g)`  where `g` 
+For `penalty_tensor, L  = sat_tensors(...)` then `SAT = penalty_tensor*(L*u - g)`  where `g` 
 is the boundary data.
 """
 function sat_tensors end
-# REVIEW: Rename `sat_op` to `penalty_tensor`?
 
 
 """
-    sat(op, grid, bc::BoundaryCondition, params...)
+    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, params...)
-    sat_op, L = sat_tensors(op, grid, bc, params...)
-    return SAT(u, g) = sat_op*(L*u - g)
+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/Grids/Grids.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/Grids/Grids.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -18,7 +18,6 @@
 export eval_on
 export componentview
 export ArrayComponentView
-export orthogonal_grid
 
 export BoundaryIdentifier
 export TensorGridBoundary

--- a/src/Grids/equidistant_grid.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -52,7 +52,6 @@
 boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
 boundary_indices(g::EquidistantGrid, id::Lower) = (1,)
 boundary_indices(g::EquidistantGrid, id::Upper) = (length(g),)
-orthogonal_grid(g::EquidistantGrid, id::Union{Lower,Upper}) = g
 
 """
     refine(g::EquidistantGrid, r::Int)

--- a/src/Grids/grid.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/Grids/grid.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -127,15 +127,6 @@
 """
 function boundary_indices end
 
-"""
-    orthogonal_grid(g::Grid, id::BoundaryIdentifier)
-
-The grid for the coordinate direction orthogonal to that of the boundary
-with given id
-"""
-function orthogonal_grid end
-
-function boundary_indices end
 
 """
     eval_on(g::Grid, f)

--- a/src/Grids/tensor_grid.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/Grids/tensor_grid.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -95,8 +95,6 @@
     return LazyTensors.concatenate_tuples(b_ind...)
 end
 
-orthogonal_grid(g::TensorGrid, id::BoundaryIdentifier) = g.grids[grid_id(id)] # REVIEW: Seems clearer to me to just inline this and remove the function definition?
-
 
 function combined_coordinate_vector_type(coordinate_types...)
     combined_coord_length = mapreduce(_ncomponents, +, coordinate_types)

--- a/src/SbpOperators/operators/standard_diagonal.toml	Wed May 29 23:28:33 2024 +0200
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Sat Jun 01 17:39:54 2024 -0700
@@ -39,7 +39,7 @@
         {s = [["2", "-1", "0"],["-3", "1",   "0"],["1","0","0"]], c = 1},
 ]
 
-Positivity.D2 = {theta_H = "1/2", theta_M = "0.3636363636", theta_R = "1.000000538455350", m_b = "2"}
+D2.positivity = {theta_M = "0.3636363636", theta_R = "1.000000538455350", m_b = "2"}
 
 [[stencil_set]]
 
@@ -137,7 +137,7 @@
     ]}
 ]
 
-Positivity.D2 = {theta_H = "17/48", theta_M = "0.2505765857", theta_R = "0.577587500088313", m_b = "4"} # REVIEW: Shouldn't the keys be reversed here (D2.positivity instead of Positivity.D2)
+D2.positivity = {theta_M = "0.2505765857", theta_R = "0.577587500088313", m_b = "4"}
 
 [[stencil_set]]
 
@@ -153,4 +153,4 @@
 e.closure = ["1"]
 d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"]
 
-Positivity.D2 = {theta_H = "13649/43200", theta_M = "0.1878687080", theta_R = "0.3697", m_b = "7"}
+D2.positivity = {theta_M = "0.1878687080", theta_R = "0.3697", m_b = "7"}

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Wed May 29 23:28:33 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -53,28 +53,26 @@
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
-# REVIEW: I think the handling of tuning parameters below should be through kwargs instead.
 
 """
-sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning)
 
 The operators required to construct the SAT for imposing a Dirichlet condition.
 `tuning` specifies the strength of the penalty. See
 
 See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
 """
-function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
-    id = bc.id
+function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning = (1., 1.))
+    id = boundary(bc)
     set  = Δ.stencil_set
     H⁻¹ = inverse_inner_product(g,set)
     Hᵧ = inner_product(boundary_grid(g, id), set)
     e = boundary_restriction(g, set, id)
     d = normal_derivative(g, set, id)
     B = positivity_decomposition(Δ, g, bc, tuning)
-    sat_op = H⁻¹∘(d' - B*e')∘Hᵧ
-    return sat_op, e
+    penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
+    return penalty_tensor, e
 end
-BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition) = BoundaryConditions.sat_tensors(Δ, g, bc, (1.,1.)) # REVIEW: Should be possible to replace this with argument default values.
 
 """
 sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
@@ -85,21 +83,23 @@
 See also: [`sat`,`NeumannCondition`](@ref).
 """
 function BoundaryConditions.sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
-    id = bc.id
+    id = boundary(bc)
     set  = Δ.stencil_set
     H⁻¹ = inverse_inner_product(g,set)
     Hᵧ = inner_product(boundary_grid(g, id), set)
     e = boundary_restriction(g, set, id)
     d = normal_derivative(g, set, id)
 
-    sat_op = -H⁻¹∘e'∘Hᵧ
-    return sat_op, d
+    penalty_tensor = -H⁻¹∘e'∘Hᵧ
+    return penalty_tensor, d
 end
 
-# REVIEW: This function assumes a TensorGrid right? In that case there should probably be a type annotation to get clearer error messages.
-function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition, tuning)
+# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
+# change bc::BoundaryCondition to id::BoundaryIdentifier
+
+function positivity_decomposition(Δ::Laplace, g::EquidistantGrid, bc::BoundaryCondition, tuning)
     pos_prop = positivity_properties(Δ)
-    h = spacing(orthogonal_grid(g, bc.id))
+    h = spacing(g)
     θ_H = pos_prop.theta_H
     τ_H = tuning[1]*ndims(g)/(h*θ_H)
     θ_R = pos_prop.theta_R
@@ -108,4 +108,19 @@
     return B
 end
 
-positivity_properties(Δ::Laplace) = parse_named_tuple(Δ.stencil_set["Positivity"]["D2"]) # REVIEW: Can this function extract theta_H from the inner product instead of storing it twice in the TOML?
+function positivity_decomposition(Δ::Laplace, g::TensorGrid, bc::BoundaryCondition, tuning)
+    pos_prop = positivity_properties(Δ)
+    h = spacing(g.grids[grid_id(boundary(bc))]) # grid spacing of the 1D grid normal to the boundary
+    θ_H = pos_prop.theta_H
+    τ_H = tuning[1]*ndims(g)/(h*θ_H)
+    θ_R = pos_prop.theta_R
+    τ_R = tuning[2]/(h*θ_R)
+    B = τ_H + τ_R
+    return B
+end
+
+function positivity_properties(Δ::Laplace)
+    D2_pos_prop = parse_named_tuple(Δ.stencil_set["D2"]["positivity"])
+    H_closure = parse_tuple(Δ.stencil_set["H"]["closure"])
+    return merge(D2_pos_prop, (theta_H = H_closure[1],))
+end

--- a/test/BoundaryConditions/boundary_condition_test.jl	Wed May 29 23:28:33 2024 +0200
+++ b/test/BoundaryConditions/boundary_condition_test.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -2,6 +2,7 @@
 
 using Sbplib.BoundaryConditions
 using Sbplib.Grids
+using Sbplib.RegionIndices
 
 @testset "BoundaryCondition" begin
     grid_1d = equidistant_grid(0.0, 1.0, 11)
@@ -13,13 +14,28 @@
 
     g = 3.14
     f(x,y,z) = x^2+y^2+z^2
-    @test DirichletCondition(g,id_l) isa BoundaryCondition{Float64}
-    @test DirichletCondition(g,id_n) isa BoundaryCondition{Float64}
-    @test NeumannCondition(f,id_b) isa BoundaryCondition{<:Function}
+    @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
 
-    @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))
+    @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	Wed May 29 23:28:33 2024 +0200
+++ b/test/BoundaryConditions/sat_test.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -10,16 +10,16 @@
 
 struct MockOp end
 
-function BoundaryConditions.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition)
-    e = boundary_restriction(g, stencil_set, id(bc))
-    L = e
+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, id(bc))
-    d = normal_derivative(g, stencil_set, id(bc))
+    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
@@ -54,10 +54,10 @@
         @test SAT_W(u, g_W) ≈ r_W atol = 1e-13
 
         dc_E = DirichletCondition(2, E)
-        SAT_E = sat(op, grid, dc_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 -> ((1. + y^2)-2.), range(0., 1., length=13))
+        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)

--- a/test/Grids/equidistant_grid_test.jl	Wed May 29 23:28:33 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -79,12 +79,6 @@
 
     end
 
-    @testset "orthogonal_grid" begin
-        g = EquidistantGrid(0:0.1:1)
-        @test orthogonal_grid(g, Lower()) == g
-        @test orthogonal_grid(g, Upper()) == g
-    end
-
     @testset "refine" begin
         g = EquidistantGrid(0:0.1:1)
         @test refine(g, 1) == g

--- a/test/Grids/tensor_grid_test.jl	Wed May 29 23:28:33 2024 +0200
+++ b/test/Grids/tensor_grid_test.jl	Sat Jun 01 17:39:54 2024 -0700
@@ -185,16 +185,6 @@
         @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Upper}()) == (11,)
     end
 
-    @testset "orthogonal_grid" begin
-        g₁ = EquidistantGrid(range(0,1,length=11))
-        g₂ = EquidistantGrid(range(2,3,length=6))
-        g₃ = EquidistantGrid(range(4,5,length=7))
-    
-        @test orthogonal_grid(TensorGrid(g₁, g₂, g₃), TensorGridBoundary{1, Lower}()) == g₁
-        @test orthogonal_grid(TensorGrid(g₁, g₂, g₃), TensorGridBoundary{2, Lower}()) == g₂
-        @test orthogonal_grid(TensorGrid(g₁, g₂, g₃), TensorGridBoundary{3, Upper}()) == g₃
-    end
-
 end
 
 @testset "combined_coordinate_vector_type" begin