--- a/diffOp.jl	Fri Jan 11 11:55:13 2019 +0100
+++ b/diffOp.jl	Fri Jan 11 11:55:48 2019 +0100
@@ -32,5 +32,47 @@
 function apply!(L::Laplace1D, u::AbstractVector, v::AbstractVector)
     h = grid.spacings(L.grid)[1]
     apply!(L.op, u, v, h)
+    u .= L.a * u
     return nothing
 end
+
+
+# Differential operator for a*d^2/dx^2 + a*d^2/dy^2
+struct Laplace2D <: DiffOp
+    grid
+    a
+    op
+end
+
+# u = L*v
+function apply!(L::Laplace2D, u::AbstractVector, v::AbstractVector)
+    u .= 0*u
+    h = grid.spacings(L.grid)
+
+    li = LinearIndices(L.grid.numberOfPointsPerDim)
+    n_x, n_y = L.grid.numberOfPointsPerDim
+
+
+    # For each x
+    temp = zeros(eltype(u), n_y)
+    for i ∈ 1:n_x
+
+        v_i = view(v, li[i,:])
+        apply!(L.op, temp, v_i, h[2])
+
+        u[li[i,:]] += temp
+    end
+
+    # For each y
+    temp = zeros(eltype(u), n_x)
+    for i ∈ 1:n_y
+        v_i = view(v, li[:,i])
+        apply!(L.op, temp, v_i, h[1])
+
+        u[li[:,i]] += temp
+    end
+
+    u .= L.a*u
+
+    return nothing
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/plotDerivative2d.jl	Fri Jan 11 11:55:48 2019 +0100
@@ -0,0 +1,16 @@
+g = sbp.grid.EquidistantGrid((100,75),((0, 0), (2pi, 3/2*pi)))
+op = sbp.readOperator("d2_4th.txt","h_4th.txt")
+Laplace = sbp.Laplace2D(g,1,op)
+
+init(x,y) = sin(x) + cos(y)
+v = sbp.grid.evalOn(g,init)
+
+u = zeros(length(v))
+
+sbp.apply!(Laplace,u,v)
+
+@show u
+@show u'*u
+
+sbp.grid.plotgridfunction(g,u)
+

--- a/sbpD2.jl	Fri Jan 11 11:55:13 2019 +0100
+++ b/sbpD2.jl	Fri Jan 11 11:55:48 2019 +0100
@@ -15,12 +15,16 @@
     end
 
     for i ∈ range(innerEnd+1, length=cSize)
-        u[i] = op.parity*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
+        u[i] = Int(op.parity)*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
     end
+
+    return nothing
 end
 
-odd = -1
-even = 1
+@enum Parity begin
+    odd = -1
+    even = 1
+end
 
 struct D2{T} <: ConstantStencilOperator
     quadratureClosure::Vector{T}
@@ -28,7 +32,7 @@
     closureStencils::Vector{Stencil} # TBD: Should this be a tuple?
     eClosure::Vector{T}
     dClosure::Vector{T}
-    parity::Int
+    parity::Parity
 end
 
 function closureSize(D::D2)::Int