--- a/src/LazyTensors/lazy_tensor_operations.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -270,6 +270,20 @@
 LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms)
 
 
+
+"""
+    inflate(tm, sz, dir)
+
+Inflate `tm` with identity tensors in all directions `d` for `d != dir`.
+
+# TODO: Describe when it is useful
+"""
+function inflate(tm::LazyTensor, sz, dir)
+    Is = IdentityTensor{eltype(tm)}.(sz)
+    parts = Base.setindex(Is, tm, dir)
+    return foldl(⊗, parts)
+end
+
 function check_domain_size(tm::LazyTensor, sz)
     if domain_size(tm) != sz
         throw(DomainSizeMismatch(tm,sz))

--- a/src/LazyTensors/tuple_manipulation.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/src/LazyTensors/tuple_manipulation.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -74,3 +74,23 @@
 flatten_tuple(t::NTuple{N, Number} where N) = t
 flatten_tuple(t::Tuple) = ((flatten_tuple.(t)...)...,) # simplify?
 flatten_tuple(ts::Vararg) = flatten_tuple(ts)
+
+
+function left_pad_tuple(t, val, N)
+    if N < length(t)
+        throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements"))
+    end
+
+    padding = ntuple(i->val, N-length(t))
+    return (padding..., t...)
+end
+
+function right_pad_tuple(t, val, N)
+    if N < length(t)
+        throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements"))
+    end
+
+    padding = ntuple(i->val, N-length(t))
+    return (t..., padding...)
+end
+

--- a/src/SbpOperators/SbpOperators.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -20,6 +20,7 @@
 export normal_derivative
 export first_derivative
 export second_derivative
+export undivided_dissipation
 
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
@@ -35,10 +36,12 @@
 include("stencil.jl")
 include("stencil_set.jl")
 include("volumeops/volume_operator.jl")
+include("volumeops/stencil_operator_distinct_closures.jl")
 include("volumeops/constant_interior_scaling_operator.jl")
 include("volumeops/derivatives/first_derivative.jl")
 include("volumeops/derivatives/second_derivative.jl")
 include("volumeops/derivatives/second_derivative_variable.jl")
+include("volumeops/derivatives/dissipation.jl")
 include("volumeops/laplace/laplace.jl")
 include("volumeops/inner_products/inner_product.jl")
 include("volumeops/inner_products/inverse_inner_product.jl")

--- a/src/SbpOperators/boundaryops/boundary_operator.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -12,19 +12,14 @@
 function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary)
     #TODO:Check that dim(boundary) <= Dim?
 
-    # Create 1D boundary operator
-    r = region(boundary)
     d = dim(boundary)
-    op = BoundaryOperator(restrict(grid, d), closure_stencil, r)
+    op = BoundaryOperator(restrict(grid, d), closure_stencil, region(boundary))
 
     # Create 1D IdentityTensors for each coordinate direction
     one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
     Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids))
 
-    # Formulate the correct outer product sequence of the identity mappings and
-    # the boundary operator
-    parts = Base.setindex(Is, op, d)
-    return foldl(⊗, parts)
+    return LazyTensors.inflate(op, size(grid), d)
 end
 
 """

--- a/src/SbpOperators/stencil.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/src/SbpOperators/stencil.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -85,6 +85,21 @@
     return w
 end
 
+function left_pad(s::Stencil, N)
+    weights = LazyTensors.left_pad_tuple(s.weights, zero(eltype(s)), N)
+    range = (first(s.range) - (N - length(s.weights))):last(s.range)
+
+    return Stencil(range, weights)
+end
+
+function right_pad(s::Stencil, N)
+    weights = LazyTensors.right_pad_tuple(s.weights, zero(eltype(s)), N)
+    range = first(s.range):(last(s.range) + (N - length(s.weights)))
+
+    return Stencil(range, weights)
+end
+
+
 
 struct NestedStencil{T,N,M}
     s::Stencil{Stencil{T,N},M}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -0,0 +1,87 @@
+function undivided_dissipation(g::EquidistantGrid, p, direction)
+    T = eltype(g)
+    interior_weights = T.(dissipation_interior_weights(p))
+
+    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
+
+undivided_dissipation(g::EquidistantGrid{1}, p) = undivided_dissipation(g, p, 1)
+
+function dissipation_interior_weights(p)
+   if p == 0
+       return (1,)
+   end
+
+   return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0)
+end
+
+midpoint(weights) = length(weights)÷2 + 1
+midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
+
+function dissipation_interior_stencil(weights)
+    return Stencil(weights..., center=midpoint(weights))
+end
+function dissipation_transpose_interior_stencil(weights)
+    if iseven(length(weights))
+        weights = map(-, weights)
+    end
+
+    return Stencil(weights..., center=midpoint_transpose(weights))
+end
+
+dissipation_lower_closure_size(weights) = midpoint(weights) - 1
+dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
+
+dissipation_lower_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=i                       ), dissipation_lower_closure_size(interior_weights))
+dissipation_upper_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=length(interior_weights)-dissipation_upper_closure_size(interior_weights)+i), dissipation_upper_closure_size(interior_weights))
+
+function dissipation_transpose_lower_closure_stencils(interior_weights)
+    closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
+
+    N = maximum(s->length(s.weights), closure)
+    return right_pad.(closure, N)
+end
+
+function dissipation_transpose_upper_closure_stencils(interior_weights)
+    closure = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights)))
+
+    N = maximum(s->length(s.weights), closure)
+    return left_pad.(closure, N)
+end
+
+
+function dissipation_transpose_lower_closure_stencil(interior_weights, i)
+    w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights))
+
+    for k ∈ i:-1:1
+        w = (w..., interior_weights[k])
+    end
+
+    return Stencil(w..., center = i)
+end
+
+function dissipation_transpose_upper_closure_stencil(interior_weights, i)
+    j = length(interior_weights)+1-i
+    w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights))
+
+    for k ∈ j:1:length(interior_weights)
+        w = (interior_weights[k], w...)
+    end
+
+    return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1)
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -0,0 +1,58 @@
+function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction)
+    op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure)
+    return LazyTensors.inflate(op, size(grid), direction)
+end
+
+"""
+    StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1}
+
+A one dimensional stencil operator with separate closures for the two boundaries.
+
+`StencilOperatorDistinctClosures` can be contrasted to `VolumeOperator` in
+that it has different closure stencils for the upper and lower boundary.
+`VolumeOperator` uses the same closure for both boundaries. Having distinct
+closures is useful for representing operators with skewed stencils like upwind
+operators.
+
+See also: [`VolumeOperator`](@ref)
+"""
+struct StencilOperatorDistinctClosures{T,K,N,M,LC<:NTuple{N,Stencil{T,L}} where L, UC<:NTuple{M,Stencil{T,L}} where L} <: LazyTensor{T,1,1}
+    inner_stencil::Stencil{T,K}
+    lower_closure::LC
+    upper_closure::UC
+    size::Tuple{Int}
+end
+
+function StencilOperatorDistinctClosures(grid::EquidistantGrid{1}, inner_stencil, lower_closure, upper_closure)
+    return StencilOperatorDistinctClosures(inner_stencil, Tuple(lower_closure), Tuple(upper_closure), size(grid))
+end
+
+lower_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = N
+upper_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = M
+
+LazyTensors.range_size(op::StencilOperatorDistinctClosures) = op.size
+LazyTensors.domain_size(op::StencilOperatorDistinctClosures) = op.size
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Lower})
+    return @inbounds apply_stencil(op.lower_closure[Int(i)], v, Int(i))
+end
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Interior})
+    return apply_stencil(op.inner_stencil, v, Int(i))
+end
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Upper})
+    stencil_index = Int(i) - (op.size[1]-upper_closure_size(op))
+    return @inbounds apply_stencil(op.upper_closure[stencil_index], v, Int(i))
+end
+
+function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i)
+    if i <= lower_closure_size(op)
+        LazyTensors.apply(op, v, Index(i, Lower))
+    elseif i > op.size[1]-upper_closure_size(op)
+        LazyTensors.apply(op, v, Index(i, Upper))
+    else
+        LazyTensors.apply(op, v, Index(i, Interior))
+    end
+end
+# TODO: Move this to LazyTensors when we have the region communication down.

--- a/src/SbpOperators/volumeops/volume_operator.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -12,19 +12,13 @@
 function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction)
     #TODO: Check that direction <= Dim?
 
-    # Create 1D volume operator in along coordinate direction
     op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity)
-    # Create 1D IdentityTensors for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
-    Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids))
-    # Formulate the correct outer product sequence of the identity mappings and
-    # the volume operator
-    parts = Base.setindex(Is, op, direction)
-    return foldl(⊗, parts)
+    return LazyTensors.inflate(op, size(grid), direction)
 end
+# TBD: Should the inflation happen here or should we remove this method and do it at the caller instead?
 
 """
-    VolumeOperator{T,N,M,K} <: TensorOperator{T,1}
+    VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
 Implements a one-dimensional constant coefficients volume operator
 """
 struct VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1}
@@ -58,3 +52,4 @@
 function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i)
     return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i)
 end
+# TODO: Move this to LazyTensors when we have the region communication down.

--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -358,3 +358,17 @@
         @test I1⊗Ã⊗I2 == InflatedTensor(I1, Ã, I2)
     end
 end
+
+@testset "inflate" begin
+    I = LazyTensors.inflate(IdentityTensor(),(3,4,5,6), 2)
+    @test I isa LazyTensor{Float64, 3,3}
+    @test range_size(I) == (3,5,6)
+    @test domain_size(I) == (3,5,6)
+
+    @test LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 1) == InflatedTensor(IdentityTensor{Float64}(),ScalingTensor(1., (4,)),IdentityTensor(4,5,6))
+    @test LazyTensors.inflate(ScalingTensor(2., (1,)),(3,4,5,6), 2) == InflatedTensor(IdentityTensor(3),ScalingTensor(2., (1,)),IdentityTensor(5,6))
+    @test LazyTensors.inflate(ScalingTensor(3., (6,)),(3,4,5,6), 4) == InflatedTensor(IdentityTensor(3,4,5),ScalingTensor(3., (6,)),IdentityTensor{Float64}())
+
+    @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 0)
+    @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 5)
+end

--- a/test/LazyTensors/tuple_manipulation_test.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/test/LazyTensors/tuple_manipulation_test.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -60,3 +60,19 @@
     @test LazyTensors.flatten_tuple((1,2,(3,(4,5)),6)) == (1,2,3,4,5,6)
     @test LazyTensors.flatten_tuple(((1,2),(3,4),(5,),6)) == (1,2,3,4,5,6)
 end
+
+@testset "left_pad_tuple" begin
+    @test LazyTensors.left_pad_tuple((1,2), 0, 2) == (1,2)
+    @test LazyTensors.left_pad_tuple((1,2), 0, 3) == (0,1,2)
+    @test LazyTensors.left_pad_tuple((3,2), 1, 6) == (1,1,1,1,3,2)
+
+    @test_throws DomainError(0, "Can't pad tuple of length 2 to 0 elements") LazyTensors.left_pad_tuple((1,2), 0, 0) == (1,2)
+end
+
+@testset "right_pad_tuple" begin
+    @test LazyTensors.right_pad_tuple((1,2), 0, 2) == (1,2)
+    @test LazyTensors.right_pad_tuple((1,2), 0, 3) == (1,2,0)
+    @test LazyTensors.right_pad_tuple((3,2), 1, 6) == (3,2,1,1,1,1)
+
+    @test_throws DomainError(0, "Can't pad tuple of length 2 to 0 elements") LazyTensors.right_pad_tuple((1,2), 0, 0) == (1,2)
+end

--- a/test/SbpOperators/stencil_test.jl	Wed Mar 23 13:32:51 2022 +0100
+++ b/test/SbpOperators/stencil_test.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -62,6 +62,23 @@
     end
 end
 
+@testset "left_pad" begin
+    @test SbpOperators.left_pad(Stencil(1,1, center = 1), 2) == Stencil(1,1, center=1)
+    @test SbpOperators.left_pad(Stencil(1,1, center = 1), 3) == Stencil(0,1,1, center=2)
+    @test SbpOperators.left_pad(Stencil(2,3, center = 2), 4) == Stencil(0,0,2,3, center=4)
+
+    @test SbpOperators.left_pad(Stencil(2.,3., center = 2), 4) == Stencil(0.,0.,2.,3., center=4)
+end
+
+@testset "right_pad" begin
+    @test SbpOperators.right_pad(Stencil(1,1, center = 1), 2) == Stencil(1,1, center=1)
+    @test SbpOperators.right_pad(Stencil(1,1, center = 1), 3) == Stencil(1,1,0, center=1)
+    @test SbpOperators.right_pad(Stencil(2,3, center = 2), 4) == Stencil(2,3,0,0, center=2)
+
+    @test SbpOperators.right_pad(Stencil(2.,3., center = 2), 4) == Stencil(2.,3.,0.,0., center=2)
+end
+
+
 @testset "NestedStencil" begin
 
     @testset "Constructors" begin
@@ -170,5 +187,4 @@
         @inferred SbpOperators.apply_stencil_backwards(s_int,   c_float, v_float, 2)
         @inferred SbpOperators.apply_stencil_backwards(s_float, c_float, v_int,   2)
     end
-
 end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -0,0 +1,217 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+using Sbplib.LazyTensors
+
+using Sbplib.SbpOperators: Stencil
+
+using Sbplib.SbpOperators: dissipation_interior_weights
+using Sbplib.SbpOperators: dissipation_interior_stencil, dissipation_transpose_interior_stencil
+using Sbplib.SbpOperators: midpoint, midpoint_transpose
+using Sbplib.SbpOperators: dissipation_lower_closure_size, dissipation_upper_closure_size
+using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
+using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
+
+"""
+    monomial(x,k)
+
+Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
+Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
+"""
+function monomial(x,k)
+    if k < 0
+        return zero(x)
+    end
+    x^k/factorial(k)
+end
+
+@testset "undivided_dissipation" begin
+    g = EquidistantGrid(20, 0., 11.)
+    D,Dᵀ = undivided_dissipation(g, 1)
+
+    @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 = only(spacing(g))
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = undivided_dissipation(g, p)
+
+            @testset "x^$k" for k ∈ 0:p
+                v  = evalOn(g, x->monomial(x,k))
+                vₚₓ = evalOn(g, x->monomial(x,k-p))
+
+                @test D*v == h^p * vₚₓ
+            end
+        end
+    end
+
+    @testset "tanspose equality" begin
+        function get_matrix(D)
+            N = only(range_size(D))
+            M = only(domain_size(D))
+
+            Dmat = zeros(N,M)
+            e = zeros(M)
+            for i ∈ 1:M
+                if i > 1
+                    e[i-1] = 0.
+                end
+                e[i] = 1.
+                Dmat[:,i] = D*e
+            end
+
+            return Dmat
+        end
+
+        g = EquidistantGrid(11, 0., 1.)
+        @testset "D_$p" for p ∈ [1,2,3,4]
+            D,Dᵀ = undivided_dissipation(g, p)
+
+            D̄  = get_matrix(D)
+            D̄ᵀ = get_matrix(Dᵀ)
+
+            @test D̄ == D̄ᵀ'
+        end
+    end
+end
+
+@testset "dissipation_interior_weights" begin
+    @test dissipation_interior_weights(1) == (-1, 1)
+    @test dissipation_interior_weights(2) == (1,-2, 1)
+    @test dissipation_interior_weights(3) == (-1, 3,-3, 1)
+    @test dissipation_interior_weights(4) == (1, -4, 6, -4, 1)
+end
+
+@testset "dissipation_interior_stencil" begin
+    @test dissipation_interior_stencil(dissipation_interior_weights(1)) == Stencil(-1, 1, center=2)
+    @test dissipation_interior_stencil(dissipation_interior_weights(2)) == Stencil( 1,-2, 1, center=2)
+    @test dissipation_interior_stencil(dissipation_interior_weights(3)) == Stencil(-1, 3,-3, 1, center=3)
+    @test dissipation_interior_stencil(dissipation_interior_weights(4)) == Stencil( 1,-4, 6,-4, 1, center=3)
+end
+
+@testset "dissipation_transpose_interior_stencil" begin
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(1)) == Stencil(1,-1, center=1)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(2)) == Stencil(1,-2, 1, center=2)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(3)) == Stencil(1,-3, 3,-1, center=2)
+    @test dissipation_transpose_interior_stencil(dissipation_interior_weights(4)) == Stencil(1,-4, 6,-4, 1, center=3)
+end
+
+@testset "midpoint" begin
+    @test midpoint((1,1)) == 2
+    @test midpoint((1,1,1)) == 2
+    @test midpoint((1,1,1,1)) == 3
+    @test midpoint((1,1,1,1,1)) == 3
+end
+
+@testset "midpoint_transpose" begin
+    @test midpoint_transpose((1,1)) == 1
+    @test midpoint_transpose((1,1,1)) == 2
+    @test midpoint_transpose((1,1,1,1)) == 2
+    @test midpoint_transpose((1,1,1,1,1)) == 3
+end
+
+@testset "dissipation_lower_closure_size" begin
+    @test dissipation_lower_closure_size((1,1)) == 1
+    @test dissipation_lower_closure_size((1,1,1)) == 1
+    @test dissipation_lower_closure_size((1,1,1,1)) == 2
+    @test dissipation_lower_closure_size((1,1,1,1,1)) == 2
+end
+
+@testset "dissipation_upper_closure_size" begin
+    @test dissipation_upper_closure_size((1,1)) == 0
+    @test dissipation_upper_closure_size((1,1,1)) == 1
+    @test dissipation_upper_closure_size((1,1,1,1)) == 1
+    @test dissipation_upper_closure_size((1,1,1,1,1)) == 2
+end
+
+@testset "dissipation_lower_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil(-1, 1, center=1),
+        ),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, center=1),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,3,-3,1, center=1),
+            Stencil(-1,3,-3,1, center=2),
+        ),
+        (1, -4, 6, -4, 1) => (
+            Stencil(1, -4, 6, -4, 1, center=1),
+            Stencil(1, -4, 6, -4, 1, center=2),
+        )
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_lower_closure_stencils(w) == closure_stencils
+    end
+end
+
+@testset "dissipation_upper_closure_stencils" begin
+    cases = (
+        (-1,1) => (),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, center=3),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,3,-3,1, center=4),
+        ),
+        (1, -4, 6, -4, 1) => (
+            Stencil(1, -4, 6, -4, 1, center=4),
+            Stencil(1, -4, 6, -4, 1, center=5),
+        )
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_upper_closure_stencils(w) == closure_stencils
+    end
+end
+
+
+@testset "dissipation_transpose_lower_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil(-1,-1, 0, center=1),
+            Stencil( 1, 1,-1, center=2),
+        ),
+        (1,-2,1) => (
+            Stencil( 1, 1, 0, 0, center=1),
+            Stencil(-2,-2, 1, 0, center=2),
+            Stencil( 1, 1,-2, 1, center=3),
+        ),
+        (-1,3,-3,1) => (
+            Stencil(-1,-1,-1, 0, 0, 0, center=1),
+            Stencil( 3, 3, 3,-1, 0, 0, center=2),
+            Stencil(-3,-3,-3, 3,-1, 0, center=3),
+            Stencil( 1, 1, 1,-3, 3,-1, center=4),
+        ),
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_transpose_lower_closure_stencils(w) == closure_stencils
+    end
+end
+
+@testset "dissipation_transpose_upper_closure_stencils" begin
+    cases = (
+        (-1,1) => (
+            Stencil( 1,-1, center = 1),
+            Stencil( 0, 1, center = 2),
+        ),
+        (1,-2,1) => (
+            Stencil( 1,-2, 1, 1, center=2),
+            Stencil( 0, 1,-2,-2, center=3),
+            Stencil( 0, 0, 1, 1, center=4),
+        ),
+        (-1,3,-3,1) => (
+            Stencil( 1,-3, 3,-1,-1, center=2),
+            Stencil( 0, 1,-3, 3, 3, center=3),
+            Stencil( 0, 0, 1,-3,-3, center=4),
+            Stencil( 0, 0, 0, 1, 1, center=5),
+        ),
+    )
+    @testset "interior_weights = $w" for (w, closure_stencils) ∈ cases
+        @test dissipation_transpose_upper_closure_stencils(w) == closure_stencils
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Thu Mar 24 12:35:14 2022 +0100
@@ -0,0 +1,83 @@
+using Test
+
+using Sbplib.SbpOperators
+using Sbplib.Grids
+# using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+
+import Sbplib.SbpOperators.Stencil
+import Sbplib.SbpOperators.StencilOperatorDistinctClosures
+import Sbplib.SbpOperators.stencil_operator_distinct_closures
+
+@testset "stencil_operator_distinct_closures" begin
+    lower_closure = (
+        Stencil(-1,1, center=1),
+    )
+
+    inner_stencil = Stencil(-2,2, center=1)
+
+    upper_closure = (
+        Stencil(-3,3, center=1),
+        Stencil(-4,4, center=2),
+    )
+
+    g₁ = EquidistantGrid(5, 0., 1.)
+    g₂ = EquidistantGrid((5,5), (0.,0.), (1.,1.))
+    h = 1/4
+
+    A₁  = stencil_operator_distinct_closures(g₁, inner_stencil, lower_closure, upper_closure, 1)
+    A₂¹ = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 1)
+    A₂² = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 2)
+
+    v₁ = evalOn(g₁, x->x)
+
+    u = [1., 2., 2., 3., 4.]*h
+    @test A₁*v₁ == u
+
+    v₂ = evalOn(g₂, (x,y)-> x + 3y)
+    @test A₂¹*v₂ == repeat(u, 1, 5)
+    @test A₂²*v₂ == repeat(3u', 5, 1)
+end
+
+@testset "StencilOperatorDistinctClosures" begin
+    g = EquidistantGrid(11, 0., 1.)
+
+    lower_closure = (
+        Stencil(-1,1, center=1),
+        Stencil(-2,2, center=2),
+    )
+
+    inner_stencil = Stencil(-3,3, center=1)
+
+    upper_closure = (
+        Stencil(4,-4,4, center=1),
+        Stencil(5,-5,5, center=2),
+        Stencil(6,-6,6, center=3),
+    )
+
+    A = StencilOperatorDistinctClosures(g, inner_stencil, lower_closure, upper_closure)
+    @test A isa LazyTensor{T,1,1} where T
+
+    @test SbpOperators.lower_closure_size(A) == 2
+    @test SbpOperators.upper_closure_size(A) == 3
+
+    @test domain_size(A) == (11,)
+    @test range_size(A) == (11,)
+
+    v = rand(11)
+    @testset "apply" begin
+        # Lower closure
+        @test LazyTensors.apply(A, v, 1) ≈ 1*(-v[1] + v[2])
+        @test LazyTensors.apply(A, v, 2) ≈ 2*(-v[1] + v[2])
+
+        # Interior
+        @test LazyTensors.apply(A, v, 3) ≈ 3*(-v[3] + v[4])
+        @test LazyTensors.apply(A, v, 4) ≈ 3*(-v[4] + v[5])
+        @test LazyTensors.apply(A, v, 8) ≈ 3*(-v[8] + v[9])
+
+        # Upper closure
+        @test LazyTensors.apply(A, v,  9) ≈ 4*(v[9] - v[10] + v[11])
+        @test LazyTensors.apply(A, v, 10) ≈ 5*(v[9] - v[10] + v[11])
+        @test LazyTensors.apply(A, v, 11) ≈ 6*(v[9] - v[10] + v[11])
+    end
+end