--- a/diffOp.jl	Mon Feb 04 09:11:53 2019 +0100
+++ b/diffOp.jl	Mon Feb 04 09:13:48 2019 +0100
@@ -37,7 +37,6 @@
     for i ∈ 1:Grid.numberOfPoints(D.grid)
         u[i] = apply(D, v, i)
     end
-
     return nothing
 end
 
@@ -63,6 +62,38 @@
 using UnsafeArrays
 
 # u = L*v
+function apply(L::Laplace{1}, v::AbstractVector, i::Int)
+    h = Grid.spacings(L.grid)[1]
+    uᵢ = L.a * apply(L.op, h, v, i)
+    return uᵢ
+end
+
+function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
+    lowerind, upperind, interiorind = stencilindices(L.grid)
+    for si ∈ lowerind
+       u[si.globalindex] = apply(L, v, si)
+    end
+    for si ∈ upperind
+       u[si.globalindex] = apply(L, v, si)
+    end
+    for si ∈ interiorind
+        u[si.globalindex] = apply(L, v, si)
+    end
+    return nothing
+end
+
+@inline function apply(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
+    h = Grid.spacings(L.grid)
+    li = LinearIndices(L.grid.numberOfPointsPerDim)
+    # 2nd x-derivative
+    @inbounds vx = uview(v, uview(li,:,si.gridindex[2]))
+    uᵢ  = apply(L.op, h[1], vx , si.gridindex[1], si)
+    # 2nd y-derivative
+    @inbounds vy = uview(v, uview(li,si.gridindex[1],:))
+    uᵢ += apply(L.op, h[2], vy, si.gridindex[2], si)
+    return uᵢ
+end
+
 function apply(L::Laplace{2}, v::AbstractVector, i::Int)
     h = Grid.spacings(L.grid)