--- a/diffOp.jl	Mon Feb 04 09:13:48 2019 +0100
+++ b/diffOp.jl	Mon Feb 04 16:09:07 2019 +0100
@@ -69,28 +69,55 @@
 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)
+    lowerind, upperind, interiorind = stencilindices(L)
+    # First x!
+    for si ∈ lowerind[1]
+       u[si.globalindex] = applyx(L, v, si)
+    end
+    for si ∈ upperind[1]
+       u[si.globalindex] = applyx(L, v, si)
     end
-    for si ∈ upperind
-       u[si.globalindex] = apply(L, v, si)
+    for si ∈ interiorind[1]
+        u[si.globalindex] = applyx(L, v, si)
     end
-    for si ∈ interiorind
-        u[si.globalindex] = apply(L, v, si)
+    # NOW y!
+    for si ∈ lowerind[2]
+       u[si.globalindex] += applyy(L, v, si)
+    end
+    for si ∈ upperind[2]
+       u[si.globalindex] += applyy(L, v, si)
+    end
+    for si ∈ interiorind[2]
+        u[si.globalindex] += applyy(L, v, si)
     end
     return nothing
 end
 
-@inline function apply(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
+@inline function applyx(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)
+    @inbounds vx = uview(v, uview(li,:,si.localindex[2]))
+    uᵢ  = apply(L.op, h[1], vx , si.localindex[1], si)
+    return uᵢ
+end
+
+@inline function applyy(L::Laplace{2}, v::AbstractVector, si::StencilIndex)
+    h = Grid.spacings(L.grid)
+    li = LinearIndices(L.grid.numberOfPointsPerDim)
     # 2nd y-derivative
-    @inbounds vy = uview(v, uview(li,si.gridindex[1],:))
-    uᵢ += apply(L.op, h[2], vy, si.gridindex[2], si)
+    # @show typeof(si.localindex)
+    # @show si.localindex
+
+    # @show typeof(si.localindex[1])
+    # @show typeof(si.localindex[2])
+
+    # @show typeof(uview(li, si.localindex[1],:))
+    # @show uview(li,si.localindex[1],:)
+
+
+    @inbounds vy = uview(v, uview(li,si.localindex[1],:))
+    uᵢ = apply(L.op, h[2], vy, si.localindex[2], si)
     return uᵢ
 end