--- a/src/SbpOperators/SbpOperators.jl	Sat Dec 05 18:12:31 2020 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Sat Dec 05 19:14:39 2020 +0100
@@ -8,6 +8,7 @@
 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")

--- a/src/SbpOperators/laplace/laplace.jl	Sat Dec 05 18:12:31 2020 +0100
+++ b/src/SbpOperators/laplace/laplace.jl	Sat Dec 05 19:14:39 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}, innerStencil, closureStencils) where Dim
+    Δ = SecondDerivative(grid, innerStencil, closureStencils, 1)
+    for d = 2:Dim
+        Δ += SecondDerivative(grid, innerStencil, closureStencils, 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	Sat Dec 05 18:12:31 2020 +0100
+++ b/src/SbpOperators/laplace/secondderivative.jl	Sat Dec 05 19:14:39 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), size(grid), 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	Sat Dec 05 19:14:39 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