--- a/TODO.md	Sun Dec 06 10:53:15 2020 +0100
+++ b/TODO.md	Wed Jan 20 17:52:55 2021 +0100
@@ -5,22 +5,22 @@
  - [ ] Skriv tester
 
 ## Coding
- - [ ] Add new Laplace opertor to DiffOps, probably named WaveEqOp(?!!?)
- - [ ] Add 1D operators (D1, D2, e, d ... ) as TensorOperators
+ - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?)
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
  - [ ] Add a quick and simple way of running all tests for all subpackages.
- - [ ] Replace getindex hack for flatteing tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`)
+ - [ ] Replace getindex hack for flattening tuples with flatten_tuple. (eg. `getindex.(range_size.(L.D2),1)`)
  - [ ] Use `@inferred` in a lot of tests.
  - [ ] Make sure we are setting tolerances in tests in a consistent way
  - [ ] Add check for correct domain sizes to lazy tensor operations using SizeMismatch
- - [ ] Write down some coding guideline or checklist for code convetions. For example i,j,... for indecies and I for multi-index
+ - [ ] Write down some coding guideline or checklist for code conventions. For example i,j,... for indices and I for multi-index
  - [ ] Add boundschecking in TensorMappingApplication
  - [ ] Start renaming things in LazyTensors
  - [ ] Clean up RegionIndices
     1. [ ] Write tests for how things should work
     2. [ ] Update RegionIndices accordingly
     3. [ ] Fix the rest of the library
- - [ ] Add posibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`?
+    Should getregion also work for getregion(::Colon,...)
+ - [ ] Add possibility to create tensor mapping application with `()`, e.g `D1(v) <=> D1*v`?
 
 ## Repo
  - [ ] Add Vidar to the authors list
@@ -31,14 +31,3 @@
  - [ ] Kolla att vi gör boundschecks överallt och att de är markerade med @boundscheck
  - [ ] Kolla att vi har @inline på rätt ställen
  - [ ] Profilera
-
-
-# Old stuff todos (Are these still relevant?)
-Borde det finns motsvarande apply_stencil för apply_quadrature,
-apply_boundary_value och apply_normal_derivative?
-
-Borde man alltid skicka in N som parameter i apply_2nd_derivative, t.ex som i
-apply_quadrature?
-
-Just nu agerar apply_normal_derivative, apply_boundary_value på inte på v som
-en vektor, utan randvärdet plockas ut utanför. Känns inte konsistent med övrig design.

--- a/src/LazyTensors/lazy_tensor_operations.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -331,7 +331,7 @@
 split_tuple((1,2,3,4),Val(3)) -> (1,2,3), (4,)
 ```
 """
-function split_tuple(t::NTuple{N},::Val{M}) where {N,M}
+function split_tuple(t::NTuple{N,Any},::Val{M}) where {N,M}
     return slice_tuple(t,Val(1), Val(M)), slice_tuple(t,Val(M+1), Val(N))
 end
 
@@ -341,7 +341,7 @@
 Same as `split_tuple(t::NTuple{N},::Val{M})` but splits the tuple in three parts. With the first
 two parts having lenght `M` and `K`.
 """
-function split_tuple(t::NTuple{N},::Val{M},::Val{K}) where {N,M,K}
+function split_tuple(t::NTuple{N,Any},::Val{M},::Val{K}) where {N,M,K}
     p1, tail = split_tuple(t, Val(M))
     p2, p3 = split_tuple(tail, Val(K))
     return p1,p2,p3

--- a/src/SbpOperators/SbpOperators.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/SbpOperators.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -5,15 +5,15 @@
 using Sbplib.Grids
 
 include("stencil.jl")
-include("constantstenciloperator.jl")
 include("d2.jl")
 include("readoperator.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("volumeops/volume_operator.jl")
+include("volumeops/derivatives/secondderivative.jl")
+include("volumeops/laplace/laplace.jl")
+include("volumeops/quadratures/quadrature.jl")
+include("volumeops/quadratures/inverse_quadrature.jl")
+include("boundaryops/boundary_operator.jl")
 include("boundaryops/boundary_restriction.jl")
+include("boundaryops/normal_derivative.jl")
 
 end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundaryops/boundary_operator.jl	Wed Jan 20 17:52:55 2021 +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}
+    #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 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
+
+"""
+    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
+
+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	Wed Jan 20 17:52:55 2021 +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 `SbpOperators.boundary_operator(...)` for more details.
+"""
+BoundaryRestriction(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary) = SbpOperators.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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -0,0 +1,18 @@
+"""
+    NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region)
+
+Creates the normal derivative boundary operator `d` as a `TensorMapping`
+
+`d` is the normal derivative of a grid function at the boundary specified by `boundary` or `region` using some `closure_stencil`.
+`d'` is the prolongation of the normal derivative of a grid function to the whole grid using the same `closure_stencil`.
+On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of
+a `BoundaryOperator`. Also see the documentation of `SbpOperators.boundary_operator(...)` for more details.
+"""
+function NormalDerivative(grid::EquidistantGrid, closure_stencil::Stencil, boundary::CartesianBoundary)
+    direction = dim(boundary)
+    h_inv = inverse_spacing(grid)[direction]
+    return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary)
+end
+NormalDerivative(grid::EquidistantGrid{1}, closure_stencil::Stencil, region::Region) = NormalDerivative(grid, closure_stencil, CartesianBoundary{1,typeof(region)}())
+export NormalDerivative

--- a/src/SbpOperators/constantstenciloperator.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +0,0 @@
-abstract type ConstantStencilOperator end
-
-# Apply for different regions Lower/Interior/Upper or Unknown region
-@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Lower})
-    return @inbounds h_inv*h_inv*apply_stencil(op.closureStencils[Int(i)], v, Int(i))
-end
-
-@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Interior})
-    return @inbounds h_inv*h_inv*apply_stencil(op.innerStencil, v, Int(i))
-end
-
-@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Index{Upper})
-    N = length(v)
-    return @inbounds h_inv*h_inv*Int(op.parity)*apply_stencil_backwards(op.closureStencils[N-Int(i)+1], v, Int(i))
-end
-
-@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i)
-    N = length(v)
-    r = getregion(Int(index), closuresize(op), N)
-    i = Index(Int(index), r)
-    return apply_2nd_derivative(op, h_inv, v, i)
-end
-export apply_2nd_derivative
-
-apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Lower}, N::Integer) where T = v*h*op.quadratureClosure[Int(i)]
-apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Upper}, N::Integer) where T = v*h*op.quadratureClosure[N-Int(i)+1]
-apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h
-
-function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i, N::Integer) where T
-    r = getregion(i, closuresize(op), N)
-    return apply_quadrature(op, h, v, Index(i, r), N)
-end
-export apply_quadrature
-
-# TODO: Evaluate if divisions affect performance
-apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Lower}, N::Integer) where T = h_inv*v/op.quadratureClosure[Int(i)]
-apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Upper}, N::Integer) where T = h_inv*v/op.quadratureClosure[N-Int(i)+1]
-apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i::Index{Interior}, N::Integer) where T = v*h_inv
-
-function apply_inverse_quadrature(op::ConstantStencilOperator, h_inv::Real, v::T, i, N::Integer) where T
-    r = getregion(i, closuresize(op), N)
-    return apply_inverse_quadrature(op, h_inv, v, Index(i, r), N)
-end
-
-export apply_inverse_quadrature
-
-function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Lower})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    h_inv*apply_stencil(op.dClosure,v,1)
-end
-
-function apply_normal_derivative_transpose(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, ::Type{Upper})
-    @boundscheck if length(v) < closuresize(op)
-        throw(BoundsError())
-    end
-    -h_inv*apply_stencil_backwards(op.dClosure,v,length(v))
-end
-
-export apply_normal_derivative_transpose
-
-function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i, N::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(Int(i)) <= N)
-        throw(BoundsError())
-    end
-    h_inv*op.dClosure[Int(i)-1]*v
-end
-
-function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i, N::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(Int(i)) <= N)
-        throw(BoundsError())
-    end
-    -h_inv*op.dClosure[N-Int(i)]*v
-end
-
-export apply_normal_derivative

--- a/src/SbpOperators/d2.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/d2.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -5,15 +5,13 @@
     even = 1
 end
 
-struct D2{T,N,M,K} <: ConstantStencilOperator
-    quadratureClosure::NTuple{M,T}
-    innerStencil::Stencil{T,N}
-    closureStencils::NTuple{M,Stencil{T,K}}
-    eClosure::Stencil{T,M}
-    dClosure::Stencil{T,M}
+struct D2{T,M}
+    innerStencil::Stencil{T}
+    closureStencils::NTuple{M,Stencil{T}}
+    eClosure::Stencil{T}
+    dClosure::Stencil{T}
+    quadratureClosure::NTuple{M,Stencil{T}}
     parity::Parity
 end
 
-function closuresize(D::D2)::Int
-    return length(D.quadratureClosure)
-end
+closuresize(D::D2{T,M}) where {T,M} = M

--- a/src/SbpOperators/laplace/laplace.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,145 +0,0 @@
-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))
-    end
-
-    return Laplace(D2)
-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
-
-# quadrature(L::Laplace) = Quadrature(L.op, L.grid)
-# inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
-# boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
-# normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
-# boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
-# export NormalDerivative
-# """
-#     NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
-#
-# Implements the boundary operator `d` as a TensorMapping
-# """
-# struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
-#     op::D2{T,N,M,K}
-#     grid::EquidistantGrid{2}
-#     bId::CartesianBoundary
-# end
-#
-# # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
-# # Can we give special treatment to TensorMappings that go to a higher dim?
-# function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer})
-#     if dim(e.bId) == 1
-#         return (UnknownDim, domain_size[1])
-#     elseif dim(e.bId) == 2
-#         return (domain_size[1], UnknownDim)
-#     end
-# end
-# LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],)
-#
-# # TODO: Not type stable D:<
-# # TODO: Make this independent of dimension
-# function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
-#     i = I[dim(d.bId)]
-#     j = I[3-dim(d.bId)]
-#     N_i = size(d.grid)[dim(d.bId)]
-#     h_inv = inverse_spacing(d.grid)[dim(d.bId)]
-#     return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
-# end
-#
-# function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
-#     u = selectdim(v,3-dim(d.bId),Int(I[1]))
-#     return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
-# end
-#
-# """
-#     BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
-#
-# Implements the boundary operator `q` as a TensorOperator
-# """
-# export BoundaryQuadrature
-# struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
-#     op::D2{T,N,M,K}
-#     grid::EquidistantGrid{2}
-#     bId::CartesianBoundary
-# end
-#
-#
-# # TODO: Make this independent of dimension
-# function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T
-#     h = spacing(q.grid)[3-dim(q.bId)]
-#     N = size(v)
-#     return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
-# end
-#
-# LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I)
-#
-#
-#
-#
-# struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
-#
-# function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
-#     e = boundary_value(L, Bid())
-#     d = normal_derivative(L, Bid())
-#     Hᵧ = boundary_quadrature(L, Bid())
-#     H⁻¹ = inverse_quadrature(L)
-#     return (-H⁻¹*e*Hᵧ*(d'*v - g))[I]
-# end
-#
-# struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
-#     tau::Float64
-# end
-#
-# function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
-#     e = boundary_value(L, Bid())
-#     d = normal_derivative(L, Bid())
-#     Hᵧ = boundary_quadrature(L, Bid())
-#     H⁻¹ = inverse_quadrature(L)
-#     return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I]
-#     # Need to handle scalar multiplication and addition of TensorMapping
-# end
-
-# function apply(s::MyWaveEq{D},  v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
-    #   return apply(s.L, v, i) +
-# 		sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau),  v, s.g_w, i) +
-# 		sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau),  v, s.g_e, i) +
-# 		sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau),  v, s.g_s, i) +
-# 		sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau),  v, s.g_n, i)
-# end

--- a/src/SbpOperators/laplace/secondderivative.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,43 +0,0 @@
-"""
-    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))
-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

--- a/src/SbpOperators/quadrature/diagonal_inner_product.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,41 +0,0 @@
-export DiagonalInnerProduct, closuresize
-"""
-    DiagonalInnerProduct{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping
-"""
-struct DiagonalInnerProduct{T,M} <: TensorMapping{T,1,1}
-    h::T
-    quadratureClosure::NTuple{M,T}
-    size::Tuple{Int}
-end
-
-function DiagonalInnerProduct(g::EquidistantGrid{1}, quadratureClosure)
-    return DiagonalInnerProduct(spacing(g)[1], quadratureClosure, size(g))
-end
-
-LazyTensors.range_size(H::DiagonalInnerProduct) = H.size
-LazyTensors.domain_size(H::DiagonalInnerProduct) = H.size
-
-function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Lower}) where T
-    return @inbounds H.h*H.quadratureClosure[Int(i)]*v[Int(i)]
-end
-
-function LazyTensors.apply(H::DiagonalInnerProduct{T},v::AbstractVector{T}, i::Index{Upper}) where T
-    N = length(v);
-    return @inbounds H.h*H.quadratureClosure[N-Int(i)+1]*v[Int(i)]
-end
-
-function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Interior}) where T
-    return @inbounds H.h*v[Int(i)]
-end
-
-function LazyTensors.apply(H::DiagonalInnerProduct{T},  v::AbstractVector{T}, i) where T
-    N = length(v);
-    r = getregion(i, closuresize(H), N)
-    return LazyTensors.apply(H, v, Index(i, r))
-end
-
-LazyTensors.apply_transpose(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(H,v,i)
-
-closuresize(H::DiagonalInnerProduct{T,M}) where {T,M} = M

--- a/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,43 +0,0 @@
-export InverseDiagonalInnerProduct, closuresize
-"""
-    InverseDiagonalInnerProduct{Dim,T<:Real,M} <: TensorMapping{T,1,1}
-
-Implements the inverse diagonal inner product operator `Hi` of as a 1D TensorOperator
-"""
-struct InverseDiagonalInnerProduct{T<:Real,M} <: TensorMapping{T,1,1}
-    h_inv::T
-    inverseQuadratureClosure::NTuple{M,T}
-    size::Tuple{Int}
-end
-
-function InverseDiagonalInnerProduct(g::EquidistantGrid{1}, quadratureClosure)
-    return InverseDiagonalInnerProduct(inverse_spacing(g)[1], 1 ./ quadratureClosure, size(g))
-end
-
-LazyTensors.range_size(Hi::InverseDiagonalInnerProduct) = Hi.size
-LazyTensors.domain_size(Hi::InverseDiagonalInnerProduct) = Hi.size
-
-
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Lower}) where T
-    return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[Int(i)]*v[Int(i)]
-end
-
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Upper}) where T
-    N = length(v);
-    return @inbounds Hi.h_inv*Hi.inverseQuadratureClosure[N-Int(i)+1]*v[Int(i)]
-end
-
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i::Index{Interior}) where T
-    return @inbounds Hi.h_inv*v[Int(i)]
-end
-
-function LazyTensors.apply(Hi::InverseDiagonalInnerProduct{T},  v::AbstractVector{T}, i) where T
-    N = length(v);
-    r = getregion(i, closuresize(Hi), N)
-    return LazyTensors.apply(Hi, v, Index(i, r))
-end
-
-LazyTensors.apply_transpose(Hi::InverseDiagonalInnerProduct{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(Hi,v,i)
-
-
-closuresize(Hi::InverseDiagonalInnerProduct{T,M}) where {T,M} =  M

--- a/src/SbpOperators/quadrature/inverse_quadrature.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-export InverseQuadrature
-"""
-    InverseQuadrature{Dim,T<:Real,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the inverse quadrature operator `Qi` of Dim dimension as a TensorMapping
-The multi-dimensional tensor operator consists of a tuple of 1D InverseDiagonalInnerProduct
-tensor operators.
-"""
-struct InverseQuadrature{Dim,T<:Real,M} <: TensorMapping{T,Dim,Dim}
-    Hi::NTuple{Dim,InverseDiagonalInnerProduct{T,M}}
-end
-
-
-function InverseQuadrature(g::EquidistantGrid{Dim}, quadratureClosure) where Dim
-    Hi = ()
-    for i ∈ 1:Dim
-        Hi = (Hi..., InverseDiagonalInnerProduct(restrict(g,i), quadratureClosure))
-    end
-
-    return InverseQuadrature(Hi)
-end
-
-LazyTensors.range_size(Hi::InverseQuadrature) = getindex.(range_size.(Hi.Hi),1)
-LazyTensors.domain_size(Hi::InverseQuadrature) = getindex.(domain_size.(Hi.Hi),1)
-
-LazyTensors.domain_size(Qi::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
-
-function LazyTensors.apply(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim}
-    error("not implemented")
-end
-
-@inline function LazyTensors.apply(Qi::InverseQuadrature{1,T}, v::AbstractVector{T}, i) where T
-    @inbounds q = apply(Qi.Hi[1], v , i)
-    return q
-end
-
-@inline function LazyTensors.apply(Qi::InverseQuadrature{2,T}, v::AbstractArray{T,2}, i,j) where T
-    # InverseQuadrature in x direction
-    @inbounds vx = view(v, :, Int(j))
-    @inbounds qx_inv = apply(Qi.Hi[1], vx , i)
-    # InverseQuadrature in y-direction
-    @inbounds vy = view(v, Int(i), :)
-    @inbounds qy_inv = apply(Qi.Hi[2], vy, j)
-    return qx_inv*qy_inv
-end
-
-LazyTensors.apply_transpose(Qi::InverseQuadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {Dim,T} = LazyTensors.apply(Qi,v,I...)

--- a/src/SbpOperators/quadrature/quadrature.jl	Sun Dec 06 10:53:15 2020 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,44 +0,0 @@
-export Quadrature
-"""
-    Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the quadrature operator `Q` of Dim dimension as a TensorMapping
-The multi-dimensional tensor operator consists of a tuple of 1D DiagonalInnerProduct H
-tensor operators.
-"""
-struct Quadrature{Dim,T<:Real,M} <: TensorMapping{T,Dim,Dim}
-    H::NTuple{Dim,DiagonalInnerProduct{T,M}}
-end
-
-function Quadrature(g::EquidistantGrid{Dim}, quadratureClosure) where Dim
-    H = ()
-    for i ∈ 1:Dim
-        H = (H..., DiagonalInnerProduct(restrict(g,i), quadratureClosure))
-    end
-
-    return Quadrature(H)
-end
-
-LazyTensors.range_size(H::Quadrature) = getindex.(range_size.(H.H),1)
-LazyTensors.domain_size(H::Quadrature) = getindex.(domain_size.(H.H),1)
-
-function LazyTensors.apply(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {T,Dim}
-    error("not implemented")
-end
-
-function LazyTensors.apply(Q::Quadrature{1,T}, v::AbstractVector{T}, i) where T
-    @inbounds q = apply(Q.H[1], v , i)
-    return q
-end
-
-function LazyTensors.apply(Q::Quadrature{2,T}, v::AbstractArray{T,2}, i, j) where T
-    # Quadrature in x direction
-    @inbounds vx = view(v, :, Int(j))
-    @inbounds qx = apply(Q.H[1], vx , i)
-    # Quadrature in y-direction
-    @inbounds vy = view(v, Int(i), :)
-    @inbounds qy = apply(Q.H[2], vy, j)
-    return qx*qy
-end
-
-LazyTensors.apply_transpose(Q::Quadrature{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Any,Dim}) where {Dim,T} = LazyTensors.apply(Q,v,I...)

--- a/src/SbpOperators/readoperator.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/src/SbpOperators/readoperator.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -9,7 +9,6 @@
 export get_stencils
 export get_tuple
 
-
 function read_D2_operator(fn; order)
     operators = TOML.parsefile(fn)["order$order"]
     D2 = operators["D2"]
@@ -23,22 +22,25 @@
     # Create boundary stencils
     boundarySize = length(D2["closure_stencils"])
     closureStencils = Vector{typeof(innerStencil)}() # TBD: is the the right way to get the correct type?
-
     for i ∈ 1:boundarySize
         closureStencils = (closureStencils..., get_stencil(operators, "D2", "closure_stencils", i; center=i))
     end
-
     # TODO: Get rid of the padding here. Any padding should be handled by the consturctor accepting the stencils.
-    quadratureClosure = pad_tuple(toml_string_array_to_tuple(Float64, H["closure"]), boundarySize)
     eClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, e["closure"]), boundarySize), center=1)
     dClosure = Stencil(pad_tuple(toml_string_array_to_tuple(Float64, d1["closure"]), boundarySize), center=1)
 
+    q_tuple = pad_tuple(toml_string_array_to_tuple(Float64, H["closure"]), boundarySize)
+    quadratureClosure = Vector{typeof(innerStencil)}()
+    for i ∈ 1:boundarySize
+        quadratureClosure = (quadratureClosure..., Stencil((q_tuple[i],), center=1))
+    end
+
     d2 = SbpOperators.D2(
-        quadratureClosure,
         innerStencil,
         closureStencils,
         eClosure,
         dClosure,
+        quadratureClosure,
         even
     )
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/derivatives/secondderivative.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -0,0 +1,20 @@
+"""
+    SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction)
+    SecondDerivative(grid::EquidistantGrid{1}, inner_stencil, closure_stencils)
+
+Creates the second-derivative operator `D2` as a `TensorMapping`
+
+`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by
+`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 `IdentityMapping`s in orthogonal coordinate dirrections.
+Also see the documentation of `SbpOperators.volume_operator(...)` for more details.
+"""
+function SecondDerivative(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils, direction) where Dim
+    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)
+end
+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/laplace/laplace.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -0,0 +1,20 @@
+"""
+    Laplace(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils)
+
+Creates the Laplace operator operator `Δ` as a `TensorMapping`
+
+`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, 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`, `Δ` is a `SecondDerivative`. On a multi-dimensional `grid`, `Δ` is the sum of
+multi-dimensional `SecondDerivative`s where the sum is carried out lazily.
+"""
+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 Δ
+end
+export Laplace

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -0,0 +1,41 @@
+
+"""
+    InverseQuadrature(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
+
+Creates the inverse `H⁻¹` of the quadrature operator as a `TensorMapping`
+
+The inverse quadrature approximates the integral operator on the grid using
+`inv_inner_stencil` in the interior and a set of stencils `inv_closure_stencils`
+for the points in the closure regions.
+
+On a one-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a multi-dimensional
+`grid`, `H` is the outer product of the 1-dimensional inverse quadrature operators in
+each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details.
+"""
+function InverseQuadrature(grid::EquidistantGrid{Dim}, inv_inner_stencil, inv_closure_stencils) where Dim
+    h⁻¹ = inverse_spacing(grid)
+    H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1)
+    for i ∈ 2:Dim
+        Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i)
+        H⁻¹ = H⁻¹∘Hᵢ⁻¹
+    end
+    return H⁻¹
+end
+export InverseQuadrature
+
+"""
+    InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
+
+Creates the inverse of the diagonal quadrature operator defined by the inner stencil
+1/h and a set of 1-element closure stencils in `closure_stencils`. Note that
+the closure stencils are those of the quadrature operator (and not the inverse).
+"""
+function InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {T,M}
+    inv_inner_stencil = Stencil(Tuple{T}(1),center=1)
+    inv_closure_stencils = reciprocal_stencil.(closure_stencils)
+    return InverseQuadrature(grid, inv_inner_stencil, inv_closure_stencils)
+end
+export InverseDiagonalQuadrature
+
+reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/quadratures/quadrature.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -0,0 +1,36 @@
+"""
+    Quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils)
+
+Creates the quadrature operator `H` as a `TensorMapping`
+
+The quadrature approximates the integral operator on the grid using
+`inner_stencil` in the interior and a set of stencils `closure_stencils`
+for the points in the closure regions.
+
+On a one-dimensional `grid`, `H` is a `VolumeOperator`. On a multi-dimensional
+`grid`, `H` is the outer product of the 1-dimensional quadrature operators in
+each coordinate direction. Also see the documentation of
+`SbpOperators.volume_operator(...)` for more details.
+"""
+function Quadrature(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
+    h = spacing(grid)
+    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
+    for i ∈ 2:Dim
+        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
+        H = H∘Hᵢ
+    end
+    return H
+end
+export Quadrature
+
+"""
+    DiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
+
+Creates the quadrature operator with the inner stencil 1/h and 1-element sized
+closure stencils (i.e the operator is diagonal)
+"""
+function DiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T}
+    inner_stencil = Stencil(Tuple{T}(1),center=1)
+    return Quadrature(grid, inner_stencil, closure_stencils)
+end
+export DiagonalQuadrature

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/volumeops/volume_operator.jl	Wed Jan 20 17:52:55 2021 +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}
+    #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 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
+
+"""
+    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
+
+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/testLazyTensors.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/test/testLazyTensors.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -487,17 +487,22 @@
         @test LazyTensors.split_tuple((1,2,3,4),Val(3)) == ((1,2,3),(4,))
         @test LazyTensors.split_tuple((1,2,3,4),Val(4)) == ((1,2,3,4),())
 
+        @test LazyTensors.split_tuple((1,2,true,4),Val(3)) == ((1,2,true),(4,))
+
         @inferred LazyTensors.split_tuple((1,2,3,4),Val(3))
+        @inferred LazyTensors.split_tuple((1,2,true,4),Val(3))
     end
 
     @testset "3 parts" begin
         @test LazyTensors.split_tuple((),Val(0),Val(0)) == ((),(),())
         @test LazyTensors.split_tuple((1,2,3),Val(1), Val(1)) == ((1,),(2,),(3,))
+        @test LazyTensors.split_tuple((1,true,3),Val(1), Val(1)) == ((1,),(true,),(3,))
 
         @test LazyTensors.split_tuple((1,2,3,4,5,6),Val(1),Val(2)) == ((1,),(2,3),(4,5,6))
         @test LazyTensors.split_tuple((1,2,3,4,5,6),Val(3),Val(2)) == ((1,2,3),(4,5),(6,))
 
         @inferred LazyTensors.split_tuple((1,2,3,4,5,6),Val(3),Val(2))
+        @inferred LazyTensors.split_tuple((1,true,3),Val(1), Val(1))
     end
 end
 

--- a/test/testSbpOperators.jl	Sun Dec 06 10:53:15 2020 +0100
+++ b/test/testSbpOperators.jl	Wed Jan 20 17:52:55 2021 +0100
@@ -7,6 +7,13 @@
 using TOML
 
 import Sbplib.SbpOperators.Stencil
+import Sbplib.SbpOperators.VolumeOperator
+import Sbplib.SbpOperators.volume_operator
+import Sbplib.SbpOperators.BoundaryOperator
+import Sbplib.SbpOperators.boundary_operator
+import Sbplib.SbpOperators.even
+import Sbplib.SbpOperators.odd
+
 
 @testset "SbpOperators" begin
 
@@ -108,240 +115,636 @@
     end
 end
 
-# @testset "apply_quadrature" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     h = 0.5
-#
-#     @test apply_quadrature(op, h, 1.0, 10, 100) == h
-#
-#     N = 10
-#     qc = op.quadratureClosure
-#     q = h.*(qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
-#     @assert length(q) == N
-#
-#     for i ∈ 1:N
-#         @test apply_quadrature(op, h, 1.0, i, N) == q[i]
-#     end
-#
-#     v = [2.,3.,2.,4.,5.,4.,3.,4.,5.,4.5]
-#     for i ∈ 1:N
-#         @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i]
-#     end
-# end
+@testset "VolumeOperator" begin
+    inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2)
+    closure_stencils = (Stencil(1/2 .* (1.,1.),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.))
+    @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 TensorMapping{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 = IdentityMapping{Float64}((11,))
+            Iy = IdentityMapping{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 TensorMapping{T,2,2} where T
+            @test op_y isa TensorMapping{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 = IdentityMapping{Float64}((11,))
+            Iy = IdentityMapping{Float64}((12,))
+            Iz = IdentityMapping{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 TensorMapping{T,3,3} where T
+            @test op_y isa TensorMapping{T,3,3} where T
+            @test op_z isa TensorMapping{T,3,3} where T
+        end
+    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
+    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
+    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]
+
+    @testset "Application" begin
+        @test op_x*v ≈ rx rtol = 1e-14
+        @test op_y*v ≈ ry rtol = 1e-14
+    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_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)]
+    end
+
+    @testset "Inferred" begin
+        @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))
+
+        @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))
+    end
+
+end
 
 @testset "SecondDerivative" begin
     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    L = 3.5
-    g = EquidistantGrid(101, 0.0, L)
-    Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils)
+    Lx = 3.5
+    Ly = 3.
+    g_1D = EquidistantGrid(121, 0.0, Lx)
+    g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
+
+    @testset "Constructors" begin
+        @testset "1D" begin
+            Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+            @test Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1)
+            @test Dₓₓ isa VolumeOperator
+        end
+        @testset "2D" begin
+            Dₓₓ = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,1)
+            D2 = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+            I = IdentityMapping{Float64}(size(g_2D)[2])
+            @test Dₓₓ == D2⊗I
+            @test Dₓₓ isa TensorMapping{T,2,2} where T
+        end
+    end
+
+    # Exact differentiation is measured point-wise. In other cases
+    # the error is measured in the l2-norm.
+    @testset "Accuracy" begin
+        @testset "1D" begin
+            l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
+            monomials = ()
+            maxOrder = 4;
+            for i = 0:maxOrder-1
+                f_i(x) = 1/factorial(i)*x^i
+                monomials = (monomials...,evalOn(g_1D,f_i))
+            end
+            v = evalOn(g_1D,x -> sin(x))
+            vₓₓ = evalOn(g_1D,x -> -sin(x))
 
-    f0(x) = 1.
-    f1(x) = x
-    f2(x) = 1/2*x^2
-    f3(x) = 1/6*x^3
-    f4(x) = 1/24*x^4
-    f5(x) = sin(x)
-    f5ₓₓ(x) = -f5(x)
+            # 2nd order interior stencil, 1nd order boundary stencil,
+            # implies that L*v should be exact for monomials up to order 2.
+            @testset "2nd order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+                @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
+                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
+            end
+
+            # 4th order interior stencil, 2nd order boundary stencil,
+            # implies that L*v should be exact for monomials up to order 3.
+            @testset "4th order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
+                # NOTE: high tolerances for checking the "exact" differentiation
+                # due to accumulation of round-off errors/cancellation errors?
+                @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
+                @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
+                @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
+                @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
+            end
+        end
+        
+        @testset "2D" begin
+            l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
+            binomials = ()
+            maxOrder = 4;
+            for i = 0:maxOrder-1
+                f_i(x,y) = 1/factorial(i)*y^i + x^i
+                binomials = (binomials...,evalOn(g_2D,f_i))
+            end
+            v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
+            v_yy = evalOn(g_2D,(x,y) -> -cos(y))
 
-    v0 = evalOn(g,f0)
-    v1 = evalOn(g,f1)
-    v2 = evalOn(g,f2)
-    v3 = evalOn(g,f3)
-    v4 = evalOn(g,f4)
-    v5 = evalOn(g,f5)
+            # 2nd order interior stencil, 1st order boundary stencil,
+            # implies that L*v should be exact for binomials up to order 2.
+            @testset "2nd order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
+                @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
+                @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
+            end
 
-    @test Dₓₓ isa TensorMapping{T,1,1} where T
-    @test Dₓₓ' isa TensorMapping{T,1,1} where T
+            # 4th order interior stencil, 2nd order boundary stencil,
+            # implies that L*v should be exact for binomials up to order 3.
+            @testset "4th order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
+                # NOTE: high tolerances for checking the "exact" differentiation
+                # due to accumulation of round-off errors/cancellation errors?
+                @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
+                @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
+                @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
+                @test Dyy*v ≈ v_yy rtol = 5e-4 norm = l2
+            end
+        end
+    end
+end
 
-    # 4th order interior stencil, 2nd order boundary stencil,
-    # implies that L*v should be exact for v - monomial up to order 3.
-    # Exact differentiation is measured point-wise. For other grid functions
+@testset "Laplace" begin
+    g_1D = EquidistantGrid(101, 0.0, 1.)
+    g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+    @testset "Constructors" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            L = Laplace(g_1D, op.innerStencil, op.closureStencils)
+            @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils)
+            @test L isa TensorMapping{T,1,1}  where T
+        end
+        @testset "3D" begin
+            L = Laplace(g_3D, op.innerStencil, op.closureStencils)
+            @test L isa TensorMapping{T,3,3} where T
+            Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1)
+            Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2)
+            Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3)
+            @test L == Dxx + Dyy + Dzz
+        end
+    end
+
+    # Exact differentiation is measured point-wise. In other cases
     # the error is measured in the l2-norm.
-    @test norm(Dₓₓ*v0) ≈ 0.0 atol=5e-10
-    @test norm(Dₓₓ*v1) ≈ 0.0 atol=5e-10
-    @test Dₓₓ*v2 ≈ v0 atol=5e-11
-    @test Dₓₓ*v3 ≈ v1 atol=5e-11
+    @testset "Accuracy" begin
+        l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2));
+        polynomials = ()
+        maxOrder = 4;
+        for i = 0:maxOrder-1
+            f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i)
+            polynomials = (polynomials...,evalOn(g_3D,f_i))
+        end
+        v = evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
+        Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
 
-    h = spacing(g)[1];
-    l2(v) = sqrt(h*sum(v.^2))
-    @test Dₓₓ*v4 ≈ v2  atol=5e-4 norm=l2
-    @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2
+        # 2nd order interior stencil, 1st order boundary stencil,
+        # implies that L*v should be exact for binomials up to order 2.
+        @testset "2nd order" begin
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            L = Laplace(g_3D,op.innerStencil,op.closureStencils)
+            @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
+            @test L*v ≈ Δv rtol = 5e-2 norm = l2
+        end
+
+        # 4th order interior stencil, 2nd order boundary stencil,
+        # implies that L*v should be exact for binomials up to order 3.
+        @testset "4th order" begin
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            L = Laplace(g_3D,op.innerStencil,op.closureStencils)
+            # NOTE: high tolerances for checking the "exact" differentiation
+            # due to accumulation of round-off errors/cancellation errors?
+            @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
+            @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
+            @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9
+            @test L*v ≈ Δv rtol = 5e-4 norm = l2
+        end
+    end
 end
 
+@testset "DiagonalQuadrature" begin
+    Lx = π/2.
+    Ly = Float64(π)
+    g_1D = EquidistantGrid(77, 0.0, Lx)
+    g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    integral(H,v) = sum(H*v)
+    @testset "Constructors" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+            inner_stencil = Stencil((1.,),center=1)
+            @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure)
+            @test H isa TensorMapping{T,1,1} where T
+        end
+        @testset "1D" begin
+            H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+            H_x = DiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
+            H_y = DiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
+            @test H == H_x⊗H_y
+            @test H isa TensorMapping{T,2,2} where T
+        end
+    end
 
-@testset "Laplace2D" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    Lx = 1.5
-    Ly = 3.2
-    g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
-    L = Laplace(g, op.innerStencil, op.closureStencils)
+    @testset "Sizes" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+            @test domain_size(H) == size(g_1D)
+            @test range_size(H) == size(g_1D)
+        end
+        @testset "2D" begin
+            H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+            @test domain_size(H) == size(g_2D)
+            @test range_size(H) == size(g_2D)
+        end
+    end
+
+    @testset "Accuracy" begin
+        @testset "1D" begin
+            v = ()
+            for i = 0:4
+                f_i(x) = 1/factorial(i)*x^i
+                v = (v...,evalOn(g_1D,f_i))
+            end
+            u = evalOn(g_1D,x->sin(x))
+
+            @testset "2nd order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                for i = 1:2
+                    @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
+                end
+                @test integral(H,u) ≈ 1. rtol = 1e-4
+            end
+
+            @testset "4th order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                for i = 1:4
+                    @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
+                end
+                @test integral(H,u) ≈ 1. rtol = 1e-8
+            end
+        end
+
+        @testset "2D" begin
+            b = 2.1
+            v = b*ones(Float64, size(g_2D))
+            u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
+            @testset "2nd order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
+                @test integral(H,u) ≈ π rtol = 1e-4
+            end
+            @testset "4th order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
+                @test integral(H,u) ≈ π rtol = 1e-8
+            end
+        end
+    end
+end
+
+@testset "InverseDiagonalQuadrature" begin
+    Lx = π/2.
+    Ly = Float64(π)
+    g_1D = EquidistantGrid(77, 0.0, Lx)
+    g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    @testset "Constructors" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure);
+            inner_stencil = Stencil((1.,),center=1)
+            closures = ()
+            for i = 1:length(op.quadratureClosure)
+                closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
+            end
+            @test Hi == InverseQuadrature(g_1D,inner_stencil,closures)
+            @test Hi isa TensorMapping{T,1,1} where T
+        end
+        @testset "2D" begin
+            Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+            Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
+            Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
+            @test Hi == Hi_x⊗Hi_y
+            @test Hi isa TensorMapping{T,2,2} where T
+        end
+    end
+
+    @testset "Sizes" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+            @test domain_size(Hi) == size(g_1D)
+            @test range_size(Hi) == size(g_1D)
+        end
+        @testset "2D" begin
+            Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+            @test domain_size(Hi) == size(g_2D)
+            @test range_size(Hi) == size(g_2D)
+        end
+    end
+
+    @testset "Accuracy" begin
+        @testset "1D" begin
+            v = evalOn(g_1D,x->sin(x))
+            u = evalOn(g_1D,x->x^3-x^2+1)
+            @testset "2nd order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+                @test Hi*H*v ≈ v rtol = 1e-15
+                @test Hi*H*u ≈ u rtol = 1e-15
+            end
+            @testset "4th order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+                @test Hi*H*v ≈ v rtol = 1e-15
+                @test Hi*H*u ≈ u rtol = 1e-15
+            end
+        end
+        @testset "2D" begin
+            v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
+            u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
+            @testset "2nd order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+                @test Hi*H*v ≈ v rtol = 1e-15
+                @test Hi*H*u ≈ u rtol = 1e-15
+            end
+            @testset "4th order" begin
+                op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
+                @test Hi*H*v ≈ v rtol = 1e-15
+                @test Hi*H*u ≈ u rtol = 1e-15
+            end
+        end
+    end
+end
+
+@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
+            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
+
+            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_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
+
+    op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}())
+    op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,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(op_l) == (11,)
+            @test domain_size(op_r) == (11,)
+
+            @test range_size(op_l) == ()
+            @test range_size(op_r) == ()
+        end
+
+        @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
+    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[]]
+        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
 
 
-    f0(x,y) = 2.
-    f1(x,y) = x+y
-    f2(x,y) = 1/2*x^2 + 1/2*y^2
-    f3(x,y) = 1/6*x^3 + 1/6*y^3
-    f4(x,y) = 1/24*x^4 + 1/24*y^4
-    f5(x,y) = sin(x) + cos(y)
-    f5ₓₓ(x,y) = -f5(x,y)
+            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
 
-    v0 = evalOn(g,f0)
-    v1 = evalOn(g,f1)
-    v2 = evalOn(g,f2)
-    v3 = evalOn(g,f3)
-    v4 = evalOn(g,f4)
-    v5 = evalOn(g,f5)
-    v5ₓₓ = evalOn(g,f5ₓₓ)
+            G_s = zeros(Float64, size(g_2D)...)
+            G_s[:,1] = 2*g_x
+            G_s[:,2] = g_x
+            G_s[:,3] = 3*g_x
 
-    @test L isa TensorMapping{T,2,2} where T
-    @test L' isa TensorMapping{T,2,2} where T
+            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
 
-    # 4th order interior stencil, 2nd order boundary stencil,
-    # 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 L*v2 ≈ v0 # Seems to be more accurate
-    @test L*v3 ≈ v1 atol=5e-10
+       @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
 
-    h = spacing(g)
-    l2(v) = sqrt(prod(h)*sum(v.^2))
-    @test L*v4 ≈ v2   atol=5e-4 norm=l2
-    @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2
-end
+    @testset "Inferred" begin
+        v = ones(Float64, 11)
+        u = fill(1.)
+
+        @inferred apply(op_l, v)
+        @inferred apply(op_r, v)
 
-@testset "DiagonalInnerProduct" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    L = 2.3
-    g = EquidistantGrid(77, 0.0, L)
-    H = DiagonalInnerProduct(g,op.quadratureClosure)
-    v = ones(Float64, size(g))
+        @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))
 
-    @test H isa TensorMapping{T,1,1} where T
-    @test H' isa TensorMapping{T,1,1} where T
-    @test sum(H*v) ≈ L
-    @test H*v == H'*v
+        @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 "Quadrature" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    Lx = 2.3
-    Ly = 5.2
-    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-
-    Q = Quadrature(g, op.quadratureClosure)
-
-    @test Q isa TensorMapping{T,2,2} where T
-    @test Q' isa TensorMapping{T,2,2} where T
-
-    v = ones(Float64, size(g))
-    @test sum(Q*v) ≈ Lx*Ly
-
-    v = 2*ones(Float64, size(g))
-    @test_broken sum(Q*v) ≈ 2*Lx*Ly
-
-    @test Q*v == Q'*v
-end
-
-@testset "InverseDiagonalInnerProduct" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    L = 2.3
-    g = EquidistantGrid(77, 0.0, L)
-    H = DiagonalInnerProduct(g, op.quadratureClosure)
-    Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure)
-    v = evalOn(g, x->sin(x))
-
-    @test Hi isa TensorMapping{T,1,1} where T
-    @test Hi' isa TensorMapping{T,1,1} where T
-    @test Hi*H*v ≈ v
-    @test Hi*v == Hi'*v
-end
-
-@testset "InverseQuadrature" begin
-    op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-    Lx = 7.3
-    Ly = 8.2
-    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-
-    Q = Quadrature(g, op.quadratureClosure)
-    Qinv = InverseQuadrature(g, op.quadratureClosure)
-    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-
-    @test Qinv isa TensorMapping{T,2,2} where T
-    @test Qinv' isa TensorMapping{T,2,2} where T
-    @test_broken Qinv*(Q*v) ≈ v
-    @test Qinv*v == Qinv'*v
-end
-
-@testset "BoundaryRestrictrion" begin
+@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{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}())
+            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{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}())
+            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 = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
+            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
 
-    e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
-    e_r = boundary_restriction(g_1D, op.eClosure, 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}())
+    @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}())
 
-    @testset "Sizes" begin
-        @testset "1D" begin
-            @test domain_size(e_l) == (11,)
-            @test domain_size(e_r) == (11,)
+            v = evalOn(g_1D,x->1+x^2)
+            u = fill(3.124)
 
-            @test range_size(e_l) == ()
-            @test range_size(e_r) == ()
+            @test (e_l*v)[] == v[1]
+            @test (e_r*v)[] == v[end]
+            @test (e_r*v)[1] == v[end]
         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 range_size(e_w) == (15,)
-            @test range_size(e_e) == (15,)
-            @test range_size(e_s) == (11,)
-            @test range_size(e_n) == (11,)
-        end
-    end
-
+            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}())
 
-    @testset "Application" begin
-        @testset "1D" begin
-            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]
-            @test e_l'*u == [u[]; zeros(10)]
-            @test e_r'*u == [zeros(10); u[]]
-        end
-
-        @testset "2D" begin
             v = rand(11, 15)
             u = fill(3.124)
 
@@ -349,194 +752,66 @@
             @test e_e*v == v[end,:]
             @test e_s*v == v[:,1]
             @test e_n*v == v[:,end]
-
-
-           g_x = rand(11)
-           g_y = rand(15)
-
-           G_w = zeros(Float64, (11,15))
-           G_w[1,:] = g_y
-
-           G_e = zeros(Float64, (11,15))
-           G_e[end,:] = g_y
-
-           G_s = zeros(Float64, (11,15))
-           G_s[:,1] = g_x
-
-           G_n = zeros(Float64, (11,15))
-           G_n[:,end] = g_x
-
-           @test e_w'*g_y == G_w
-           @test e_e'*g_y == G_e
-           @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)
+end
 
-        @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))
+@testset "NormalDerivative" begin
+    g_1D = EquidistantGrid(11, 0.0, 1.0)
+    g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
+    @testset "Constructors" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            d_l = NormalDerivative(g_1D, op.dClosure, Lower())
+            @test d_l == NormalDerivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
+            @test d_l isa BoundaryOperator{T,Lower} where T
+            @test d_l isa TensorMapping{T,0,1} where T
+        end
+        @testset "2D" begin
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+            Ix = IdentityMapping{Float64}((size(g_2D)[1],))
+            Iy = IdentityMapping{Float64}((size(g_2D)[2],))
+            d_l = NormalDerivative(restrict(g_2D,1),op.dClosure,Lower())
+            d_r = NormalDerivative(restrict(g_2D,2),op.dClosure,Upper())
+            @test d_w ==  d_l⊗Iy
+            @test d_n ==  Ix⊗d_r
+            @test d_w isa TensorMapping{T,1,2} where T
+            @test d_n isa TensorMapping{T,1,2} where T
+        end
+    end
+    @testset "Accuracy" begin
+        v = evalOn(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y)
+        v∂x = evalOn(g_2D, (x,y)-> 2*x + y)
+        v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x)
+        # TODO: Test for higher order polynomials?
+        @testset "2nd order" begin
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
+            d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
+            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,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))
+            @test d_w*v ≈ v∂x[1,:] atol = 1e-13
+            @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
+            @test d_s*v ≈ v∂y[:,1] atol = 1e-13
+            @test d_n*v ≈ -v∂y[:,end] atol = 1e-13
+        end
+
+        @testset "4th order" begin
+            op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+            d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
+            d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
+            d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
+            d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
+
+            @test d_w*v ≈ v∂x[1,:] atol = 1e-13
+            @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
+            @test d_s*v ≈ v∂y[:,1] atol = 1e-13
+            @test d_n*v ≈ -v∂y[:,end] atol = 1e-13
+        end
     end
+end
 
 end
-#
-# @testset "NormalDerivative" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
-#
-#     d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
-#     d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
-#     d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
-#     d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
-#
-#
-#     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-#     v∂x = evalOn(g, (x,y)-> 2*x + y)
-#     v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
-#
-#     @test d_w  isa TensorMapping{T,2,1} where T
-#     @test d_w' isa TensorMapping{T,1,2} where T
-#
-#     @test domain_size(d_w, (3,2)) == (2,)
-#     @test domain_size(d_e, (3,2)) == (2,)
-#     @test domain_size(d_s, (3,2)) == (3,)
-#     @test domain_size(d_n, (3,2)) == (3,)
-#
-#     @test size(d_w'*v) == (6,)
-#     @test size(d_e'*v) == (6,)
-#     @test size(d_s'*v) == (5,)
-#     @test size(d_n'*v) == (5,)
-#
-#     @test d_w'*v .≈ v∂x[1,:]
-#     @test d_e'*v .≈ v∂x[5,:]
-#     @test d_s'*v .≈ v∂y[:,1]
-#     @test d_n'*v .≈ v∂y[:,6]
-#
-#
-#     d_x_l = zeros(Float64, 5)
-#     d_x_u = zeros(Float64, 5)
-#     for i ∈ eachindex(d_x_l)
-#         d_x_l[i] = op.dClosure[i-1]
-#         d_x_u[i] = -op.dClosure[length(d_x_u)-i]
-#     end
-#
-#     d_y_l = zeros(Float64, 6)
-#     d_y_u = zeros(Float64, 6)
-#     for i ∈ eachindex(d_y_l)
-#         d_y_l[i] = op.dClosure[i-1]
-#         d_y_u[i] = -op.dClosure[length(d_y_u)-i]
-#     end
-#
-#     function prod_matrix(x,y)
-#         G = zeros(Float64, length(x), length(y))
-#         for I ∈ CartesianIndices(G)
-#             G[I] = x[I[1]]*y[I[2]]
-#         end
-#
-#         return G
-#     end
-#
-#     g_x = [1,2,3,4.0,5]
-#     g_y = [5,4,3,2,1.0,11]
-#
-#     G_w = prod_matrix(d_x_l, g_y)
-#     G_e = prod_matrix(d_x_u, g_y)
-#     G_s = prod_matrix(g_x, d_y_l)
-#     G_n = prod_matrix(g_x, d_y_u)
-#
-#
-#     @test size(d_w*g_y) == (UnknownDim,6)
-#     @test size(d_e*g_y) == (UnknownDim,6)
-#     @test size(d_s*g_x) == (5,UnknownDim)
-#     @test size(d_n*g_x) == (5,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken d_w*g_y .≈ G_w
-#     @test_broken d_e*g_y .≈ G_e
-#     @test_broken d_s*g_x .≈ G_s
-#     @test_broken d_n*g_x .≈ G_n
-# end
-#
-# @testset "BoundaryQuadrature" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
-#
-#     H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
-#     H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
-#     H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
-#     H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-#
-#     function get_quadrature(N)
-#         qc = op.quadratureClosure
-#         q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
-#         @assert length(q) == N
-#         return q
-#     end
-#
-#     v_w = v[1,:]
-#     v_e = v[10,:]
-#     v_s = v[:,1]
-#     v_n = v[:,11]
-#
-#     q_x = spacing(g)[1].*get_quadrature(10)
-#     q_y = spacing(g)[2].*get_quadrature(11)
-#
-#     @test H_w isa TensorOperator{T,1} where T
-#
-#     @test domain_size(H_w, (3,)) == (3,)
-#     @test domain_size(H_n, (3,)) == (3,)
-#
-#     @test range_size(H_w, (3,)) == (3,)
-#     @test range_size(H_n, (3,)) == (3,)
-#
-#     @test size(H_w*v_w) == (11,)
-#     @test size(H_e*v_e) == (11,)
-#     @test size(H_s*v_s) == (10,)
-#     @test size(H_n*v_n) == (10,)
-#
-#     @test H_w*v_w .≈ q_y.*v_w
-#     @test H_e*v_e .≈ q_y.*v_e
-#     @test H_s*v_s .≈ q_x.*v_s
-#     @test H_n*v_n .≈ q_x.*v_n
-#
-#     @test H_w'*v_w == H_w'*v_w
-#     @test H_e'*v_e == H_e'*v_e
-#     @test H_s'*v_s == H_s'*v_s
-#     @test H_n'*v_n == H_n'*v_n
-# end
-
-end