--- a/EquidistantGrid.jl	Fri Jan 25 15:20:40 2019 +0100
+++ b/EquidistantGrid.jl	Thu Feb 21 16:27:28 2019 +0100
@@ -5,93 +5,52 @@
 # the domain is defined as (-1,1)x(0,2).
 
 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
-    numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then
-
+    size::NTuple{Dim, Int} # First coordinate direction stored first
     limit_lower::NTuple{Dim, T}
     limit_upper::NTuple{Dim, T}
+    inverse_spacing::NTuple{Dim, T} # The reciprocal of the grid spacing
 
     # General constructor
-    function EquidistantGrid(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
-        @assert all(nPointsPerDim.>0)
+    function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
+        @assert all(size.>0)
         @assert all(limit_upper.-limit_lower .!= 0)
-        return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper)
+        inverse_spacing = (size.-1)./abs.(limit_upper.-limit_lower)
+        return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing)
     end
+end
 
-    # # 1D constructor which can be called as EquidistantGrid(m, (xl,xr))
-    # function EquidistantGrid(nPointsPerDim::Integer, lims::NTuple{2,Real})
-    #     return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
-    # end
-
+function Base.eachindex(grid::EquidistantGrid)
+    CartesianIndices(grid.size)
 end
 
 # Returns the number of dimensions of an EquidistantGrid.
 #
 # @Input: grid - an EquidistantGrid
-# @Return: numberOfPoints - The number of dimensions
-function numberOfDimensions(grid::EquidistantGrid)
-    return length(grid.numberOfPointsPerDim)
-end
-
-# Computes the total number of points of an EquidistantGrid.
-#
-# @Input: grid - an EquidistantGrid
-# @Return: numberOfPoints - The total number of points
-function numberOfPoints(grid::EquidistantGrid)
-    return prod(grid.numberOfPointsPerDim)
+# @Return: dimension - The dimension of the grid
+function dimension(grid::EquidistantGrid)
+    return length(grid.size)
 end
 
-# Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
-# between two points for each coordinate direction.
+# Returns the spacing of the grid
 #
-# @Input: grid - an EquidistantGrid
-# @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
-function spacings(grid::EquidistantGrid)
-    return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
+function spacing(grid::EquidistantGrid)
+    return 1.0./grid.inverse_spacing
 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)
+# Computes the points of an EquidistantGrid as an array of tuples with
+# the same dimension as the grid.
 #
 # @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?
+    indices = Tuple.(CartesianIndices(grid.size))
+    h = spacing(grid)
+    return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
 end
 
 function pointsalongdim(grid::EquidistantGrid, dim::Integer)
-    @assert dim<=numberOfDimensions(grid)
+    @assert dim<=dimension(grid)
     @assert dim>0
-    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.numberOfPointsPerDim[dim])
+    points = range(grid.limit_lower[dim],stop=grid.limit_lower[dim],length=grid.size[dim])
 end
-
-using PyPlot, PyCall
-
-function plotgridfunction(grid::EquidistantGrid, gridfunction)
-    if numberOfDimensions(grid) == 1
-        plot(pointsalongdim(grid,1), gridfunction, linewidth=2.0)
-    elseif numberOfDimensions(grid) == 2
-        mx = grid.numberOfPointsPerDim[1]
-        my = grid.numberOfPointsPerDim[2]
-        X = repeat(pointsalongdim(grid,1),1,my)
-        Y = permutedims(repeat(pointsalongdim(grid,2),1,mx))
-        plot_surface(X,Y,reshape(gridfunction,mx,my));
-    else
-        error(string("Plot not implemented for dimension ", string(numberOfDimensions(grid))))
-    end
-end