--- a/DiffOps/Manifest.toml	Wed Jun 26 17:53:15 2019 +0200
+++ b/DiffOps/Manifest.toml	Wed Jun 26 17:54:32 2019 +0200
@@ -17,6 +17,11 @@
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[LazyTensors]]
+path = "../LazyTensors"
+uuid = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
+version = "0.1.0"
+
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 

--- a/DiffOps/Project.toml	Wed Jun 26 17:53:15 2019 +0200
+++ b/DiffOps/Project.toml	Wed Jun 26 17:54:32 2019 +0200
@@ -5,6 +5,7 @@
 
 [deps]
 Grids = "960fdf28-97ed-11e9-2a74-bd90bf2fab5a"
+LazyTensors = "62fbed2c-918d-11e9-279b-eb3a325b37d3"
 RegionIndices = "5d527584-97f1-11e9-084c-4540c7ecf219"
 SbpOperators = "204357d8-97fd-11e9-05e7-010897a14cd0"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"

--- a/DiffOps/src/DiffOps.jl	Wed Jun 26 17:53:15 2019 +0200
+++ b/DiffOps/src/DiffOps.jl	Wed Jun 26 17:54:32 2019 +0200
@@ -3,6 +3,7 @@
 using RegionIndices
 using SbpOperators
 using Grids
+using LazyTensors
 
 """
     DiffOp

--- a/DiffOps/src/laplace.jl	Wed Jun 26 17:53:15 2019 +0200
+++ b/DiffOps/src/laplace.jl	Wed Jun 26 17:54:32 2019 +0200
@@ -1,17 +1,18 @@
-struct NormalDerivative{N,M,K}
-	op::D2{Float64,N,M,K}
+"""
+    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
 	bId::CartesianBoundary
 end
-
-function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer)
-	u = selectdim(v,3-dim(d.bId),I)
-	return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
-end
+export NormalDerivative
 
 # Not correct abstraction level
 # TODO: Not type stable D:<
-function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer})
+function apply(d::NormalDerivative, v::AbstractArray, I::CartesianIndex{2})
 	i = I[dim(d.bId)]
 	j = I[3-dim(d.bId)]
 	N_i = d.grid.size[dim(d.bId)]
@@ -30,13 +31,30 @@
 	end
 end
 
-struct BoundaryValue{N,M,K}
-	op::D2{Float64,N,M,K}
+function apply_transpose(d::NormalDerivative, v::AbstractArray, I::CartesianIndex{1})
+    u = selectdim(v,3-dim(d.bId),I)
+    return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
+end
+
+
+"""
+    BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+
+Implements the boundary operator `e` as a TensorMapping
+"""
+struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
+	op::D2{T,N,M,K}
 	grid::EquidistantGrid
 	bId::CartesianBoundary
 end
+export BoundaryValue
 
-function apply(e::BoundaryValue, v::AbstractArray, I::Tuple{Integer,Integer})
+# 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?
+LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
+LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],);
+
+function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::CartesianIndex{2})
 	i = I[dim(e.bId)]
 	j = I[3-dim(e.bId)]
 	N_i = e.grid.size[dim(e.bId)]
@@ -55,11 +73,13 @@
 	end
 end
 
-function apply_transpose(e::BoundaryValue, v::AbstractArray, I::Integer)
+function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::CartesianIndex{1})
 	u = selectdim(v,3-dim(e.bId),I)
 	return apply_e(e.op, u, region(e.bId))
 end
 
+
+
 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
     grid::EquidistantGrid{Dim,T}
     a::T

--- a/DiffOps/test/runtests.jl	Wed Jun 26 17:53:15 2019 +0200
+++ b/DiffOps/test/runtests.jl	Wed Jun 26 17:54:32 2019 +0200
@@ -1,4 +1,37 @@
+using Test
 using DiffOps
-using Test
+using Grids
+using SbpOperators
+using RegionIndices
+using LazyTensors
 
 @test_broken false
+
+@testset "BoundaryValue" begin
+    op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
+    g = EquidistantGrid((3,3), (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 = [
+        1 2 3;
+        4 5 6;
+        7 8 9.0;
+        10 11 12;
+    ]
+
+    @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 collect(e_w'*v) == [1,4,7.0]
+
+end