--- a/README.md	Wed Sep 30 21:09:21 2020 +0200
+++ b/README.md	Wed Sep 30 21:09:35 2020 +0200
@@ -1,1 +1,15 @@
 # Sbplib
+
+## Running tests
+To run all tests simply run
+```
+(@v1.5) pkg> activate .
+(Sbplib) pkg> test
+```
+
+If you want to run tests from a specific file in `test/`, you can do
+```
+julia> using Pkg
+julia> Pkg.test(test_args=["testLazyTensors"])
+```
+This works by using the `@includetests` macro from the [TestSetExtensions](https://github.com/ssfrr/TestSetExtensions.jl) package. For more information, see their documentation.

--- a/src/Grids/EquidistantGrid.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/Grids/EquidistantGrid.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -66,6 +66,20 @@
     return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
 end
 
+"""
+    restrict(::EquidistantGrid, dim)
+
+Pick out given dimensions from the grid and return a grid for them
+"""
+function restrict(grid::EquidistantGrid, dim)
+    size = grid.size[dim]
+    limit_lower = grid.limit_lower[dim]
+    limit_upper = grid.limit_upper[dim]
+
+    return EquidistantGrid(size, limit_lower, limit_upper)
+end
+export restrict
+
 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
     @assert dim<=dimension(grid)
     @assert dim>0

--- a/src/LazyTensors/lazy_array.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/LazyTensors/lazy_array.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -17,6 +17,30 @@
 Base.getindex(lca::LazyConstantArray{T,D}, I::Vararg{Int,D}) where {T,D} = lca.val
 
 """
+    LazyFunctionArray{F<:Function,T, D} <: LazyArray{T,D}
+
+A lazy array where each element is defined by a function f(i,j,...)
+"""
+struct LazyFunctionArray{F<:Function,T, D} <: LazyArray{T,D}
+    f::F
+    size::NTuple{D,Int}
+end
+export LazyFunctionArray
+
+function LazyFunctionArray(f::F, size::NTuple{D,Int}) where {F<:Function,D}
+    T = typeof(f(ones(D)...))
+    return LazyFunctionArray{F,T,D}(f,size)
+end
+
+Base.size(lfa::LazyFunctionArray) = lfa.size
+
+function Base.getindex(lfa::LazyFunctionArray{F,T,D}, I::Vararg{Int,D}) where {F,T,D}
+    @boundscheck checkbounds(lfa, I...)
+    return lfa.f(I...)
+end
+
+
+"""
     LazyElementwiseOperation{T,D,Op} <: LazyArray{T,D}
 Struct allowing for lazy evaluation of elementwise operations on AbstractArrays.
 

--- a/src/LazyTensors/lazy_tensor_operations.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -11,12 +11,13 @@
     t::TensorMapping{T,R,D}
     o::AbstractArray{T,D}
 end
+# TODO: Do boundschecking on creation!
 export LazyTensorMappingApplication
 
 Base.:*(tm::TensorMapping{T,R,D}, o::AbstractArray{T,D}) where {T,R,D} = LazyTensorMappingApplication(tm,o)
 Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg{Index,R}) where {T,R,D} = apply(ta.t, ta.o, I...)
 Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg{Int,R}) where {T,R,D} = apply(ta.t, ta.o, Index{Unknown}.(I)...)
-Base.size(ta::LazyTensorMappingApplication{T,R,D}) where {T,R,D} = range_size(ta.t,size(ta.o))
+Base.size(ta::LazyTensorMappingApplication) = range_size(ta.t)
 # TODO: What else is needed to implement the AbstractArray interface?
 
 # # We need the associativity to be a→b→c = a→(b→c), which is the case for '→'
@@ -41,16 +42,15 @@
 export LazyTensorMappingTranspose
 
 # # TBD: Should this be implemented on a type by type basis or through a trait to provide earlier errors?
+# Jonatan 2020-09-25: Is the problem that you can take the transpose of any TensorMapping even if it doesn't implement `apply_transpose`?
 Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm)
 Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm
 
 apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Index,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...)
 apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmt.tm, v, I...)
 
-range_size(tmt::LazyTensorMappingTranspose{T,R,D}, d_size::NTuple{R,Integer}) where {T,R,D} = domain_size(tmt.tm, d_size)
-domain_size(tmt::LazyTensorMappingTranspose{T,R,D}, r_size::NTuple{D,Integer}) where {T,R,D} = range_size(tmt.tm, r_size)
-
-
+range_size(tmt::LazyTensorMappingTranspose) = domain_size(tmt.tm)
+domain_size(tmt::LazyTensorMappingTranspose) = range_size(tmt.tm)
 
 
 struct LazyTensorMappingBinaryOperation{Op,T,R,D,T1<:TensorMapping{T,R,D},T2<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R}
@@ -61,12 +61,13 @@
         return new{Op,T,R,D,T1,T2}(tm1,tm2)
     end
 end
+# TODO: Boundschecking in constructor.
 
 apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...)
 apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...)
 
-range_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}, domain_size::NTuple{D,Integer}) where {Op,T,R,D} = range_size(tmBinOp.tm1, domain_size)
-domain_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}, range_size::NTuple{R,Integer}) where {Op,T,R,D} = domain_size(tmBinOp.tm2, range_size)
+range_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}) where {Op,T,R,D} = range_size(tmBinOp.tm1)
+domain_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}) where {Op,T,R,D} = domain_size(tmBinOp.tm1)
 
 Base.:+(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:+,T,R,D}(tm1,tm2)
 Base.:-(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:-,T,R,D}(tm1,tm2)

--- a/src/LazyTensors/tensor_mapping.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/LazyTensors/tensor_mapping.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -6,11 +6,11 @@
 
     apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T}
 
-The size of range tensor should be dependent on the size of the domain tensor
-and the type should implement the methods
+The size of the range and domain that the operator works with should be returned by
+the functions
 
-    range_size(::TensorMapping{T,R,D}, domain_size::NTuple{D,Integer}) where {T,R,D}
-    domain_size(::TensorMapping{T,R,D}, range_size::NTuple{R,Integer}) where {T,R,D}
+    range_size(::TensorMapping)
+    domain_size(::TensorMapping)
 
 to allow querying for one or the other.
 
@@ -49,35 +49,19 @@
 export range_dim, domain_dim
 
 """
-    range_size(M::TensorMapping, domain_size)
+    range_size(M::TensorMapping)
 
-Return the resulting range size for the mapping applied to a given domain_size
+Return the range size for the mapping.
 """
 function range_size end
 
 """
-    domain_size(M::TensorMapping, range_size)
+    domain_size(M::TensorMapping)
 
-Return the resulting domain size for the mapping applied to a given range_size
+Return the domain size for the mapping.
 """
 function domain_size end
 
-"""
-    Dummy type for representing dimensions of tensormappings when domain_size is unknown
-"""
-struct UnknownDim end
-export range_size, domain_size, TensorMappingDim, UnknownDim
+export range_size, domain_size
 
 # TODO: Think about boundschecking!
-
-
-"""
-    TensorOperator{T,D}
-
-A `TensorMapping{T,D,D}` where the range and domain tensor have the same number of
-dimensions and the same size.
-"""
-abstract type TensorOperator{T,D} <: TensorMapping{T,D,D} end
-export TensorOperator
-domain_size(::TensorOperator{T,D}, range_size::NTuple{D,Integer}) where {T,D} = range_size
-range_size(::TensorOperator{T,D}, domain_size::NTuple{D,Integer}) where {T,D} = domain_size

--- a/src/SbpOperators/SbpOperators.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -2,6 +2,7 @@
 
 using Sbplib.RegionIndices
 using Sbplib.LazyTensors
+using Sbplib.Grids
 
 include("stencil.jl")
 include("constantstenciloperator.jl")
@@ -13,4 +14,5 @@
 include("quadrature/quadrature.jl")
 include("quadrature/inverse_diagonal_inner_product.jl")
 include("quadrature/inverse_quadrature.jl")
+
 end # module

--- a/src/SbpOperators/laplace/laplace.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/laplace/laplace.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -1,18 +1,27 @@
 export Laplace
 """
-    Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
+    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} <: TensorOperator{T,Dim}
+struct Laplace{Dim,T,N,M,K} <: TensorMapping{T,Dim,Dim}
     D2::NTuple{Dim,SecondDerivative{T,N,M,K}}
-    #TODO: Write a good constructor
 end
 
-LazyTensors.domain_size(L::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
+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{Index,Dim}) where {T,Dim}
     error("not implemented")
@@ -36,8 +45,6 @@
     return uᵢ
 end
 
-LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} = LazyTensors.apply(L, v, I...)
-
 # 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)

--- a/src/SbpOperators/laplace/secondderivative.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/laplace/secondderivative.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -4,16 +4,22 @@
 in 1D dimension
 """
 
-struct SecondDerivative{T,N,M,K} <: TensorOperator{T,1}
+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}}
     parity::Parity
-    #TODO: Write a nice constructor
+    size::NTuple{1,Int}
 end
 export SecondDerivative
 
-LazyTensors.domain_size(D2::SecondDerivative, range_size::NTuple{1,Integer}) = range_size
+function SecondDerivative(grid::EquidistantGrid{1}, innerStencil, closureStencils)
+    h_inv = grid.inverse_spacing[1]
+    return SecondDerivative(h_inv, innerStencil, closureStencils, even, size(grid))
+end
+
+LazyTensors.range_size(D2::SecondDerivative) = D2.size
+LazyTensors.domain_size(D2::SecondDerivative) = D2.size
 
 #TODO: The 1D tensor mappings should not have to dispatch on 1D tuples if we write LazyTensor.apply for vararg right?!?!
 #      Currently have to index the Tuple{Index} in each method in order to call the stencil methods which is ugly.
@@ -40,6 +46,4 @@
     return LazyTensors.apply(D2, v, I)
 end
 
-LazyTensors.apply_transpose(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index) where {T} = LazyTensors.apply(D2, v, I)
-
 closuresize(D2::SecondDerivative{T,N,M,K}) where {T<:Real,N,M,K} = M

--- a/src/SbpOperators/quadrature/diagonal_inner_product.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/quadrature/diagonal_inner_product.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -4,25 +4,30 @@
 
 Implements the diagnoal norm operator `H` of Dim dimension as a TensorMapping
 """
-struct DiagonalInnerProduct{T,M} <: TensorOperator{T,1}
-    h::T # The grid spacing could be included in the stencil already. Preferable?
-    closure::NTuple{M,T}
-    #TODO: Write a nice constructor
+struct DiagonalInnerProduct{T,M} <: TensorMapping{T,1,1}
+    h::T
+    quadratureClosure::NTuple{M,T}
+    size::Tuple{Int}
 end
 
-LazyTensors.domain_size(H::DiagonalInnerProduct, range_size::NTuple{1,Integer}) = range_size
+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) where T
     return @inbounds apply(H, v, I)
 end
 
 function LazyTensors.apply(H::DiagonalInnerProduct{T}, v::AbstractVector{T}, I::Index{Lower}) where T
-    return @inbounds H.h*H.closure[Int(I)]*v[Int(I)]
+    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.closure[N-Int(I)+1]v[Int(I)]
+    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

--- a/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/quadrature/inverse_diagonal_inner_product.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -1,22 +1,30 @@
 export InverseDiagonalInnerProduct, closuresize
 """
-    InverseDiagonalInnerProduct{Dim,T<:Real,M} <: TensorOperator{T,1}
+    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} <: TensorOperator{T,1}
-    h_inv::T # The reciprocl grid spacing could be included in the stencil already. Preferable?
-    closure::NTuple{M,T}
-    #TODO: Write a nice constructor
+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.closure[Int(I)]*v[Int(I)]
+    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.closure[N-Int(I)+1]v[Int(I)]
+    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

--- a/src/SbpOperators/quadrature/inverse_quadrature.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/quadrature/inverse_quadrature.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -2,14 +2,27 @@
 """
     InverseQuadrature{Dim,T<:Real,M,K} <: TensorMapping{T,Dim,Dim}
 
-Implements the inverse quadrature operator `Qi` of Dim dimension as a TensorOperator
+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} <: TensorOperator{T,Dim}
+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{Index,Dim}) where {T,Dim}

--- a/src/SbpOperators/quadrature/quadrature.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/src/SbpOperators/quadrature/quadrature.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -6,11 +6,21 @@
 The multi-dimensional tensor operator consists of a tuple of 1D DiagonalInnerProduct H
 tensor operators.
 """
-struct Quadrature{Dim,T<:Real,M} <: TensorOperator{T,Dim}
+struct Quadrature{Dim,T<:Real,M} <: TensorMapping{T,Dim,Dim}
     H::NTuple{Dim,DiagonalInnerProduct{T,M}}
 end
 
-LazyTensors.domain_size(Q::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
+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{Index,Dim}) where {T,Dim}
     error("not implemented")

--- a/test/testDiffOps.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/test/testDiffOps.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -6,268 +6,193 @@
 using Sbplib.LazyTensors
 
 @testset "DiffOps" begin
-
-@testset "Laplace2D" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    Lx = 3.5
-    Ly = 7.2
-    g = EquidistantGrid((42,41), (0.0, 0.0), (Lx,Ly))
-    L = Laplace(g, 1., op)
-    H = quadrature(L)
-
-    f0(x::Float64,y::Float64) = 2.
-    f1(x::Float64,y::Float64) = x+y
-    f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2
-    f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3
-    f4(x::Float64,y::Float64) = 1/24*x^4 + 1/24*y^4
-    f5(x::Float64,y::Float64) = sin(x) + cos(y)
-    f5ₓₓ(x::Float64,y::Float64) = -f5(x,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ₓₓ)
-
-    @test L isa TensorOperator{T,2} where T
-    @test L' isa TensorMapping{T,2,2} where T
-
-    # TODO: Should perhaps set tolerance level for isapporx instead?
-    #       Are these tolerance levels resonable or should tests be constructed
-    #       differently?
-    equalitytol = 0.5*1e-10
-    accuracytol = 0.5*1e-3
-    # 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 all(abs.(collect(L*v0)) .<= equalitytol)
-    @test all(abs.(collect(L*v1)) .<= equalitytol)
-    @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
-    @test all(abs.((collect(L*v3) - v1)) .<= equalitytol)
-    e4 = collect(L*v4) - v2
-    e5 = collect(L*v5) - v5ₓₓ
-    @test sum(collect(H*e4.^2)) <= accuracytol
-    @test sum(collect(H*e5.^2)) <= accuracytol
-end
-
-@testset "Quadrature" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    Lx = 2.3
-    Ly = 5.2
-    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-    H = Quadrature(op,g)
-    v = ones(Float64, size(g))
-
-    @test H isa TensorOperator{T,2} where T
-    @test H' isa TensorMapping{T,2,2} where T
-    @test sum(collect(H*v)) ≈ (Lx*Ly)
-    @test collect(H*v) == collect(H'*v)
-end
-
-@testset "InverseQuadrature" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    Lx = 7.3
-    Ly = 8.2
-    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
-    H = Quadrature(op,g)
-    Hinv = InverseQuadrature(op,g)
-    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-
-    @test Hinv isa TensorOperator{T,2} where T
-    @test Hinv' isa TensorMapping{T,2,2} where T
-    @test collect(Hinv*H*v)  ≈ v
-    @test collect(Hinv*v) == collect(Hinv'*v)
-end
-
-@testset "BoundaryValue" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
-
-    e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
-    e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
-    e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
-    e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
-
-    v = zeros(Float64, 4, 5)
-    v[:,5] = [1, 2, 3,4]
-    v[:,4] = [1, 2, 3,4]
-    v[:,3] = [4, 5, 6, 7]
-    v[:,2] = [7, 8, 9, 10]
-    v[:,1] = [10, 11, 12, 13]
-
-    @test e_w  isa TensorMapping{T,2,1} where T
-    @test e_w' isa TensorMapping{T,1,2} where T
-
-    @test domain_size(e_w, (3,2)) == (2,)
-    @test domain_size(e_e, (3,2)) == (2,)
-    @test domain_size(e_s, (3,2)) == (3,)
-    @test domain_size(e_n, (3,2)) == (3,)
-
-    @test size(e_w'*v) == (5,)
-    @test size(e_e'*v) == (5,)
-    @test size(e_s'*v) == (4,)
-    @test size(e_n'*v) == (4,)
-
-    @test collect(e_w'*v) == [10,7,4,1.0,1]
-    @test collect(e_e'*v) == [13,10,7,4,4.0]
-    @test collect(e_s'*v) == [10,11,12,13.0]
-    @test collect(e_n'*v) == [1,2,3,4.0]
-
-    g_x = [1,2,3,4.0]
-    g_y = [5,4,3,2,1.0]
-
-    G_w = zeros(Float64, (4,5))
-    G_w[1,:] = g_y
-
-    G_e = zeros(Float64, (4,5))
-    G_e[4,:] = g_y
-
-    G_s = zeros(Float64, (4,5))
-    G_s[:,1] = g_x
-
-    G_n = zeros(Float64, (4,5))
-    G_n[:,5] = g_x
-
-    @test size(e_w*g_y) == (UnknownDim,5)
-    @test size(e_e*g_y) == (UnknownDim,5)
-    @test size(e_s*g_x) == (4,UnknownDim)
-    @test size(e_n*g_x) == (4,UnknownDim)
+#
+# @testset "BoundaryValue" begin
+#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
+#
+#     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
+#     e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
+#     e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
+#     e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
+#
+#     v = zeros(Float64, 4, 5)
+#     v[:,5] = [1, 2, 3,4]
+#     v[:,4] = [1, 2, 3,4]
+#     v[:,3] = [4, 5, 6, 7]
+#     v[:,2] = [7, 8, 9, 10]
+#     v[:,1] = [10, 11, 12, 13]
+#
+#     @test e_w  isa TensorMapping{T,2,1} where T
+#     @test e_w' isa TensorMapping{T,1,2} where T
+#
+#     @test domain_size(e_w, (3,2)) == (2,)
+#     @test domain_size(e_e, (3,2)) == (2,)
+#     @test domain_size(e_s, (3,2)) == (3,)
+#     @test domain_size(e_n, (3,2)) == (3,)
+#
+#     @test size(e_w'*v) == (5,)
+#     @test size(e_e'*v) == (5,)
+#     @test size(e_s'*v) == (4,)
+#     @test size(e_n'*v) == (4,)
+#
+#     @test collect(e_w'*v) == [10,7,4,1.0,1]
+#     @test collect(e_e'*v) == [13,10,7,4,4.0]
+#     @test collect(e_s'*v) == [10,11,12,13.0]
+#     @test collect(e_n'*v) == [1,2,3,4.0]
+#
+#     g_x = [1,2,3,4.0]
+#     g_y = [5,4,3,2,1.0]
+#
+#     G_w = zeros(Float64, (4,5))
+#     G_w[1,:] = g_y
+#
+#     G_e = zeros(Float64, (4,5))
+#     G_e[4,:] = g_y
+#
+#     G_s = zeros(Float64, (4,5))
+#     G_s[:,1] = g_x
+#
+#     G_n = zeros(Float64, (4,5))
+#     G_n[:,5] = g_x
+#
+#     @test size(e_w*g_y) == (UnknownDim,5)
+#     @test size(e_e*g_y) == (UnknownDim,5)
+#     @test size(e_s*g_x) == (4,UnknownDim)
+#     @test size(e_n*g_x) == (4,UnknownDim)
+#
+#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
+#     @test_broken collect(e_w*g_y) == G_w
+#     @test_broken collect(e_e*g_y) == G_e
+#     @test_broken collect(e_s*g_x) == G_s
+#     @test_broken collect(e_n*g_x) == G_n
+# end
+#
+# @testset "NormalDerivative" begin
+#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     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 collect(d_w'*v) ≈ v∂x[1,:]
+#     @test collect(d_e'*v) ≈ v∂x[5,:]
+#     @test collect(d_s'*v) ≈ v∂y[:,1]
+#     @test collect(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 collect(d_w*g_y) ≈ G_w
+#     @test_broken collect(d_e*g_y) ≈ G_e
+#     @test_broken collect(d_s*g_x) ≈ G_s
+#     @test_broken collect(d_n*g_x) ≈ G_n
+# end
+#
+# @testset "BoundaryQuadrature" begin
+#     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+#     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 collect(H_w*v_w) ≈ q_y.*v_w
+#     @test collect(H_e*v_e) ≈ q_y.*v_e
+#     @test collect(H_s*v_s) ≈ q_x.*v_s
+#     @test collect(H_n*v_n) ≈ q_x.*v_n
+#
+#     @test collect(H_w'*v_w) == collect(H_w'*v_w)
+#     @test collect(H_e'*v_e) == collect(H_e'*v_e)
+#     @test collect(H_s'*v_s) == collect(H_s'*v_s)
+#     @test collect(H_n'*v_n) == collect(H_n'*v_n)
+# end
 
-    # These tests should be moved to where they are possible (i.e we know what the grid should be)
-    @test_broken collect(e_w*g_y) == G_w
-    @test_broken collect(e_e*g_y) == G_e
-    @test_broken collect(e_s*g_x) == G_s
-    @test_broken collect(e_n*g_x) == G_n
 end
-
-@testset "NormalDerivative" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    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 collect(d_w'*v) ≈ v∂x[1,:]
-    @test collect(d_e'*v) ≈ v∂x[5,:]
-    @test collect(d_s'*v) ≈ v∂y[:,1]
-    @test collect(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 collect(d_w*g_y) ≈ G_w
-    @test_broken collect(d_e*g_y) ≈ G_e
-    @test_broken collect(d_s*g_x) ≈ G_s
-    @test_broken collect(d_n*g_x) ≈ G_n
-end
-
-@testset "BoundaryQuadrature" begin
-    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    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 collect(H_w*v_w) ≈ q_y.*v_w
-    @test collect(H_e*v_e) ≈ q_y.*v_e
-    @test collect(H_s*v_s) ≈ q_x.*v_s
-    @test collect(H_n*v_n) ≈ q_x.*v_n
-
-    @test collect(H_w'*v_w) == collect(H_w'*v_w)
-    @test collect(H_e'*v_e) == collect(H_e'*v_e)
-    @test collect(H_s'*v_s) == collect(H_s'*v_s)
-    @test collect(H_n'*v_n) == collect(H_n'*v_n)
-end
-
-end
\ No newline at end of file

--- a/test/testGrids.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/test/testGrids.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -4,9 +4,23 @@
 @testset "Grids" begin
 
 @testset "EquidistantGrid" begin
-    @test EquidistantGrid(4,0,1) isa EquidistantGrid
-    @test dimension(EquidistantGrid(4,0,1)) == 1
-    @test EquidistantGrid(4,0,1) == EquidistantGrid((4,),(0,),(1,))
+    @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid
+    @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid
+    @test dimension(EquidistantGrid(4,0.0,1.0)) == 1
+    @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,))
+
+    g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))
+    @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0)
+    @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0)
+
+    g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
+    @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0)
+    @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0)
+    @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0)
+    @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0))
+    @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
+    @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0))
+    @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0))
 end
 
 end

--- a/test/testLazyTensors.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/test/testLazyTensors.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -12,20 +12,14 @@
     @test apply(DummyMapping{Int,2,3}(), zeros(Int, (0,0,0)),(Index{Unknown}(0),Index{Unknown}(0))) == :apply
 end
 
-@testset "Generic Operator methods" begin
-    struct DummyOperator{T,D} <: TensorOperator{T,D} end
-    @test range_size(DummyOperator{Int,2}(), (3,5)) == (3,5)
-    @test domain_size(DummyOperator{Float64, 3}(), (3,3,1)) == (3,3,1)
-end
-
 @testset "Mapping transpose" begin
     struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
 
     LazyTensors.apply(m::DummyMapping{T,R,D}, v, I::Vararg{Index{<:Region},R}) where {T,R,D} = :apply
     LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::Vararg{Index{<:Region},D}) where {T,R,D} = :apply_transpose
 
-    LazyTensors.range_size(m::DummyMapping{T,R,D}, domain_size::NTuple{D,Integer}) where {T,R,D} = :range_size
-    LazyTensors.domain_size(m::DummyMapping{T,R,D}, range_size::NTuple{R,Integer}) where {T,R,D} = :domain_size
+    LazyTensors.range_size(m::DummyMapping{T,R,D}) where {T,R,D} = :range_size
+    LazyTensors.domain_size(m::DummyMapping{T,R,D}) where {T,R,D} = :domain_size
 
     m = DummyMapping{Float64,2,3}()
     I = Index{Unknown}(0)
@@ -35,19 +29,21 @@
     @test apply(m'',zeros(Float64,(0,0,0)), I, I) == :apply
     @test apply_transpose(m', zeros(Float64,(0,0,0)), I, I) == :apply
 
-    @test range_size(m', (0,0)) == :domain_size
-    @test domain_size(m', (0,0,0)) == :range_size
+    @test range_size(m') == :domain_size
+    @test domain_size(m') == :range_size
 end
 
 @testset "TensorApplication" begin
-    struct DummyMapping{T,R,D} <: TensorMapping{T,R,D} end
+    struct SizeDoublingMapping{T,R,D} <: TensorMapping{T,R,D}
+        domain_size::NTuple{D,Int}
+    end
 
-    LazyTensors.apply(m::DummyMapping{T,R,D}, v, i::Vararg{Index{<:Region},R}) where {T,R,D} = (:apply,v,i)
-    LazyTensors.range_size(m::DummyMapping{T,R,D}, domain_size::NTuple{D,Integer}) where {T,R,D} = 2 .* domain_size
-    LazyTensors.domain_size(m::DummyMapping{T,R,D}, range_size::NTuple{R,Integer}) where {T,R,D} = range_size.÷2
+    LazyTensors.apply(m::SizeDoublingMapping{T,R,D}, v, i::Vararg{Index{<:Region},R}) where {T,R,D} = (:apply,v,i)
+    LazyTensors.range_size(m::SizeDoublingMapping) = 2 .* m.domain_size
+    LazyTensors.domain_size(m::SizeDoublingMapping) = m.domain_size
 
 
-    m = DummyMapping{Int, 1, 1}()
+    m = SizeDoublingMapping{Int, 1, 1}((3,))
     v = [0,1,2]
     @test m*v isa AbstractVector{Int}
     @test size(m*v) == 2 .*size(v)
@@ -63,28 +59,31 @@
     @test_broken BoundsError == (m*m*v)[0]
     @test_broken BoundsError == (m*m*v)[7]
 
-    m = DummyMapping{Int, 2, 1}()
+    m = SizeDoublingMapping{Int, 2, 1}((3,))
     @test_throws MethodError m*ones(Int,2,2)
     @test_throws MethodError m*m*v
 
-    m = DummyMapping{Float64, 2, 2}()
+    m = SizeDoublingMapping{Float64, 2, 2}((3,3))
     v = ones(3,3)
     I = (Index{Lower}(1),Index{Interior}(2));
     @test size(m*v) == 2 .*size(v)
     @test (m*v)[I] == (:apply,v,I)
 
-    struct ScalingOperator{T,D} <: TensorOperator{T,D}
+    struct ScalingOperator{T,D} <: TensorMapping{T,D,D}
         λ::T
+        size::NTuple{D,Int}
     end
 
     LazyTensors.apply(m::ScalingOperator{T,D}, v, I::Vararg{Index,D}) where {T,D} = m.λ*v[I]
+    LazyTensors.range_size(m::ScalingOperator) = m.size
+    LazyTensors.domain_size(m::ScalingOperator) = m.size
 
-    m = ScalingOperator{Int,1}(2)
+    m = ScalingOperator{Int,1}(2,(3,))
     v = [1,2,3]
     @test m*v isa AbstractVector
     @test m*v == [2,4,6]
 
-    m = ScalingOperator{Int,2}(2)
+    m = ScalingOperator{Int,2}(2,(2,2))
     v = [[1 2];[3 4]]
     @test m*v == [[2 4];[6 8]]
     I = (Index{Upper}(2),Index{Lower}(1))
@@ -94,14 +93,16 @@
 @testset "TensorMapping binary operations" begin
     struct ScalarMapping{T,R,D} <: TensorMapping{T,R,D}
         λ::T
+        range_size::NTuple{R,Int}
+        domain_size::NTuple{D,Int}
     end
 
     LazyTensors.apply(m::ScalarMapping{T,R,D}, v, I::Vararg{Index{<:Region}}) where {T,R,D} = m.λ*v[I...]
-    LazyTensors.range_size(m::ScalarMapping, domain_size) = domain_size
-    LazyTensors.domain_size(m::ScalarMapping, range_sizes) = range_sizes
+    LazyTensors.range_size(m::ScalarMapping) = m.domain_size
+    LazyTensors.domain_size(m::ScalarMapping) = m.range_size
 
-    A = ScalarMapping{Float64,1,1}(2.0)
-    B = ScalarMapping{Float64,1,1}(3.0)
+    A = ScalarMapping{Float64,1,1}(2.0, (3,), (3,))
+    B = ScalarMapping{Float64,1,1}(3.0, (3,), (3,))
 
     v = [1.1,1.2,1.3]
     for i ∈ eachindex(v)
@@ -112,8 +113,8 @@
         @test ((A-B)*v)[i] == 2*v[i] - 3*v[i]
     end
 
-    @test range_size(A+B, (3,)) == range_size(A, (3,)) == range_size(B,(3,))
-    @test domain_size(A+B, (3,)) == domain_size(A, (3,)) == domain_size(B,(3,))
+    @test range_size(A+B) == range_size(A) == range_size(B)
+    @test domain_size(A+B) == domain_size(A) == domain_size(B)
 end
 
 @testset "LazyArray" begin
@@ -192,4 +193,23 @@
     @test_throws DimensionMismatch v1 + v2
 end
 
-end
\ No newline at end of file
+@testset "LazyFunctionArray" begin
+    @test LazyFunctionArray(i->i^2, (3,)) == [1,4,9]
+    @test LazyFunctionArray((i,j)->i*j, (3,2)) == [
+        1 2;
+        2 4;
+        3 6;
+    ]
+
+    @test size(LazyFunctionArray(i->i^2, (3,))) == (3,)
+    @test size(LazyFunctionArray((i,j)->i*j, (3,2))) == (3,2)
+
+    @inferred LazyFunctionArray(i->i^2, (3,))[2]
+
+    @test_throws BoundsError LazyFunctionArray(i->i^2, (3,))[4]
+    @test_throws BoundsError LazyFunctionArray((i,j)->i*j, (3,2))[4,2]
+    @test_throws BoundsError LazyFunctionArray((i,j)->i*j, (3,2))[2,3]
+
+end
+
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/testRegionIndices.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -0,0 +1,6 @@
+using Sbplib.RegionIndices
+using Test
+
+@testset "RegionIndices" begin
+	@test_broken false
+end

--- a/test/testRegionInidices.jl	Wed Sep 30 21:09:21 2020 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-using Sbplib.RegionIndices
-using Test
-
-@testset "RegionIndices" begin
-	@test_broken false
-end

--- a/test/testSbpOperators.jl	Wed Sep 30 21:09:21 2020 +0200
+++ b/test/testSbpOperators.jl	Wed Sep 30 21:09:35 2020 +0200
@@ -33,7 +33,7 @@
     g = EquidistantGrid((101,), (0.0,), (L,))
     h_inv = inverse_spacing(g)
     h = 1/h_inv[1];
-    Dₓₓ = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity)
+    Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils)
 
     f0(x::Float64) = 1.
     f1(x::Float64) = x
@@ -50,7 +50,7 @@
     v4 = evalOn(g,f4)
     v5 = evalOn(g,f5)
 
-    @test Dₓₓ isa TensorOperator{T,1} where T
+    @test Dₓₓ isa TensorMapping{T,1,1} where T
     @test Dₓₓ' isa TensorMapping{T,1,1} where T
 
     # TODO: Should perhaps set tolerance level for isapporx instead?
@@ -77,12 +77,8 @@
     Lx = 1.5
     Ly = 3.2
     g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
+    L = Laplace(g, op.innerStencil, op.closureStencils)
 
-    h_inv = inverse_spacing(g)
-    h = spacing(g)
-    D_xx = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity)
-    D_yy = SecondDerivative(h_inv[2],op.innerStencil,op.closureStencils,op.parity)
-    L = Laplace((D_xx,D_yy))
 
     f0(x::Float64,y::Float64) = 2.
     f1(x::Float64,y::Float64) = x+y
@@ -100,7 +96,7 @@
     v5 = evalOn(g,f5)
     v5ₓₓ = evalOn(g,f5ₓₓ)
 
-    @test L isa TensorOperator{T,2} where T
+    @test L isa TensorMapping{T,2,2} where T
     @test L' isa TensorMapping{T,2,2} where T
 
     # TODO: Should perhaps set tolerance level for isapporx instead?
@@ -118,6 +114,8 @@
     @test all(abs.((collect(L*v3) - v1)) .<= equalitytol)
     e4 = collect(L*v4) - v2
     e5 = collect(L*v5) - v5ₓₓ
+
+    h = spacing(g)
     @test sqrt(prod(h)*sum(collect(e4.^2))) <= accuracytol
     @test sqrt(prod(h)*sum(collect(e5.^2))) <= accuracytol
 end
@@ -126,11 +124,10 @@
     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
     L = 2.3
     g = EquidistantGrid((77,), (0.0,), (L,))
-    h = spacing(g)
-    H = DiagonalInnerProduct(h[1],op.quadratureClosure)
+    H = DiagonalInnerProduct(g,op.quadratureClosure)
     v = ones(Float64, size(g))
 
-    @test H isa TensorOperator{T,1} where T
+    @test H isa TensorMapping{T,1,1} where T
     @test H' isa TensorMapping{T,1,1} where T
     @test sum(collect(H*v)) ≈ L
     @test collect(H*v) == collect(H'*v)
@@ -142,14 +139,11 @@
     Ly = 5.2
     g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
 
-    h = spacing(g)
-    Hx = DiagonalInnerProduct(h[1],op.quadratureClosure);
-    Hy = DiagonalInnerProduct(h[2],op.quadratureClosure);
-    Q = Quadrature((Hx,Hy))
+    Q = Quadrature(g, op.quadratureClosure)
 
     v = ones(Float64, size(g))
 
-    @test Q isa TensorOperator{T,2} where T
+    @test Q isa TensorMapping{T,2,2} where T
     @test Q' isa TensorMapping{T,2,2} where T
     @test sum(collect(Q*v)) ≈ (Lx*Ly)
     @test collect(Q*v) == collect(Q'*v)
@@ -159,14 +153,11 @@
     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
     L = 2.3
     g = EquidistantGrid((77,), (0.0,), (L,))
-    h = spacing(g)
-    H = DiagonalInnerProduct(h[1],op.quadratureClosure)
-
-    h_i = inverse_spacing(g)
-    Hi = InverseDiagonalInnerProduct(h_i[1],1 ./ op.quadratureClosure)
+    H = DiagonalInnerProduct(g, op.quadratureClosure)
+    Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure)
     v = evalOn(g, x->sin(x))
 
-    @test Hi isa TensorOperator{T,1} where T
+    @test Hi isa TensorMapping{T,1,1} where T
     @test Hi' isa TensorMapping{T,1,1} where T
     @test collect(Hi*H*v)  ≈ v
     @test collect(Hi*v) == collect(Hi'*v)
@@ -178,20 +169,13 @@
     Ly = 8.2
     g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
 
-    h = spacing(g)
-    Hx = DiagonalInnerProduct(h[1], op.quadratureClosure);
-    Hy = DiagonalInnerProduct(h[2], op.quadratureClosure);
-    Q = Quadrature((Hx,Hy))
-
-    hi = inverse_spacing(g)
-    Hix = InverseDiagonalInnerProduct(hi[1], 1 ./ op.quadratureClosure);
-    Hiy = InverseDiagonalInnerProduct(hi[2], 1 ./ op.quadratureClosure);
-    Qinv = InverseQuadrature((Hix,Hiy))
+    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 TensorOperator{T,2} where T
+    @test Qinv isa TensorMapping{T,2,2} where T
     @test Qinv' isa TensorMapping{T,2,2} where T
-    @test collect(Qinv*Q*v)  ≈ v
+    @test collect(Qinv*(Q*v)) ≈ v
     @test collect(Qinv*v) == collect(Qinv'*v)
 end
 #