--- a/src/LazyTensors/lazy_tensor_operations.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -269,6 +269,29 @@
 LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms)
 
 
+
+"""
+    inflate(tm::LazyTensor, sz, dir)
+
+Inflate `tm` such that it gets the size `sz` in all directions except `dir`.
+Here `sz[dir]` is ignored and replaced with the range and domains size of
+`tm`.
+
+An example of when this operation is useful is when extending a one
+dimensional difference operator `D` to a 2D grid of a ceratin size. In that
+case we could have
+
+```julia
+Dx = inflate(D, (10,10), 1)
+Dy = inflate(D, (10,10), 2)
+```
+"""
+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/SbpOperators/boundaryops/boundary_operator.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -1,32 +1,3 @@
-"""
-    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 `IdentityTensors` 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, 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)
-
-    # Create 1D IdentityTensors for each coordinate direction
-    one_d_grids = restrict.(Ref(grid), Tuple(dims(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)
-end
-
 """
     BoundaryOperator{T,R,N} <: LazyTensor{T,0,1}
 
@@ -41,8 +12,6 @@
     size::Int
 end
 
-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)
 
@@ -55,6 +24,7 @@
 
 """
     closure_size(::BoundaryOperator)
+
 The size of the closure stencil.
 """
 closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N

--- a/src/SbpOperators/boundaryops/boundary_restriction.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -8,11 +8,13 @@
 On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of
 a `BoundaryOperator`.
 
-See also: [`boundary_operator`](@ref).
+See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function boundary_restriction(grid, closure_stencil, boundary)
     converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil)
-    return SbpOperators.boundary_operator(grid, converted_stencil, boundary)
+
+    op = BoundaryOperator(restrict(grid, dim(boundary)), converted_stencil, region(boundary))
+    return LazyTensors.inflate(op, size(grid), dim(boundary))
 end
 
 """

--- a/src/SbpOperators/boundaryops/normal_derivative.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -8,12 +8,14 @@
 On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
 a `BoundaryOperator`.
 
-See also: [`boundary_operator`](@ref).
+See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function normal_derivative(grid, closure_stencil, boundary)
     direction = dim(boundary)
     h_inv = inverse_spacing(grid)[direction]
-    return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
+
+    op = BoundaryOperator(restrict(grid, dim(boundary)), scale(closure_stencil,h_inv), region(boundary))
+    return LazyTensors.inflate(op, size(grid), dim(boundary))
 end
 
 """

--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -7,14 +7,16 @@
 `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
 for the points in the closure regions.
 
-On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the outer product of the
-one-dimensional operator with the `IdentityTensor`s in orthogonal coordinate dirrections.
+On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the inflation of
+a `VolumeOperator`.
 
-See also: [`volume_operator`](@ref).
+See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
     h_inv = inverse_spacing(grid)[direction]
-    return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction)
+
+    D₁ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd)
+    return LazyTensors.inflate(D₁, size(grid), direction)
 end
 
 

--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -7,14 +7,16 @@
 `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils`
 for the points in the closure regions.
 
-On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the outer product of the
-one-dimensional operator with the `IdentityTensor`s in orthogonal coordinate dirrections.
+On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the inflation of
+a `VolumeOperator`.
 
-See also: [`volume_operator`](@ref).
+See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref).
 """
 function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction)
     h_inv = inverse_spacing(grid)[direction]
-    return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction)
+
+    D₂ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even)
+    return LazyTensors.inflate(D₂, size(grid), direction)
 end
 
 

--- a/src/SbpOperators/volumeops/volume_operator.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -1,30 +1,6 @@
-"""
-    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
-returned operator is the appropriate outer product of a one-dimensional
-operators and `IdentityTensor`s, e.g for `Dim=3` the volume operator in the
-y-direction is `I⊗op⊗I`.
-"""
-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(dims(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)
-end
-
 """
     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}
@@ -59,3 +35,4 @@
     r = getregion(i, closure_size(op), op.size[1])
     return LazyTensors.apply(op, v, Index(i, r))
 end
+# TODO: Move this to LazyTensors when we have the region communication down.

--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -366,3 +366,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/SbpOperators/boundaryops/boundary_operator_test.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -6,7 +6,7 @@
 using Sbplib.RegionIndices
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.BoundaryOperator
-import Sbplib.SbpOperators.boundary_operator
+
 
 @testset "BoundaryOperator" begin
     closure_stencil = Stencil(2.,1.,3.; center = 1)
@@ -14,120 +14,50 @@
     g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
 
     @testset "Constructors" begin
-        @testset "1D" begin
-            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 LazyTensor{T,0,1} where T
+        @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
+    end
 
-            op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
-            @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper())
-            @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
-            @test op_r isa LazyTensor{T,0,1} where T
-        end
-
-        @testset "2D" begin
-            e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}())
-            @test e_w isa InflatedTensor
-            @test e_w isa LazyTensor{T,1,2} where T
-        end
-    end
-    op_l, op_r = boundary_operator.(Ref(g_1D), Ref(closure_stencil), boundary_identifiers(g_1D))
-    op_w, op_e, op_s, op_n = boundary_operator.(Ref(g_2D), Ref(closure_stencil), boundary_identifiers(g_2D))
+    op_l = BoundaryOperator(g_1D, closure_stencil, Lower())
+    op_r = BoundaryOperator(g_1D, closure_stencil, Upper())
 
     @testset "Sizes" begin
-        @testset "1D" begin
-            @test domain_size(op_l) == (11,)
-            @test domain_size(op_r) == (11,)
-
-            @test range_size(op_l) == ()
-            @test range_size(op_r) == ()
-        end
+        @test domain_size(op_l) == (11,)
+        @test domain_size(op_r) == (11,)
 
-        @testset "2D" begin
-            @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(op_w) == (15,)
-            @test range_size(op_e) == (15,)
-            @test range_size(op_s) == (11,)
-            @test range_size(op_n) == (11,)
-        end
+        @test range_size(op_l) == ()
+        @test range_size(op_r) == ()
     end
 
     @testset "Application" begin
-        @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[]]
-
-            v = evalOn(g_1D, x->1. +x*im)
-            @test (op_l*v)[] isa ComplexF64
+        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[]]
 
-            u = fill(1. +im)
-            @test (op_l'*u)[1] isa ComplexF64
-            @test (op_l'*u)[5] isa ComplexF64
-            @test (op_l'*u)[11] isa ComplexF64
-        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
+        v = evalOn(g_1D, x->1. +x*im)
+        @test (op_l*v)[] isa ComplexF64
 
-            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
+        u = fill(1. +im)
+        @test (op_l'*u)[1] isa ComplexF64
+        @test (op_l'*u)[5] isa ComplexF64
+        @test (op_l'*u)[11] isa ComplexF64
 
-            @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
+        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
 
-       @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
+        @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
 
     @testset "Inferred" begin

--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -6,6 +6,8 @@
 
 import Sbplib.SbpOperators.VolumeOperator
 
+# TODO: Refactor these test to look more like the tests in first_derivative_test.jl.
+
 @testset "SecondDerivative" begin
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
     stencil_set = read_stencil_set(operator_path; order=4)

--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Fri Oct 21 21:41:43 2022 +0200
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Tue Nov 29 22:11:01 2022 +0100
@@ -7,120 +7,76 @@
 
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.VolumeOperator
-import Sbplib.SbpOperators.volume_operator
 import Sbplib.SbpOperators.odd
 import Sbplib.SbpOperators.even
 
+
 @testset "VolumeOperator" begin
     inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
-    closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2))
-    g_1D = EquidistantGrid(11,0.,1.)
-    g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
-    g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
+    closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2))
+    g = EquidistantGrid(11,0.,1.)
+
     @testset "Constructors" begin
-        @testset "1D" begin
-            op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)
-            @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even)
-            @test op == volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
-            @test op isa LazyTensor{T,1,1} where T
-        end
-        @testset "2D" begin
-            op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
-            Ix = IdentityTensor{Float64}((11,))
-            Iy = IdentityTensor{Float64}((12,))
-            @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy
-            @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)
-            @test op_x isa LazyTensor{T,2,2} where T
-            @test op_y isa LazyTensor{T,2,2} where T
-        end
-        @testset "3D" begin
-            op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
-            op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
-            Ix = IdentityTensor{Float64}((11,))
-            Iy = IdentityTensor{Float64}((12,))
-            Iz = IdentityTensor{Float64}((10,))
-            @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy⊗Iz
-            @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)⊗Iz
-            @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),even)
-            @test op_x isa LazyTensor{T,3,3} where T
-            @test op_y isa LazyTensor{T,3,3} where T
-            @test op_z isa LazyTensor{T,3,3} where T
-        end
+        op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)
+        @test op == VolumeOperator(g,inner_stencil,closure_stencils,even)
+        @test op isa LazyTensor{T,1,1} where T
     end
 
     @testset "Sizes" begin
-        @testset "1D" begin
-            op = volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
-            @test range_size(op) == domain_size(op) == size(g_1D)
-        end
-
-        @testset "2D" begin
-            op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
-            @test range_size(op_y) == domain_size(op_y) ==
-                  range_size(op_x) == domain_size(op_x) == size(g_2D)
-        end
-        @testset "3D" begin
-            op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
-            op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
-            op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
-            @test range_size(op_z) == domain_size(op_z) ==
-                  range_size(op_y) == domain_size(op_y) ==
-                  range_size(op_x) == domain_size(op_x) == size(g_3D)
-        end
+        op = VolumeOperator(g,inner_stencil,closure_stencils,even)
+        @test range_size(op) == domain_size(op) == size(g)
     end
 
-    op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-    op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2)
-    v = zeros(size(g_2D))
-    Nx = size(g_2D)[1]
-    Ny = size(g_2D)[2]
-    for i = 1:Nx
-        v[i,:] .= i
+
+    op_even = VolumeOperator(g, inner_stencil, closure_stencils, even)
+    op_odd =  VolumeOperator(g, inner_stencil, closure_stencils, odd)
+
+    N = size(g)[1]
+    v = rand(N)
+
+    r_even = copy(v)
+    r_odd  = copy(v)
+
+    r_even[1] = (v[1] + v[2])/2
+    r_odd[1]  = (v[1] + v[2])/2
+
+    r_even[2] = 2v[1] + v[2]
+    r_odd[2]  = 2v[1] + v[2]
+
+    for i ∈ 3:N-2
+        r_even[i] = (v[i-1] + 2v[i] + v[i+1])/4
+        r_odd[i]  = (v[i-1] + 2v[i] + v[i+1])/4
     end
-    rx = copy(v)
-    rx[1,:] .= 1.5
-    rx[Nx,:] .= (2*Nx-1)/2
-    ry = copy(v)
-    ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny]
+
+    r_even[N-1] =  v[N-1] + 2v[N]
+    r_odd[N-1]  = -v[N-1] - 2v[N]
+
+    r_even[N] =  (v[N-1] + v[N])/2
+    r_odd[N]  = -(v[N-1] + v[N])/2
+
 
     @testset "Application" begin
-        @test op_x*v ≈ rx rtol = 1e-14
-        @test op_y*v ≈ ry rtol = 1e-14
+        @test op_even*v ≈ r_even
+        @test op_odd*v  ≈ r_odd
 
-        @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64
+        @test (op_even*rand(ComplexF64,size(g)))[2] isa ComplexF64
     end
 
     @testset "Regions" begin
-        @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14
-        @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14
-        @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14
-        @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14
-        @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14
-
-        @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)]
-        @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)]
+        @test (op_even*v)[Index(1,Lower)]    ≈ r_even[1]
+        @test (op_even*v)[Index(2,Lower)]    ≈ r_even[2]
+        @test (op_even*v)[Index(6,Interior)] ≈ r_even[6]
+        @test (op_even*v)[Index(10,Upper)]   ≈ r_even[10]
+        @test (op_even*v)[Index(11,Upper)]   ≈ r_even[11]
 
-        @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14
-        @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14
-
-        @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)]
-        @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)]
+        @test_throws BoundsError (op_even*v)[Index(3,Lower)]
+        @test_throws BoundsError (op_even*v)[Index(9,Upper)]
     end
 
     @testset "Inferred" begin
-        @test_skip @inferred apply(op_x, v,1,1)
-        @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower))
-        @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower))
-        @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower))
-        @test_skip @inferred apply(op_y, v,1,1)
-        @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower))
-        @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior))
-        @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper))
+        @inferred apply(op_even, v, 1)
+        @inferred apply(op_even, v, Index(1,Lower))
+        @inferred apply(op_even, v, Index(6,Interior))
+        @inferred apply(op_even, v, Index(11,Upper))
     end
 end