--- a/src/Grids/equidistant_grid.jl	Mon May 15 22:37:33 2023 +0200
+++ b/src/Grids/equidistant_grid.jl	Mon May 15 22:37:55 2023 +0200
@@ -14,10 +14,9 @@
     points::R
 end
 
-Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
-
 # Indexing interface
 Base.getindex(g::EquidistantGrid, i) = g.points[i]
+Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
 Base.firstindex(g::EquidistantGrid) = firstindex(g.points)
 Base.lastindex(g::EquidistantGrid) = lastindex(g.points)
 

--- a/src/Grids/grid.jl	Mon May 15 22:37:33 2023 +0200
+++ b/src/Grids/grid.jl	Mon May 15 22:37:55 2023 +0200
@@ -71,10 +71,6 @@
 function boundary_grid end
 # TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff?
 
-
-# TODO: Make sure that all grids implement all of the above.
-
-
 """
     eval_on(g::Grid, f)
 

--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl	Mon May 15 22:37:33 2023 +0200
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Mon May 15 22:37:55 2023 +0200
@@ -14,8 +14,11 @@
 function undivided_skewed04 end
 
 function undivided_skewed04(g::TensorGrid, p, direction)
-    op = undivided_skewed04(g.grids[direction], p)
-    return LazyTensors.inflate(op, size(g), direction)
+    D,Dᵀ = undivided_skewed04(g.grids[direction], p)
+    return (
+        LazyTensors.inflate(D, size(g), direction),
+        LazyTensors.inflate(Dᵀ, size(g), direction),
+    )
 end
 
 function undivided_skewed04(g::EquidistantGrid, p)

--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Mon May 15 22:37:33 2023 +0200
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Mon May 15 22:37:55 2023 +0200
@@ -47,8 +47,6 @@
                 @test D*v == h^p * vₚₓ
             end
         end
-
-        # TODO: Add 2D tests
     end
 
     @testset "transpose equality" begin
@@ -79,6 +77,19 @@
             @test D̄ == D̄ᵀ'
         end
     end
+
+    @testset "2D" begin
+        N = 20
+        g = equidistant_grid((N,2N), (0,0), (2,1))
+        h = spacing.(g.grids)
+
+        D,Dᵀ = undivided_skewed04(g, 3, 2)
+
+        v = eval_on(g, x->monomial(x[1],4)*monomial(x[2],3))
+        d³vdy³ = eval_on(g, x->monomial(x[1],4)*monomial(x[2],0))
+
+        @test D*v ≈ h[2]^3*d³vdy³
+    end
 end
 
 @testset "dissipation_interior_weights" begin