abstract type StencilIndex end

function Base.getindex(si::StencilIndex, i::Int)
    return si.index[i]
end

struct LowerClosureIndex <: StencilIndex
    index::CartesianIndex
end

struct UpperClosureIndex  <: StencilIndex
    index::CartesianIndex
end

struct InteriorIndex  <: StencilIndex
    index::CartesianIndex
end

# TODO: This should take a Stencil or DiffOp so that we can extract all the
# indices in the closures.
# 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 ∈ ci
        I = Tuple(i)
        if any(islower, I)
            push!(lowerclosure, LowerClosureIndex(i))
            # TODO: Corner points should be in both Lower and Upper?
            #       Should they have a separate type?
        elseif any(isupper, I)
            push!(upperclosure, UpperClosureIndex(i))
        else
            push!(interior, InteriorIndex(i))
        end
    end
    return lowerclosure, upperclosure, interior
end