--- a/DiffOps/src/laplace.jl	Wed Jun 26 20:23:23 2019 +0200
+++ b/DiffOps/src/laplace.jl	Wed Jun 26 20:24:21 2019 +0200
@@ -12,7 +12,7 @@
 
 # Not correct abstraction level
 # TODO: Not type stable D:<
-function apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
+function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
 	i = I[dim(d.bId)]
 	j = I[3-dim(d.bId)]
 	N_i = d.grid.size[dim(d.bId)]
@@ -31,7 +31,7 @@
 	end
 end
 
-function apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
+function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
     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
@@ -57,19 +57,13 @@
 function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int})
 	i = I[dim(e.bId)]
 	j = I[3-dim(e.bId)]
-	N_i = e.grid.size[dim(e.bId)]
-
-	r = getregion(i, closureSize(e.op), N_i)
+	N_i = size(e.grid)[dim(e.bId)]
 
-	if r != region(e.bId)
-		return 0
-	end
-
-	if r == Lower
-		# Note, closures are indexed by offset. Fix this D:<
+	if region(e.bId) == Lower
+		# NOTE: closures are indexed by offset. Fix this D:<
 		return e.op.eClosure[i-1]*v[j]
-	elseif r == Upper
-		return e.op.eClosure[N_i-j]*v[j]
+	elseif region(e.bId) == Upper
+		return e.op.eClosure[N_i-i]*v[j]
 	end
 end
 

--- a/DiffOps/test/runtests.jl	Wed Jun 26 20:23:23 2019 +0200
+++ b/DiffOps/test/runtests.jl	Wed Jun 26 20:24:21 2019 +0200
@@ -41,4 +41,24 @@
     @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 collect(e_w*g_y) == G_w
+    @test collect(e_e*g_y) == G_e
+    @test collect(e_s*g_x) == G_s
+    @test collect(e_n*g_x) == G_n
+
 end