--- a/EquidistantGrid.jl	Fri Jan 25 15:26:47 2019 +0100
+++ b/EquidistantGrid.jl	Mon Feb 04 17:51:36 2019 +0100
@@ -49,6 +49,10 @@
     return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
 end
 
+function Base.eachindex(grid::EquidistantGrid)
+    CartesianIndices(grid.numberOfPointsPerDim)
+end
+
 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
 # such that points in the first coordinate direction varies first, then the second
 # and lastely the third (if applicable)
@@ -56,22 +60,10 @@
 # @Input: grid - an EquidistantGrid
 # @Return: points - the points of the grid.
 function points(grid::EquidistantGrid)
-    dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
-
-    points = Vector{typeof(dx̄)}(undef, numberOfPoints(grid))
-    # Compute the points based on their Cartesian indices and the signed
-    # grid spacings
-    cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
-    for i ∈ 1:numberOfPoints(grid)
-        ci = Tuple(cartesianIndices[i]) .-1
-        points[i] = grid.limit_lower .+ dx̄.*ci
-    end
-
-    # TBD: Keep? this? How do we want to represent points in 1D?
-    if numberOfDimensions(grid) == 1
-        points = broadcast(x -> x[1], points)
-    end
-    return points
+    # TODO: Make this return an abstract array?
+    physical_domain_size = (grid.limit_upper .- grid.limit_lower)
+    indices = Tuple.(CartesianIndices(grid.numberOfPointsPerDim))
+    return broadcast(I -> grid.limit_lower .+ physical_domain_size.*(I.-1), indices)
 end
 
 function pointsalongdim(grid::EquidistantGrid, dim::Integer)

--- a/diffOp.jl	Fri Jan 25 15:26:47 2019 +0100
+++ b/diffOp.jl	Mon Feb 04 17:51:36 2019 +0100
@@ -1,5 +1,6 @@
 abstract type DiffOp end
 
+# TBD: The "error("not implemented")" thing seems to be hiding good error information. How to fix that? Different way of saying that these should be implemented?
 function apply(D::DiffOp, v::AbstractVector, i::Int)
     error("not implemented")
 end
@@ -32,10 +33,12 @@
     error("not implemented")
 end
 
-# DiffOp must have a grid!!!
-function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector)
-    for i ∈ 1:Grid.numberOfPoints(D.grid)
-        u[i] = apply(D, v, i)
+abstract type DiffOpCartesian{Dim} <: DiffOp end
+
+# DiffOp must have a grid of dimension Dim!!!
+function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
+    for I ∈ eachindex(D.grid)
+        u[I] = apply(D, v, I)
     end
 
     return nothing
@@ -47,12 +50,16 @@
     return u
 end
 
-struct Laplace{Dim,T<:Real,N,M,K} <: DiffOp
+struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
     grid::Grid.EquidistantGrid{Dim,T}
     a::T
     op::D2{Float64,N,M,K}
 end
 
+function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
+    error("not implemented")
+end
+
 # u = L*v
 function apply(L::Laplace{1}, v::AbstractVector, i::Int)
     h = Grid.spacings(L.grid)[1]
@@ -63,19 +70,15 @@
 using UnsafeArrays
 
 # u = L*v
-function apply(L::Laplace{2}, v::AbstractVector, i::Int)
+function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, I::CartesianIndex{2})
     h = Grid.spacings(L.grid)
 
-    li = LinearIndices(L.grid.numberOfPointsPerDim)
-    ci = CartesianIndices(L.grid.numberOfPointsPerDim)
-    I = ci[i]
-
     # 2nd x-derivative
-    @inbounds vx = uview(v, uview(li,:,I[2]))
-    uᵢ  = apply(L.op, h[1], vx , I[1])
+    @inbounds vx = uview(v, :, I[2])
+    @inbounds uᵢ  = apply(L.op, h[1], vx , I[1])
     # 2nd y-derivative
-    @inbounds vy = uview(v, uview(li,I[1],:))
-    uᵢ += apply(L.op, h[2], vy, I[2])
+    @inbounds vy = uview(v, I[1], :)
+    @inbounds uᵢ += apply(L.op, h[2], vy, I[2])
 
     return uᵢ
 end