--- a/src/SbpOperators/SbpOperators.jl	Tue Mar 22 11:27:45 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Tue Mar 22 12:33:58 2022 +0100
@@ -9,7 +9,7 @@
 export normal_derivative
 export first_derivative
 export second_derivative
-export dissipation
+export undevided_dissipation
 
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors

--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl	Tue Mar 22 11:27:45 2022 +0100
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Tue Mar 22 12:33:58 2022 +0100
@@ -1,11 +1,26 @@
-function dissipation(g::EquidistantGrid, p, direction)
-    h_inv = inverse_spacing(g)[direction]
+function undevided_dissipation(g::EquidistantGrid, p, direction)
+    T = eltype(g)
+    interior_weights = T.(dissipation_interior_weights(p))
 
-    # D = volume_operator(g,CenteredStencil(1),(CenteredStencil(1)), )
-    return nothing, nothing
+    D  = stencil_operator_distinct_closures(
+        g,
+        dissipation_interior_stencil(interior_weights),
+        dissipation_lower_closure_stencils(interior_weights),
+        dissipation_upper_closure_stencils(interior_weights),
+        direction,
+    )
+    Dᵀ = stencil_operator_distinct_closures(
+        g,
+        dissipation_transpose_interior_stencil(interior_weights),
+        dissipation_transpose_lower_closure_stencils(interior_weights),
+        dissipation_transpose_upper_closure_stencils(interior_weights),
+        direction,
+    )
+
+    return D, Dᵀ
 end
 
-dissipation(g::EquidistantGrid{1}, p) = dissipation(g, p, 1)
+undevided_dissipation(g::EquidistantGrid{1}, p) = undevided_dissipation(g, p, 1)
 
 function dissipation_interior_weights(p)
    if p == 0

--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Mar 22 11:27:45 2022 +0100
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Mar 22 12:33:58 2022 +0100
@@ -26,27 +26,36 @@
     x^k/factorial(k)
 end
 
-@testset "dissipation" begin
+@testset "undevided_dissipation" begin
     g = EquidistantGrid(20, 0., 11.)
-    D,Dᵀ = dissipation(g, 1)
+    D,Dᵀ = undevided_dissipation(g, 1)
 
-    @test_broken D isa LazyTensor{Float64,1,1} where T
-    @test_broken Dᵀ isa LazyTensor{Float64,1,1} where T
+    @test D isa LazyTensor{Float64,1,1} where T
+    @test Dᵀ isa LazyTensor{Float64,1,1} where T
 
      @testset "Accuracy conditions" begin
         N = 20
         g = EquidistantGrid(N, 0//1,2//1)
+        h = spacing(g)[1]
         @testset "D_$p" for p ∈ [1,2,3,4]
-            D,Dᵀ = dissipation(g, p)
+            D,Dᵀ = undevided_dissipation(g, p)
 
-            @testset "x^$k" for k ∈ 0:1
-                v = evalOn(g, x->monomial(x,k))
+            @testset "x^$k" for k ∈ 0:p
+                v  = evalOn(g, x->monomial(x,k))
+                vₚₓ = evalOn(g, x->monomial(x,k-p))
 
-                x, = points(g)[10]
-                @test_broken (D*v)[10] == monomial(x,k-1)
+                @test D*v == h^p * vₚₓ
             end
 
-            # Test Dᵀ works backwards and interior forwards
+
+            @testset "x^$k" for k ∈ 0:p
+                v  = evalOn(g, x->monomial(x,k))
+                vₚₓ = evalOn(g, x->monomial(x,k-p))
+
+                for i ∈ SbpOperators.lower_closure_size(Dᵀ)+1:N-SbpOperators.upper_closure_size(Dᵀ)
+                    @test (Dᵀ*v)[i] == h^p * vₚₓ[i]
+                end
+            end
         end
     end
 end