--- a/StencilIndex.jl	Mon Feb 04 09:13:48 2019 +0100
+++ b/StencilIndex.jl	Mon Feb 04 16:09:07 2019 +0100
@@ -1,45 +1,55 @@
-abstract type StencilIndex end
+abstract type Region end
+struct Interior <: Region end
+struct Lower    <: Region end
+struct Upper    <: Region end
 
-function Base.getindex(si::StencilIndex, i::Int)
-    return si.gridindex[i]
+struct StencilIndex{R<:Region, T<:Integer}
+    localindex::CartesianIndex
+    globalindex::T
+
+    StencilIndex{R}(li::CartesianIndex, gi::T) where {R<:Region,T<:Integer} = new{R, T}(li, gi)
+    StencilIndex(li::CartesianIndex, gi::T, ::Type{R}) where {R<:Region,T<:Integer} = StencilIndex{R}(li, gi)
+
+    # Index(t::Tuple{T, Type{R}}) where {R<:Region,T<:Integer} = Index{t[2]}(t[1])
+    # Above doesn't work, below does but is less type strict
+    #Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2]}(t[1])
 end
 
-struct LowerClosureIndex <: StencilIndex
-    globalindex::Integer
-    gridindex::CartesianIndex
-end
-
-struct UpperClosureIndex <: StencilIndex
-    globalindex::Integer
-    gridindex::CartesianIndex
-end
-
-struct InteriorIndex <: StencilIndex
-    globalindex::Integer
-    gridindex::CartesianIndex
+function Base.getindex(si::StencilIndex, i::Int)
+     return si.localindex[i]
 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.
+#Index(t::Vararg{Tuple{T, DataType}}) where T = Index.(t)
 # TODO: Where to place this function?
-function stencilindices(grid::Grid.EquidistantGrid)
-    lowerclosure = Vector{LowerClosureIndex}(undef, 0)
-    upperclosure = Vector{UpperClosureIndex}(undef, 0)
-    interior = Vector{InteriorIndex}(undef, 0)
-    # TODO: Fix such that the indices of the entire closure width is included.
-    islower = x -> (x == 1)
-    isupper = x -> (x in grid.numberOfPointsPerDim)
-    ci = CartesianIndices(grid.numberOfPointsPerDim)
-    for i ∈ 1:Grid.numberOfPoints(grid)
-        I = Tuple(ci[i])
-        if any(islower, I)
-            push!(lowerclosure, LowerClosureIndex(i,ci[i]))
-        elseif any(isupper, I)
-            push!(upperclosure, UpperClosureIndex(i,ci[i]))
-        else
-            push!(interior, InteriorIndex(i,ci[i]))
+
+function stencilindices(diffOp)
+    N = diffOp.grid.numberOfPointsPerDim
+
+    lowerclosure = Vector{Vector{StencilIndex{Lower, Int64}}}(undef,0)
+    upperclosure = Vector{Vector{StencilIndex{Upper, Int64}}}(undef,0)
+    interior = Vector{Vector{StencilIndex{Interior, Int64}}}(undef,0)
+    cSize = closureSize(diffOp.op)
+    ci = CartesianIndices(diffOp.grid.numberOfPointsPerDim)
+
+    # TODO: Loop over all points or one loop for each region?
+    for j = 1:Grid.numberOfDimensions(diffOp.grid)
+        templ = Vector{StencilIndex{Lower,Int64}}(undef, 0)
+        tempu = Vector{StencilIndex{Upper,Int64}}(undef, 0)
+        tempi = Vector{StencilIndex{Interior,Int64}}(undef, 0)
+        for i ∈ 1:Grid.numberOfPoints(diffOp.grid)
+            val = ci[i][j]
+            if val ∈ range(1; length=cSize)
+                push!(templ, StencilIndex{Lower}(ci[i],i))
+            elseif val ∈ range(N[j] - cSize+1, length=cSize)
+                push!(tempu, StencilIndex{Upper}(ci[i],i))
+            else
+                push!(tempi, StencilIndex{Interior}(ci[i],i))
+            end
         end
+        push!(lowerclosure,templ)
+        push!(upperclosure,tempu)
+        push!(interior,tempi)
     end
     return lowerclosure, upperclosure, interior
 end
+