--- a/DiffOps/src/DiffOps.jl	Thu Jan 09 10:53:03 2020 +0100
+++ b/DiffOps/src/DiffOps.jl	Thu Jan 09 10:54:24 2020 +0100
@@ -98,7 +98,5 @@
 
 
 include("laplace.jl")
-export Laplace
-
 
 end # module

--- a/DiffOps/src/laplace.jl	Thu Jan 09 10:53:03 2020 +0100
+++ b/DiffOps/src/laplace.jl	Thu Jan 09 10:54:24 2020 +0100
@@ -1,36 +1,33 @@
-struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
+struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
     grid::EquidistantGrid{Dim,T}
-    a::T
-    op::D2{Float64,N,M,K}
+    a::T # TODO: Better name?
+    op::D2{T,N,M,K}
 end
+export Laplace
 
-function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
+LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where Dim = size(L.grid)
+
+function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
     error("not implemented")
 end
 
 # u = L*v
-function apply(L::Laplace{1}, v::AbstractVector, i::Int)
-    uᵢ = L.a * SbpOperators.apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], v, i)
+function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
+    uᵢ = L.a*apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], v, I[1])
     return uᵢ
 end
 
-@inline function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
+
+@inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = L.a*SbpOperators.apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], vx , I[1])
+    @inbounds uᵢ = L.a*apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], vx , I[1])
     # 2nd y-derivative
     @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += L.a*SbpOperators.apply_2nd_derivative(L.op, inverse_spacing(L.grid)[2], vy, I[2])
-    # NOTE: the package qualifier 'SbpOperators' can problably be removed once all "applying" objects use LazyTensors
+    @inbounds uᵢ += L.a*apply_2nd_derivative(L.op, inverse_spacing(L.grid)[2], vy, I[2])
     return uᵢ
 end
 
-# Slow but maybe convenient?
-function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
-    I = Index{Unknown}.(Tuple(i))
-    apply(L, v, I)
-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)

--- a/DiffOps/test/runtests.jl	Thu Jan 09 10:53:03 2020 +0100
+++ b/DiffOps/test/runtests.jl	Thu Jan 09 10:54:24 2020 +0100
@@ -5,6 +5,44 @@
 using RegionIndices
 using LazyTensors
 
+@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((21,42), (0.0, 0.0), (Lx,Ly))
+    L = Laplace(g, 1., op)
+
+    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) = x^5 + 2*y^3 + 3/2*x^2 + y + 7
+    f5ₓₓ(x::Float64,y::Float64) = 20*x^3 + 12*y + 3
+
+    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?
+    equalitytol = 0.5*1e-11
+    accuracytol = 1e-4
+    @test all(collect(L*v0) .<= equalitytol)
+    @test all(collect(L*v1) .<= equalitytol)
+    @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
+    @test all((collect(L*v3) - v1) .<= equalitytol)
+    @show maximum(abs.(collect(L*v4) - v2))
+    @show maximum(abs.(collect(L*v5) - v5ₓₓ))
+    @test_broken all((collect(L*v4) - v2) .<= accuracytol) #TODO: Should be equality?
+    @test_broken all((collect(L*v5) - v5ₓₓ) .<= accuracytol) #TODO: This is just wrong
+end
+
 @testset "Quadrature" begin
     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
     Lx = 2.3
@@ -13,7 +51,7 @@
     H = Quadrature(op,g)
     v = ones(Float64, size(g))
 
-    @test H isa TensorMapping{T,2,2} where T
+    @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)
@@ -28,7 +66,7 @@
     Hinv = InverseQuadrature(op,g)
     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
 
-    @test Hinv isa TensorMapping{T,2,2} where T
+    @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)
@@ -141,7 +179,7 @@
         d_y_u[i] = -op.dClosure[length(d_y_u)-i]
     end
 
-    function ❓(x,y)
+    function prod_matrix(x,y)
         G = zeros(Float64, length(x), length(y))
         for I ∈ CartesianIndices(G)
             G[I] = x[I[1]]*y[I[2]]
@@ -153,10 +191,10 @@
     g_x = [1,2,3,4.0,5]
     g_y = [5,4,3,2,1.0,11]
 
-    G_w = ❓(d_x_l, g_y)
-    G_e = ❓(d_x_u, g_y)
-    G_s = ❓(g_x, d_y_l)
-    G_n = ❓(g_x, d_y_u)
+    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) == (5,6)