--- a/diffOp.jl	Fri Jan 25 10:10:30 2019 +0100
+++ b/diffOp.jl	Fri Jan 25 13:40:15 2019 +0100
@@ -48,8 +48,9 @@
 end
 
 # u = L*v
+using UnsafeArrays
 @inline function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
-    u .= 0*u
+    u .= 0*u # Fix this?
     h = Grid.spacings(L.grid)
 
     li = LinearIndices(L.grid.numberOfPointsPerDim)
@@ -59,20 +60,21 @@
     # For each x
     temp = zeros(eltype(u), n_y)
     for i ∈ 1:n_x
-
-        @inbounds v_i = view(v, li[i,:])
-        @inbounds apply!(L.op, temp, v_i, h[2])
-
-        @inbounds u[li[i,:]] += temp
+        @inbounds indices = uview(li,i,:)
+        @inbounds apply!(L.op, temp, uview(v, indices), h[2])
+        for i ∈ eachindex(indices)
+            @inbounds u[indices[i]] += temp[i]
+        end
     end
 
     # For each y
     temp = zeros(eltype(u), n_x)
     for i ∈ 1:n_y
-        @inbounds v_i = view(v, li[:,i])
-        @inbounds apply!(L.op, temp, v_i, h[1])
-
-        @inbounds u[li[:,i]] += temp
+        @inbounds indices = uview(li,:,i)
+        @inbounds apply!(L.op, temp, uview(v, indices), h[1])
+        for i ∈ eachindex(indices)
+            @inbounds u[indices[i]] += temp[i]
+        end
     end
 
     u .= L.a*u

--- a/stencil.jl	Fri Jan 25 10:10:30 2019 +0100
+++ b/stencil.jl	Fri Jan 25 13:40:15 2019 +0100
@@ -17,18 +17,17 @@
 
 # Provides index into the Stencil based on offset for the root element
 function Base.getindex(s::Stencil, i::Int)
-    # TBD: Rearrange to mark with @boundscheck?
-    if s.range[1] <= i <= s.range[2]
-        return s.weights[1 + i - s.range[1]]
-    else
+    @boundscheck if i < s.range[1] || s.range[2] < i
         return eltype(s.weights)(0)
     end
+    return s.weights[1 + i - s.range[1]]
 end
 
-@inline function apply(s::Stencil, v::AbstractVector, i::Int)
+Base.@propagate_inbounds function apply(s::Stencil, v::AbstractVector, i::Int)
     w = zero(eltype(v))
     for j ∈ s.range[1]:s.range[2]
-        @inbounds w += s[j]*v[i+j] # TBD: Make this without boundschecks?
+        @inbounds weight = s[j]
+        w += weight*v[i+j]
     end
     return w
 end