--- a/EquidistantGrid.jl	Mon Feb 04 09:11:53 2019 +0100
+++ b/EquidistantGrid.jl	Mon Feb 04 09:13:48 2019 +0100
@@ -80,7 +80,7 @@
     points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim])
 end
 
-# TODO: Move to own plotting module.
+# # TODO: Move to own plotting module.
 using PyPlot, PyCall
 
 function plotgridfunction(grid::EquidistantGrid, gridfunction)

--- a/StencilIndex.jl	Mon Feb 04 09:11:53 2019 +0100
+++ b/StencilIndex.jl	Mon Feb 04 09:13:48 2019 +0100
@@ -1,21 +1,25 @@
 abstract type StencilIndex end
 
 function Base.getindex(si::StencilIndex, i::Int)
-    return si.index[i]
+    return si.gridindex[i]
 end
 
 struct LowerClosureIndex <: StencilIndex
-    index::CartesianIndex
+    globalindex::Integer
+    gridindex::CartesianIndex
 end
 
-struct UpperClosureIndex  <: StencilIndex
-    index::CartesianIndex
+struct UpperClosureIndex <: StencilIndex
+    globalindex::Integer
+    gridindex::CartesianIndex
 end
 
-struct InteriorIndex  <: StencilIndex
-    index::CartesianIndex
+struct InteriorIndex <: StencilIndex
+    globalindex::Integer
+    gridindex::CartesianIndex
 end
 
+# TODO: The design of StencilIndex is wrong. Use Jonatans design instead.
 # TODO: This should take a Stencil or DiffOp so that we can extract all the
 # indices in the closures.
 # TODO: Where to place this function?
@@ -27,16 +31,14 @@
     islower = x -> (x == 1)
     isupper = x -> (x in grid.numberOfPointsPerDim)
     ci = CartesianIndices(grid.numberOfPointsPerDim)
-    for i ∈ ci
-        I = Tuple(i)
+    for i ∈ 1:Grid.numberOfPoints(grid)
+        I = Tuple(ci[i])
         if any(islower, I)
-            push!(lowerclosure, LowerClosureIndex(i))
-            # TODO: Corner points should be in both Lower and Upper?
-            #       Should they have a separate type?
+            push!(lowerclosure, LowerClosureIndex(i,ci[i]))
         elseif any(isupper, I)
-            push!(upperclosure, UpperClosureIndex(i))
+            push!(upperclosure, UpperClosureIndex(i,ci[i]))
         else
-            push!(interior, InteriorIndex(i))
+            push!(interior, InteriorIndex(i,ci[i]))
         end
     end
     return lowerclosure, upperclosure, interior

--- 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)
 

--- a/sbpD2.jl	Mon Feb 04 09:11:53 2019 +0100
+++ b/sbpD2.jl	Mon Feb 04 09:13:48 2019 +0100
@@ -15,15 +15,16 @@
     return uᵢ
 end
 
-Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::InteriorIndex)
-    return apply(op.innerStencil, v, i)/h^2
+Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::InteriorIndex)
+    return apply(op.innerStencil, v,  i)/h^2
 end
 
-Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::LowerClosureIndex)
+Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::LowerClosureIndex)
     return apply(op.closureStencils[i], v, i)/h^2
 end
 
-Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::UpperClosureIndex)
+Base.@propagate_inbounds function apply(op::ConstantStencilOperator, h::Real, v::AbstractVector, i::Int, ::UpperClosureIndex)
+    N = length(v)
     return Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2 #TODO: Write an applybackwards instead?
 end
 

--- a/stencil.jl	Mon Feb 04 09:11:53 2019 +0100
+++ b/stencil.jl	Mon Feb 04 09:13:48 2019 +0100
@@ -19,6 +19,7 @@
 
 Base.@propagate_inbounds function apply(s::Stencil, v::AbstractVector, i::Int)
     w = zero(eltype(v))
+    @show s.range[1]:s.range[2]
     for j ∈ s.range[1]:s.range[2]
         @inbounds weight = s[j]
         w += weight*v[i+j]