--- a/DiffOps/src/laplace.jl	Wed Jun 26 21:22:36 2019 +0200
+++ b/DiffOps/src/laplace.jl	Thu Jun 27 12:49:58 2019 +0200
@@ -22,17 +22,13 @@
 	j = I[3-dim(d.bId)]
 	N_i = size(d.grid)[dim(d.bId)]
 
-	if region(d.bId) == Lower
-		# Note, closures are indexed by offset. Fix this D:<
-		return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
-	elseif region(d.bId) == Upper
-		return -d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-i]*v[j]
-	end
+    h_inv = d.grid.inverse_spacing[dim(d.bId)]
+    return apply_d(d.op, h_inv, v[j], N_i, i, region(d.bId))
 end
 
 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
     u = selectdim(v,3-dim(d.bId),I[1])
-    return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
+    return apply_d_T(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
 end
 
 
@@ -58,17 +54,12 @@
 	j = I[3-dim(e.bId)]
 	N_i = size(e.grid)[dim(e.bId)]
 
-	if region(e.bId) == Lower
-		# NOTE: closures are indexed by offset. Fix this D:<
-		return e.op.eClosure[i-1]*v[j]
-	elseif region(e.bId) == Upper
-		return e.op.eClosure[N_i-i]*v[j]
-	end
+    return apply_e(e.op, v[j], N_i, i, region(e.bId))
 end
 
 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
 	u = selectdim(v,3-dim(e.bId),I[1])
-	return apply_e(e.op, u, region(e.bId))
+	return apply_e_T(e.op, u, region(e.bId))
 end