--- a/src/SbpOperators/SbpOperators.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -8,12 +8,14 @@
 include("constantstenciloperator.jl")
 include("d2.jl")
 include("readoperator.jl")
+include("volumeops/volume_operator.jl")
 include("laplace/secondderivative.jl")
 include("laplace/laplace.jl")
 include("quadrature/diagonal_inner_product.jl")
 include("quadrature/quadrature.jl")
 include("quadrature/inverse_diagonal_inner_product.jl")
 include("quadrature/inverse_quadrature.jl")
+include("boundaryops/boundary_operator.jl")
 include("boundaryops/boundary_restriction.jl")
 
 end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -0,0 +1,89 @@
+"""
+    boundary_operator(grid,closure_stencil,boundary)
+
+Creates a boundary operator on a `Dim`-dimensional grid for the
+specified `boundary`. The action of the operator is determined by `closure_stencil`.
+
+When `Dim=1`, the corresponding `BoundaryOperator` tensor mapping is returned.
+When `Dim>1`, the `BoundaryOperator` `op` is inflated by the outer product
+of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
+the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
+"""
+function boundary_operator(grid::EquidistantGrid{Dim,T}, closure_stencil::Stencil{T}, boundary::CartesianBoundary) where {Dim,T}
+    # Create 1D boundary operator
+    r = region(boundary)
+    d = dim(boundary)
+    op = BoundaryOperator(restrict(grid, d), closure_stencil, r)
+
+    # Create 1D IdentityMappings for each coordinate direction
+    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
+    Is = IdentityMapping{T}.(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)
+end
+export boundary_operator
+
+"""
+    BoundaryOperator{T,R,N} <: TensorMapping{T,0,1}
+
+Implements the boundary operator `op` for 1D as a `TensorMapping`
+
+`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 prolongation of a zero dimensional array to the whole grid using the same closure stencil.
+"""
+struct BoundaryOperator{T,R<:Region,N} <: TensorMapping{T,0,1}
+    stencil::Stencil{T,N}
+    size::Int
+end
+export BoundaryOperator
+
+BoundaryOperator{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryOperator{T,R,N}(stencil, size)
+
+"""
+    BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region)
+
+Constructs the BoundaryOperator with stencil `closure_stencil` for a one-dimensional `grid`, restricting to
+to the boundary specified by `region`.
+"""
+function BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil::Stencil{T,N}, region::Region) where {T,N}
+    return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1])
+end
+
+"""
+    closure_size(::BoundaryOperator)
+The size of the closure stencil.
+"""
+closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N
+
+LazyTensors.range_size(op::BoundaryOperator) = ()
+LazyTensors.domain_size(op::BoundaryOperator) = (op.size,)
+
+function LazyTensors.apply(op::BoundaryOperator{T,Lower}, v::AbstractVector{T}) where T
+    apply_stencil(op.stencil,v,1)
+end
+
+function LazyTensors.apply(op::BoundaryOperator{T,Upper}, v::AbstractVector{T}) where T
+    apply_stencil_backwards(op.stencil,v,op.size)
+end
+
+function LazyTensors.apply_transpose(op::BoundaryOperator{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T
+    return op.stencil[Int(i)-1]*v[]
+end
+
+function LazyTensors.apply_transpose(op::BoundaryOperator{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T
+    return op.stencil[op.size[1] - Int(i)]*v[]
+end
+
+# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods.
+function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i::Index) where T
+    return zero(T)
+end
+
+function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T
+    r = getregion(i, closure_size(op), op.size)
+    apply_transpose(op, v, Index(i,r))
+end

--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -1,81 +1,15 @@
-"""
-    boundary_restriction(grid,closureStencil,boundary)
-
-Creates a boundary restriction operator on a `Dim`-dimensional grid for the
-specified `boundary`.
-
-When `Dim=1`, the corresponding `BoundaryRestriction` tensor mapping is returned.
-When `Dim>1`, the `BoundaryRestriction` `e` is inflated by the outer product
-of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
-the boundary restriction operator in the y-direction direction is `Ix⊗e⊗Iz`.
-"""
-function boundary_restriction(grid::EquidistantGrid{Dim,T}, closureStencil::Stencil{T,M}, boundary::CartesianBoundary) where {Dim,T,M}
-    # Create 1D boundary restriction operator
-    r = region(boundary)
-    d = dim(boundary)
-    e = BoundaryRestriction(restrict(grid, d), closureStencil, r)
-
-    # Create 1D IdentityMappings for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
-    Is = IdentityMapping{T}.(size.(one_d_grids))
-
-    # Formulate the correct outer product sequence of the identity mappings and
-    # the boundary restriction operator
-    parts = Base.setindex(Is, e, d)
-    return foldl(⊗, parts)
-end
-
-export boundary_restriction
-
-"""
-    BoundaryRestriction{T,R,N} <: TensorMapping{T,0,1}
-
-Implements the boundary operator `e` for 1D as a `TensorMapping`
-
-`e` is the restriction of a grid function to the boundary using some `closureStencil`.
-The boundary to restrict to is determined by `R`.
-
-`e'` is the prolongation of a zero dimensional array to the whole grid using the same `closureStencil`.
 """
-struct BoundaryRestriction{T,R<:Region,N} <: TensorMapping{T,0,1}
-    stencil::Stencil{T,N}
-    size::Int
-end
-export BoundaryRestriction
-
-BoundaryRestriction{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryRestriction{T,R,N}(stencil, size)
+    BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
 
-function BoundaryRestriction(grid::EquidistantGrid{1}, closureStencil::Stencil{T,N}, region::Region) where {T,N}
-    return BoundaryRestriction{T,typeof(region),N}(closureStencil,size(grid)[1])
-end
-
-closure_size(::BoundaryRestriction{T,R,N}) where {T,R,N} = N
-
-LazyTensors.range_size(e::BoundaryRestriction) = ()
-LazyTensors.domain_size(e::BoundaryRestriction) = (e.size,)
-
-function LazyTensors.apply(e::BoundaryRestriction{T,Lower}, v::AbstractVector{T}) where T
-    apply_stencil(e.stencil,v,1)
-end
+Creates the boundary restriction operator `e` as a `TensorMapping`
 
-function LazyTensors.apply(e::BoundaryRestriction{T,Upper}, v::AbstractVector{T}) where T
-    apply_stencil_backwards(e.stencil,v,e.size)
-end
-
-function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T
-    return e.stencil[Int(i)-1]*v[]
-end
+`e` is the restriction of a grid function to the boundary specified by `boundary` or `region` using some `closure_stencil`.
+`e'` is the prolongation of a grid function on the boundary to the whole grid using the same `closure_stencil`.
+On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
+a `BoundaryOperator`. Also see the documentation of `boundary_operator(...)` for more details.
+"""
+BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = boundary_operator(grid, closure_stencil, boundary)
+BoundaryRestriction(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = BoundaryRestriction(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
 
-function LazyTensors.apply_transpose(e::BoundaryRestriction{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T
-    return e.stencil[e.size[1] - Int(i)]*v[]
-end
-
-# Catch all combinations of Lower, Upper and Interior not caught by the two previous methods.
-function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i::Index) where T
-    return zero(T)
-end
-
-function LazyTensors.apply_transpose(e::BoundaryRestriction{T}, v::AbstractArray{T,0}, i) where T
-    r = getregion(i, closure_size(e), e.size)
-    apply_transpose(e, v, Index(i,r))
-end
+export BoundaryRestriction

--- a/src/SbpOperators/laplace/laplace.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/laplace/laplace.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -1,49 +1,18 @@
-export Laplace
-"""
-    Laplace{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the Laplace operator `L` in Dim dimensions as a tensor operator
-The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
-tensor operators.
-"""
-#export quadrature, inverse_quadrature, boundary_quadrature, boundary_value, normal_derivative
-struct Laplace{Dim,T,N,M,K} <: TensorMapping{T,Dim,Dim}
-    D2::NTuple{Dim,SecondDerivative{T,N,M,K}}
-end
-
-function Laplace(g::EquidistantGrid{Dim}, innerStencil, closureStencils) where Dim
-    D2 = ()
-    for i ∈ 1:Dim
-        D2 = (D2..., SecondDerivative(restrict(g,i), innerStencil, closureStencils))
+# """
+#     Laplace{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+#
+# Implements the Laplace operator `L` in Dim dimensions as a tensor operator
+# The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
+# tensor operators.
+# """
+function Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
+    Δ = SecondDerivative(grid, inner_stencil, closure_stencils, 1)
+    for d = 2:Dim
+        Δ += SecondDerivative(grid, inner_stencil, closure_stencils, d)
     end
-
-    return Laplace(D2)
+    return Δ
 end
-
-LazyTensors.range_size(L::Laplace) = getindex.(range_size.(L.D2),1)
-LazyTensors.domain_size(L::Laplace) = getindex.(domain_size.(L.D2),1)
-
-function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim}
-    error("not implemented")
-end
-
-# u = L*v
-function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, i) where T
-    @inbounds u = LazyTensors.apply(L.D2[1],v,i)
-    return u
-end
-
-function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, i, j) where T
-    # 2nd x-derivative
-    @inbounds vx = view(v, :, Int(j))
-    @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , i)
-
-    # 2nd y-derivative
-    @inbounds vy = view(v, Int(i), :)
-    @inbounds uᵢ += LazyTensors.apply(L.D2[2], vy , j)
-
-    return uᵢ
-end
+export Laplace
 
 # quadrature(L::Laplace) = Quadrature(L.op, L.grid)
 # inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)

--- a/src/SbpOperators/laplace/secondderivative.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/laplace/secondderivative.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -1,43 +1,6 @@
-"""
-    SecondDerivative{T<:Real,N,M,K} <: TensorOperator{T,1}
-Implements the Laplace tensor operator `L` with constant grid spacing and coefficients
-in 1D dimension
-"""
-
-struct SecondDerivative{T,N,M,K} <: TensorMapping{T,1,1}
-    h_inv::T # The grid spacing could be included in the stencil already. Preferable?
-    innerStencil::Stencil{T,N}
-    closureStencils::NTuple{M,Stencil{T,K}}
-    size::NTuple{1,Int}
-end
-export SecondDerivative
-
-function SecondDerivative(grid::EquidistantGrid{1}, innerStencil, closureStencils)
-    h_inv = inverse_spacing(grid)[1]
-    return SecondDerivative(h_inv, innerStencil, closureStencils, size(grid))
+function SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
+    h_inv = inverse_spacing(grid)[direction]
+    return volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
 end
-
-LazyTensors.range_size(D2::SecondDerivative) = D2.size
-LazyTensors.domain_size(D2::SecondDerivative) = D2.size
-
-# Apply for different regions Lower/Interior/Upper or Unknown region
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Lower}) where T
-    return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.closureStencils[Int(i)], v, Int(i))
-end
-
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Interior}) where T
-    return @inbounds D2.h_inv*D2.h_inv*apply_stencil(D2.innerStencil, v, Int(i))
-end
-
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i::Index{Upper}) where T
-    N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v))
-    return @inbounds D2.h_inv*D2.h_inv*apply_stencil_backwards(D2.closureStencils[N-Int(i)+1], v, Int(i))
-end
-
-function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, i) where T
-    N = length(v)  # TODO: Use domain_size here instead?
-    r = getregion(i, closuresize(D2), N)
-    return LazyTensors.apply(D2, v, Index(i, r))
-end
-
-closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} = M
+SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils) = SecondDerivative(grid,inner_stencil,closure_stencils,1)
+export SecondDerivative

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -0,0 +1,60 @@
+"""
+    volume_operator(grid,inner_stencil,closure_stencils,parity,direction)
+Creates a volume operator on a `Dim`-dimensional grid acting along the
+specified coordinate `direction`. The action of the operator is determined by the
+stencils `inner_stencil` and `closure_stencils`.
+When `Dim=1`, the corresponding `VolumeOperator` tensor mapping is returned.
+When `Dim>1`, the `VolumeOperator` `op` is inflated by the outer product
+of `IdentityMappings` in orthogonal coordinate directions, e.g for `Dim=3`,
+the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`.
+"""
+function volume_operator(grid::EquidistantGrid{Dim,T}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, parity, direction) where {Dim,T,M}
+    # Create 1D volume operator in along coordinate direction
+    op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity)
+    # Create 1D IdentityMappings for each coordinate direction
+    one_d_grids = restrict.(Ref(grid), Tuple(1:Dim))
+    Is = IdentityMapping{T}.(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)
+end
+export volume_operator
+
+"""
+    VolumeOperator{T,N,M,K} <: TensorOperator{T,1}
+Implements a one-dimensional constant coefficients volume operator
+"""
+struct VolumeOperator{T,N,M,K} <: TensorMapping{T,1,1}
+    inner_stencil::Stencil{T,N}
+    closure_stencils::NTuple{M,Stencil{T,K}}
+    size::NTuple{1,Int}
+    parity::Parity
+end
+export VolumeOperator
+
+function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity)
+    return VolumeOperator(inner_stencil, closure_stencils, size(grid), parity)
+end
+
+closure_size(::VolumeOperator{T,N,M}) where {T,N,M} = M
+
+LazyTensors.range_size(op::VolumeOperator) = op.size
+LazyTensors.domain_size(op::VolumeOperator) = op.size
+
+function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T
+    return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i))
+end
+
+function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T
+    return apply_stencil(op.inner_stencil, v, Int(i))
+end
+
+function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T
+    return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i))
+end
+
+function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T
+    r = getregion(i, closure_size(op), op.size[1])
+    return LazyTensors.apply(op, v, Index(i, r))
+end

--- a/test/testSbpOperators.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/test/testSbpOperators.jl	Mon Dec 07 12:07:29 2020 +0100
@@ -200,10 +200,10 @@
     # implies that L*v should be exact for v - monomial up to order 3.
     # Exact differentiation is measured point-wise. For other grid functions
     # the error is measured in the H-norm.
-    @test norm(L*v0) ≈ 0 atol=5e-10
-    @test norm(L*v1) ≈ 0 atol=5e-10
+    @test norm(L*v0) ≈ 0 atol=1e-9
+    @test norm(L*v1) ≈ 0 atol=1e-9
     @test L*v2 ≈ v0 # Seems to be more accurate
-    @test L*v3 ≈ v1 atol=5e-10
+    @test L*v3 ≈ v1 atol=1e-9
 
     h = spacing(g)
     l2(v) = sqrt(prod(h)*sum(v.^2))
@@ -274,58 +274,58 @@
     @test Qinv*v == Qinv'*v
 end
 
-@testset "BoundaryRestrictrion" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+@testset "BoundaryOperator" begin
+    closure_stencil = Stencil((0,2), (2.,1.,3.))
     g_1D = EquidistantGrid(11, 0.0, 1.0)
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "Constructors" begin
         @testset "1D" begin
-            e_l = BoundaryRestriction{Lower}(op.eClosure,size(g_1D)[1])
-            @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower())
-            @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
-            @test e_l isa TensorMapping{T,0,1} where T
+            op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1])
+            @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower())
+            @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}())
+            @test op_l isa TensorMapping{T,0,1} where T
 
-            e_r = BoundaryRestriction{Upper}(op.eClosure,size(g_1D)[1])
-            @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper())
-            @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
-            @test e_r isa TensorMapping{T,0,1} where T
+            op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
+            @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper())
+            @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
+            @test op_r isa TensorMapping{T,0,1} where T
         end
 
         @testset "2D" begin
-            e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}())
             @test e_w isa InflatedTensorMapping
             @test e_w isa TensorMapping{T,1,2} where T
         end
     end
 
-    e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
-    e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+    op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}())
+    op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}())
 
-    e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
-    e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
-    e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
-    e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
+    op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}())
+    op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}())
+    op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}())
+    op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}())
 
     @testset "Sizes" begin
         @testset "1D" begin
-            @test domain_size(e_l) == (11,)
-            @test domain_size(e_r) == (11,)
+            @test domain_size(op_l) == (11,)
+            @test domain_size(op_r) == (11,)
 
-            @test range_size(e_l) == ()
-            @test range_size(e_r) == ()
+            @test range_size(op_l) == ()
+            @test range_size(op_r) == ()
         end
 
         @testset "2D" begin
-            @test domain_size(e_w) == (11,15)
-            @test domain_size(e_e) == (11,15)
-            @test domain_size(e_s) == (11,15)
-            @test domain_size(e_n) == (11,15)
+            @test domain_size(op_w) == (11,15)
+            @test domain_size(op_e) == (11,15)
+            @test domain_size(op_s) == (11,15)
+            @test domain_size(op_n) == (11,15)
 
-            @test range_size(e_w) == (15,)
-            @test range_size(e_e) == (15,)
-            @test range_size(e_s) == (11,)
-            @test range_size(e_n) == (11,)
+            @test range_size(op_w) == (15,)
+            @test range_size(op_e) == (15,)
+            @test range_size(op_s) == (11,)
+            @test range_size(op_n) == (11,)
         end
     end
 
@@ -334,6 +334,126 @@
         @testset "1D" begin
             v = evalOn(g_1D,x->1+x^2)
             u = fill(3.124)
+            @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3]
+            @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2]
+            @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2]
+            @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
+            @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
+        end
+
+        @testset "2D" begin
+            v = rand(size(g_2D)...)
+            u = fill(3.124)
+            @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14
+            @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14
+            @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14
+            @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14
+
+
+           g_x = rand(size(g_2D)[1])
+           g_y = rand(size(g_2D)[2])
+
+           G_w = zeros(Float64, size(g_2D)...)
+           G_w[1,:] = 2*g_y
+           G_w[2,:] = g_y
+           G_w[3,:] = 3*g_y
+
+           G_e = zeros(Float64, size(g_2D)...)
+           G_e[end,:] = 2*g_y
+           G_e[end-1,:] = g_y
+           G_e[end-2,:] = 3*g_y
+
+           G_s = zeros(Float64, size(g_2D)...)
+           G_s[:,1] = 2*g_x
+           G_s[:,2] = g_x
+           G_s[:,3] = 3*g_x
+
+           G_n = zeros(Float64, size(g_2D)...)
+           G_n[:,end] = 2*g_x
+           G_n[:,end-1] = g_x
+           G_n[:,end-2] = 3*g_x
+
+           @test op_w'*g_y == G_w
+           @test op_e'*g_y == G_e
+           @test op_s'*g_x == G_s
+           @test op_n'*g_x == G_n
+       end
+
+       @testset "Regions" begin
+           u = fill(3.124)
+           @test (op_l'*u)[Index(1,Lower)] == 2*u[]
+           @test (op_l'*u)[Index(2,Lower)] == u[]
+           @test (op_l'*u)[Index(6,Interior)] == 0
+           @test (op_l'*u)[Index(10,Upper)] == 0
+           @test (op_l'*u)[Index(11,Upper)] == 0
+
+           @test (op_r'*u)[Index(1,Lower)] == 0
+           @test (op_r'*u)[Index(2,Lower)] == 0
+           @test (op_r'*u)[Index(6,Interior)] == 0
+           @test (op_r'*u)[Index(10,Upper)] == u[]
+           @test (op_r'*u)[Index(11,Upper)] == 2*u[]
+       end
+    end
+
+    @testset "Inferred" begin
+        v = ones(Float64, 11)
+        u = fill(1.)
+
+        @inferred apply(op_l, v)
+        @inferred apply(op_r, v)
+
+        @inferred apply_transpose(op_l, u, 4)
+        @inferred apply_transpose(op_l, u, Index(1,Lower))
+        @inferred apply_transpose(op_l, u, Index(2,Lower))
+        @inferred apply_transpose(op_l, u, Index(6,Interior))
+        @inferred apply_transpose(op_l, u, Index(10,Upper))
+        @inferred apply_transpose(op_l, u, Index(11,Upper))
+
+        @inferred apply_transpose(op_r, u, 4)
+        @inferred apply_transpose(op_r, u, Index(1,Lower))
+        @inferred apply_transpose(op_r, u, Index(2,Lower))
+        @inferred apply_transpose(op_r, u, Index(6,Interior))
+        @inferred apply_transpose(op_r, u, Index(10,Upper))
+        @inferred apply_transpose(op_r, u, Index(11,Upper))
+    end
+
+end
+
+@testset "BoundaryRestriction" begin
+    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+    g_1D = EquidistantGrid(11, 0.0, 1.0)
+    g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
+
+    @testset "Constructors" begin
+        @testset "1D" begin
+            e_l = BoundaryRestriction(g_1D,op.eClosure,Lower())
+            @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
+            @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
+            @test e_l isa BoundaryOperator{T,Lower} where T
+            @test e_l isa TensorMapping{T,0,1} where T
+
+            e_r = BoundaryRestriction(g_1D,op.eClosure,Upper())
+            @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
+            @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper())
+            @test e_r isa BoundaryOperator{T,Upper} where T
+            @test e_r isa TensorMapping{T,0,1} where T
+        end
+
+        @testset "2D" begin
+            e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            @test e_w isa InflatedTensorMapping
+            @test e_w isa TensorMapping{T,1,2} where T
+        end
+    end
+
+    @testset "Application" begin
+        @testset "1D" begin
+            e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
+            e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
+
+            v = evalOn(g_1D,x->1+x^2)
+            u = fill(3.124)
+
             @test (e_l*v)[] == v[1]
             @test (e_r*v)[] == v[end]
             @test (e_r*v)[1] == v[end]
@@ -342,6 +462,11 @@
         end
 
         @testset "2D" begin
+            e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
+            e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
+            e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
+            e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
+
             v = rand(11, 15)
             u = fill(3.124)
 
@@ -371,45 +496,7 @@
            @test e_s'*g_x == G_s
            @test e_n'*g_x == G_n
        end
-
-       @testset "Regions" begin
-           u = fill(3.124)
-           @test (e_l'*u)[Index(1,Lower)] == 3.124
-           @test (e_l'*u)[Index(2,Lower)] == 0
-           @test (e_l'*u)[Index(6,Interior)] == 0
-           @test (e_l'*u)[Index(10,Upper)] == 0
-           @test (e_l'*u)[Index(11,Upper)] == 0
-
-           @test (e_r'*u)[Index(1,Lower)] == 0
-           @test (e_r'*u)[Index(2,Lower)] == 0
-           @test (e_r'*u)[Index(6,Interior)] == 0
-           @test (e_r'*u)[Index(10,Upper)] == 0
-           @test (e_r'*u)[Index(11,Upper)] == 3.124
-       end
     end
-
-    @testset "Inferred" begin
-        v = ones(Float64, 11)
-        u = fill(1.)
-
-        @inferred apply(e_l, v)
-        @inferred apply(e_r, v)
-
-        @inferred apply_transpose(e_l, u, 4)
-        @inferred apply_transpose(e_l, u, Index(1,Lower))
-        @inferred apply_transpose(e_l, u, Index(2,Lower))
-        @inferred apply_transpose(e_l, u, Index(6,Interior))
-        @inferred apply_transpose(e_l, u, Index(10,Upper))
-        @inferred apply_transpose(e_l, u, Index(11,Upper))
-
-        @inferred apply_transpose(e_r, u, 4)
-        @inferred apply_transpose(e_r, u, Index(1,Lower))
-        @inferred apply_transpose(e_r, u, Index(2,Lower))
-        @inferred apply_transpose(e_r, u, Index(6,Interior))
-        @inferred apply_transpose(e_r, u, Index(10,Upper))
-        @inferred apply_transpose(e_r, u, Index(11,Upper))
-    end
-
 end
 #
 # @testset "NormalDerivative" begin