diff -r 3d8bfb695497 -r c061d1bddba5 diffOp.jl
--- a/diffOp.jl	Thu Jan 10 16:48:49 2019 +0100
+++ b/diffOp.jl	Thu Jan 10 16:49:20 2019 +0100
@@ -1,6 +1,6 @@
 abstract type DiffOp end
 
-function apply(D::DiffOp, v::AbstractVector)
+function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
     error("not implemented")
 end
 
@@ -35,17 +35,6 @@
 
     h = scaling(L.grid)
 
-    for i ∈ 1:N
-        u[i] = apply(L.op.closureStencils[i], v, i)/h^2
-    end
-
-    for i ∈ N+1:M-N
-        u[i] = apply(L.op.innerStencil, i)/h^2
-    end
-
-    for i ∈ M:-1:M-N+1
-        u[i] = apply(flip(L.op.closureStencils[M-i+1]), v, i)/h^2
-    end
-
+    apply!(L.op, u, v, grid.spacings(L.grid)[1], 1, L.grid.numberOfPointsPerDim, stride=1)
     return nothing
 end

diff -r 3d8bfb695497 -r c061d1bddba5 sbpD2.jl
--- a/sbpD2.jl	Thu Jan 10 16:48:49 2019 +0100
+++ b/sbpD2.jl	Thu Jan 10 16:49:20 2019 +0100
@@ -1,9 +1,33 @@
-struct D2{T}
+abstract type ConstantStencilOperator end
+
+function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real, start::Int, N::Int; stride::Int=1)
+    cSize = closureSize(op)
+
+    for i ∈ range(start; length=cSize, step=stride)
+        u[i] = apply!(op.closureStencils[i], v, i; stride=stride)/h^2
+    end
+
+    innerStart = start + cSize*stride
+    innerEnd = N - cSize*stride-1
+    for i ∈ range(innerStart, stop=innerEnd, step=stride)
+        u[i] = apply(op.innerStencil, v, i; stride=stride)/h^2
+    end
+
+    for i ∈ range(innerEnd+1, length=cSize, step=cSize)
+        u[i] = op.parity*apply(flip(op.closureStencils[M-i+1]), v, i; stride=stride)/h^2
+    end
+end
+
+odd = -1
+even = 1
+
+struct D2{T} <: ConstantStencilOperator
     quadratureClosure::Vector{T}
     innerStencil::Stencil
     closureStencils::Vector{Stencil} # TBD: Should this be a tuple?
     eClosure::Vector{T}
     dClosure::Vector{T}
+    parity::Int
 end
 
 function closureSize(D::D2)::Int
@@ -38,6 +62,7 @@
         closureStencils,
         stringToVector(Float64, d["e"][1]),
         stringToVector(Float64, d["d1"][1]),
+        even
     )
 
     return d2