--- a/Notes.md	Sat Jun 29 17:06:27 2024 +0200
+++ b/Notes.md	Sun Jul 07 15:15:12 2024 -0700
@@ -129,8 +129,24 @@
  - [ ] Figure out how to treat the borrowing parameters of operators. Include in into the struct? Expose via function dispatched on the operator type and grid?
 
 ## Identifiers for regions
-The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probably be included in the grid module. This allows new grid types to come with their own regions.
-We implement this by refactoring RegionIndices to be agnostic to the region types and then moving the actual types to Grids.
+The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probably be included in the grid module. This allows new grid types to come with their own regions. We implement this by refactoring RegionIndices to be agnostic to the region types and then moving the actual types to Grids.
+
+`Region`s denote the parts of an operator for which closures and interior stencils are applied respectively. Boundaries/interfaces are included in the closure regions but or are separate concepts. The currently used `Upper`, `Lower` should perhaps instead be subtypes of `BoundaryIdentifier`. Perhaps we should also change the name of `Region` to `StencilRegion`.
+
+Dispatching on regions:
+ * Use region indices `Index{Region}` as presently. This is simply an integer parametrized on the region. 
+ * Use tag dispatch in the apply functions, i.e. `apply(op,u,I,::Region)` where `I` now is a regular `CartesianIndex`. Drawback with this approch is that regular indexing (op*u)[I] won't be aware of regions.
+
+The structure of the apply function may then either be (using the tag dispatch alternative)
+```julia
+function apply(op,u,I,r::Region)
+    if r == Lower()
+        return apply_lower(op,u,I)
+    ...
+    end
+end
+```
+with specific implementations of `apply_lower`. Here we need to make sure that the top level apply function actually is removed by the compiler. This should be possible since `r` is known in compile time. Alternatively one directly specializes `apply(op,u,I,::Lower)` and then introduces a toplevel function `apply(op,u,I) = apply(op,u,I,get_region(op,I))`. The second alterative seems more Julian to me.
 
 ## Regions and tensormappings
 - [ ] Use a trait to indicate if a LazyTensor uses indices with regions.

--- a/benchmark/benchmarks.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/benchmark/benchmarks.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -3,7 +3,6 @@
 using Sbplib
 using Sbplib.Grids
 using Sbplib.SbpOperators
-using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 
 using LinearAlgebra
@@ -168,20 +167,20 @@
 Dxx = second_derivative(g2, stencil_set, 1)
 Dyy = second_derivative(g2, stencil_set, 2)
 
-e₁ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,Lower}())
-e₁ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,Upper}())
-e₂ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,Lower}())
-e₂ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,Upper}())
+e₁ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,LowerBoundary}())
+e₁ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,UpperBoundary}())
+e₂ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,LowerBoundary}())
+e₂ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,UpperBoundary}())
 
-d₁ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{1,Lower}())
-d₁ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{1,Upper}())
-d₂ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{2,Lower}())
-d₂ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{2,Upper}())
+d₁ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{1,LowerBoundary}())
+d₁ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{1,UpperBoundary}())
+d₂ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{2,LowerBoundary}())
+d₂ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{2,UpperBoundary}())
 
-H₁ₗ = inner_product(boundary_grid(g2, CartesianBoundary{1,Lower}()), stencil_set)
-H₁ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{1,Upper}()), stencil_set)
-H₂ₗ = inner_product(boundary_grid(g2, CartesianBoundary{2,Lower}()), stencil_set)
-H₂ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{2,Upper}()), stencil_set)
+H₁ₗ = inner_product(boundary_grid(g2, CartesianBoundary{1,LowerBoundary}()), stencil_set)
+H₁ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{1,UpperBoundary}()), stencil_set)
+H₂ₗ = inner_product(boundary_grid(g2, CartesianBoundary{2,LowerBoundary}()), stencil_set)
+H₂ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{2,UpperBoundary}()), stencil_set)
 
 SUITE["boundary_terms"]["pre_composition"] = @benchmarkable $u2 .= $(H⁻¹∘e₁ₗ'∘H₁ₗ∘d₁ₗ)*$v2
 SUITE["boundary_terms"]["composition"]     = @benchmarkable $u2 .= ($H⁻¹∘$e₁ₗ'∘$H₁ₗ∘$d₁ₗ)*$v2

--- a/src/Grids/Grids.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/src/Grids/Grids.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -1,6 +1,5 @@
 module Grids
 
-using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 using StaticArrays
 
@@ -23,6 +22,8 @@
 export BoundaryIdentifier
 export TensorGridBoundary
 export CartesianBoundary
+export LowerBoundary
+export UpperBoundary
 
 export TensorGrid
 export ZeroDimGrid
@@ -32,9 +33,6 @@
 export spacing
 export equidistant_grid
 
-
-abstract type BoundaryIdentifier end
-
 include("grid.jl")
 include("tensor_grid.jl")
 include("equidistant_grid.jl")

--- a/src/Grids/equidistant_grid.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -49,12 +49,30 @@
 
 min_spacing(g::EquidistantGrid) = spacing(g)
 
+"""
+    LowerBoundary <: BoundaryIdentifier
 
-boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
-boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
-boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
-boundary_indices(g::EquidistantGrid, id::Lower) = (1,)
-boundary_indices(g::EquidistantGrid, id::Upper) = (length(g),)
+Boundary identifier for the the lower (left) boundary of an EquidistantGrid
+
+See also: [`BoundaryIdentifier`](@ref)
+"""
+struct LowerBoundary <: BoundaryIdentifier end
+
+"""
+    UpperBoundary <: BoundaryIdentifier
+
+Boundary identifier for the the upper (right) boundary of an EquidistantGrid
+
+See also: [`BoundaryIdentifier`](@ref)
+"""
+struct UpperBoundary <: BoundaryIdentifier end
+
+
+boundary_identifiers(::EquidistantGrid) = (LowerBoundary(), UpperBoundary())
+boundary_grid(g::EquidistantGrid, id::LowerBoundary) = ZeroDimGrid(g[begin])
+boundary_grid(g::EquidistantGrid, id::UpperBoundary) = ZeroDimGrid(g[end])
+boundary_indices(g::EquidistantGrid, id::LowerBoundary) = (firstindex(g),)
+boundary_indices(g::EquidistantGrid, id::UpperBoundary) = (lastindex(g),)
 
 """
     refine(g::EquidistantGrid, r::Int)

--- a/src/Grids/grid.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/src/Grids/grid.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -109,6 +109,13 @@
 function coarsen end
 
 """
+   BoundaryIdentifier
+
+An identifier for a boundary of a grid.
+"""
+abstract type BoundaryIdentifier end
+
+"""
     boundary_identifiers(g::Grid)
 
 Identifiers for all the boundaries of `g`.

--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -1,27 +1,27 @@
 """
-    BoundaryOperator{T,R,N} <: LazyTensor{T,0,1}
+    BoundaryOperator{T,B,N} <: LazyTensor{T,0,1}
 
 Implements the boundary operator `op` for 1D as a `LazyTensor`
 
 `op` is the restriction of a grid function to the boundary using some closure
-`Stencil{T,N}`. The boundary to restrict to is determined by `R`. `op'` is the
+`Stencil{T,N}`. The boundary to restrict to is determined by `B`. `op'` is the
 prolongation of a zero dimensional array to the whole grid using the same
 closure stencil.
 """
-struct BoundaryOperator{T,R<:Region,N} <: LazyTensor{T,0,1}
+struct BoundaryOperator{T,B<:BoundaryIdentifier,N} <: LazyTensor{T,0,1}
     stencil::Stencil{T,N}
     size::Int
 end
 
 """
-    BoundaryOperator(grid::EquidistantGrid, closure_stencil, region)
+    BoundaryOperator(grid::EquidistantGrid, closure_stencil, boundary)
 
 Constructs the BoundaryOperator with stencil `closure_stencil` for a
 `EquidistantGrid` `grid`, restricting to to the boundary specified by
-`region`.
+`boundary`.
 """
-function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, region::Region) where {T,N}
-    return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1])
+function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, boundary::BoundaryIdentifier) where {T,N}
+    return BoundaryOperator{T,typeof(boundary),N}(closure_stencil,size(grid)[1])
 end
 
 """
@@ -29,24 +29,24 @@
 
 The size of the closure stencil.
 """
-closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N
+closure_size(::BoundaryOperator{T,B,N}) where {T,B,N} = N
 
 LazyTensors.range_size(op::BoundaryOperator) = ()
 LazyTensors.domain_size(op::BoundaryOperator) = (op.size,)
 
-function LazyTensors.apply(op::BoundaryOperator{<:Any,Lower}, v::AbstractVector)
+function LazyTensors.apply(op::BoundaryOperator{<:Any,LowerBoundary}, v::AbstractVector)
     apply_stencil(op.stencil,v,1)
 end
 
-function LazyTensors.apply(op::BoundaryOperator{<:Any,Upper}, v::AbstractVector)
+function LazyTensors.apply(op::BoundaryOperator{<:Any,UpperBoundary}, v::AbstractVector)
     apply_stencil_backwards(op.stencil,v,op.size)
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Lower}, v::AbstractArray{<:Any,0}, i::Index{Lower})
+function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,LowerBoundary}, v::AbstractArray{<:Any,0}, i::Index{Lower})
     return op.stencil[Int(i)-1]*v[]
 end
 
-function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Upper}, v::AbstractArray{<:Any,0}, i::Index{Upper})
+function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,UpperBoundary}, v::AbstractArray{<:Any,0}, i::Index{Upper})
     return op.stencil[op.size[1] - Int(i)]*v[]
 end
 

--- a/test/Grids/equidistant_grid_test.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -1,6 +1,5 @@
 using Sbplib.Grids
 using Test
-using Sbplib.RegionIndices
 using Sbplib.LazyTensors
 
 
@@ -63,24 +62,24 @@
 
     @testset "boundary_identifiers" begin
         g = EquidistantGrid(0:0.1:10)
-        @test boundary_identifiers(g) == (Lower(), Upper())
+        @test boundary_identifiers(g) == (LowerBoundary(), UpperBoundary())
         @inferred boundary_identifiers(g)
     end
 
     @testset "boundary_grid" begin
         g = EquidistantGrid(0:0.1:1)
-        @test boundary_grid(g, Lower()) == ZeroDimGrid(0.0)
-        @test boundary_grid(g, Upper()) == ZeroDimGrid(1.0)
+        @test boundary_grid(g, LowerBoundary()) == ZeroDimGrid(0.0)
+        @test boundary_grid(g, UpperBoundary()) == ZeroDimGrid(1.0)
     end
 
     @testset "boundary_indices" begin
         g = EquidistantGrid(0:0.1:1)
-        @test boundary_indices(g, Lower()) == (1,)
-        @test boundary_indices(g, Upper()) == (11,)
+        @test boundary_indices(g, LowerBoundary()) == (1,)
+        @test boundary_indices(g, UpperBoundary()) == (11,)
 
         g = EquidistantGrid(2:0.1:10)
-        @test boundary_indices(g, Lower()) == (1,)
-        @test boundary_indices(g, Upper()) == (81,)
+        @test boundary_indices(g, LowerBoundary()) == (1,)
+        @test boundary_indices(g, UpperBoundary()) == (81,)
 
     end
 

--- a/test/Grids/tensor_grid_test.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/test/Grids/tensor_grid_test.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -1,7 +1,6 @@
 using Test
 using Sbplib.Grids
 using StaticArrays
-using Sbplib.RegionIndices
 
 @testset "TensorGrid" begin
     g₁ = EquidistantGrid(range(0,1,length=11))
@@ -170,13 +169,13 @@
     end
 
     @testset "boundary_identifiers" begin
-        @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (Lower,Upper,Lower,Upper))
-        @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,Lower}(),TensorGridBoundary{1,Upper}())
+        @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (LowerBoundary,UpperBoundary,LowerBoundary,UpperBoundary))
+        @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,LowerBoundary}(),TensorGridBoundary{1,UpperBoundary}())
     end
 
     @testset "boundary_grid" begin
-        @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂)
-        @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄)
+        @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, UpperBoundary}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂)
+        @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, UpperBoundary}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄)
     end
 
     @testset "boundary_indices" begin
@@ -184,14 +183,14 @@
         g₂ = EquidistantGrid(range(2,3,length=6))
         g₄ = ZeroDimGrid(@SVector[1,2])
 
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Lower}()) == (1,:)
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == (11,:)
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Lower}()) == (:,1)
-        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Upper}()) == (:,6)
-        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Lower}()) == (1,)
-        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == (11,)
-        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Lower}()) == (1,)
-        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Upper}()) == (11,)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, LowerBoundary}()) == (1,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, UpperBoundary}()) == (11,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, LowerBoundary}()) == (:,1)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, UpperBoundary}()) == (:,6)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, LowerBoundary}()) == (1,)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, UpperBoundary}()) == (11,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, LowerBoundary}()) == (1,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, UpperBoundary}()) == (11,)
     end
 end
 

--- a/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -15,8 +15,8 @@
     g = 3.14
     f(x,y,z) = x^2+y^2+z^2
     @testset "Constructors" begin
-        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower}
-        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,Lower}}
+        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,LowerBoundary}
+        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,LowerBoundary}}
     end
 
     @testset "boundary" begin

--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -13,12 +13,12 @@
     g_1D = EquidistantGrid(range(0,1,length=11))
 
     @testset "Constructors" begin
-        @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T
-        @test BoundaryOperator(g_1D, closure_stencil, Upper()) isa LazyTensor{T,0,1} where T
+        @test BoundaryOperator(g_1D, closure_stencil, LowerBoundary()) isa LazyTensor{T,0,1} where T
+        @test BoundaryOperator(g_1D, closure_stencil, UpperBoundary()) isa LazyTensor{T,0,1} where T
     end
 
-    op_l = BoundaryOperator(g_1D, closure_stencil, Lower())
-    op_r = BoundaryOperator(g_1D, closure_stencil, Upper())
+    op_l = BoundaryOperator(g_1D, closure_stencil, LowerBoundary())
+    op_r = BoundaryOperator(g_1D, closure_stencil, UpperBoundary())
 
     @testset "Sizes" begin
         @test domain_size(op_l) == (11,)

--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -14,19 +14,19 @@
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
-            e_l = boundary_restriction(g_1D,stencil_set,Lower())
-            @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Lower())
-            @test e_l isa BoundaryOperator{T,Lower} where T
+            e_l = boundary_restriction(g_1D,stencil_set,LowerBoundary())
+            @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),LowerBoundary())
+            @test e_l isa BoundaryOperator{T,LowerBoundary} where T
             @test e_l isa LazyTensor{T,0,1} where T
 
-            e_r = boundary_restriction(g_1D,stencil_set,Upper())
-            @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Upper())
-            @test e_r isa BoundaryOperator{T,Upper} where T
+            e_r = boundary_restriction(g_1D,stencil_set,UpperBoundary())
+            @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),UpperBoundary())
+            @test e_r isa BoundaryOperator{T,UpperBoundary} where T
             @test e_r isa LazyTensor{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}())
+            e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,UpperBoundary}())
             @test e_w isa InflatedTensor
             @test e_w isa LazyTensor{T,1,2} where T
         end

--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sat Jun 29 17:06:27 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sun Jul 07 15:15:12 2024 -0700
@@ -12,19 +12,19 @@
     @testset "normal_derivative" begin
     	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            d_l = normal_derivative(g_1D, stencil_set, Lower())
-            @test d_l == normal_derivative(g_1D, stencil_set, Lower())
-            @test d_l isa BoundaryOperator{T,Lower} where T
+            d_l = normal_derivative(g_1D, stencil_set, LowerBoundary())
+            @test d_l == normal_derivative(g_1D, stencil_set, LowerBoundary())
+            @test d_l isa BoundaryOperator{T,LowerBoundary} where T
             @test d_l isa LazyTensor{T,0,1} where T
         end
         @testset "2D" begin
-            d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
-            d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}())
+            d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}())
+            d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,UpperBoundary}())
             Ix = IdentityTensor{Float64}((size(g_2D)[1],))
             Iy = IdentityTensor{Float64}((size(g_2D)[2],))
-            d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower())
-            d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper())
-            @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}())
+            d_l = normal_derivative(g_2D.grids[1], stencil_set, LowerBoundary())
+            d_r = normal_derivative(g_2D.grids[2], stencil_set, UpperBoundary())
+            @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}())
             @test d_w ==  d_l⊗Iy
             @test d_n ==  Ix⊗d_r
             @test d_w isa LazyTensor{T,1,2} where T