--- a/grid.jl	Mon Dec 17 16:00:02 2018 +0100
+++ b/grid.jl	Tue Dec 18 13:27:25 2018 +0100
@@ -16,39 +16,62 @@
 
 abstract type BoundaryId end
 
-# Move to seperate file.
+# TODO: Move to seperate file.
+# Prefer to use UInt here, but printing UInt returns hex.
 struct EquidistantGrid <: Grid
-    nPointsPerDim::Vector{Int}
-    limits::Vector{Pair{Real, Real}}
-    function EquidistantGrid(nPointsPerDim, lims)
-        @assert length(lims) == length(nPointsPerDim)
+    numberOfPointsPerDim::Tuple
+    limits::NTuple{2,Tuple} # Stores the points at the lower and upper corner of the domain.
+                            # e.g (-1,0) and (1,2) for a domain of size (-1,1)x(0,2)
+
+    # General constructor
+    function EquidistantGrid(nPointsPerDim::Tuple, lims::NTuple{2,Tuple})
+        @assert length(nPointsPerDim) > 0
+        @assert count(x -> x > 0, nPointsPerDim) == length(nPointsPerDim)
+        @assert length(lims[1]) == length(nPointsPerDim)
+        @assert length(lims[2]) == length(nPointsPerDim)
+        # TODO: Assert that the same values are not passed in both lims[1] and lims[2]
+        #       i.e the domain length is positive for all dimensions
         return new(nPointsPerDim, lims)
     end
+    # 1D constructor which can be called as EquidistantGrid(m, (x_l,x_r))
+    function EquidistantGrid(nPointsPerDim::Int, lims::NTuple{2,Int})
+        return EquidistantGrid((nPointsPerDim,), ((lims[1],),(lims[2],)))
+    end
+
 end
 
 function numberOfDimensions(grid::EquidistantGrid)
-    return length(grid.nPointsPerDim)
+    return length(grid.numberOfPointsPerDim)
 end
 
 function numberOfPoints(grid::EquidistantGrid)
-    numberOfPoints = grid.nPointsPerDim[1];
-    for i = 2:length(grid.nPointsPerDim);
-        numberOfPoints = numberOfPoints*grid.nPointsPerDim[i]
+    numberOfPoints = grid.numberOfPointsPerDim[1];
+    for i = 2:length(grid.numberOfPointsPerDim);
+        numberOfPoints = numberOfPoints*grid.numberOfPointsPerDim[i]
     end
     return numberOfPoints
 end
 
+# TODO: Decide if spacings should be positive or if it is allowed to be negative
+#       If defined as positive, then need to do something extra when calculating the
+#       points. The current implementation works for
+function spacings(grid::EquidistantGrid)
+    h = Vector{Real}(undef, numberOfDimensions(grid))
+    for i ∈ eachindex(h)
+        h[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
+    end
+    return Tuple(h)
+end
+
 function points(grid::EquidistantGrid)
-    points = Vector{Real}(undef, numberOfPoints(grid))
-    for i = 1:numberOfDimensions(grid)
-        lims = limitsForDimension(grid,i)
-        points = range(lims.first, stop=lims.second, length=grid.nPointsPerDim[i])
+    nPoints = numberOfPoints(grid)
+    points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, nPoints)
+    cartesianIndices = CartesianIndices(grid.numberOfPointsPerDim)
+    for i ∈ 1:nPoints
+        ci = Tuple(cartesianIndices[i]) .-1
+        points[i] = grid.limits[1] .+ spacings(grid).*ci
     end
     return points
 end
 
-function limitsForDimension(grid::EquidistantGrid, dim::Int)
-    return grid.limits[dim]
 end
-
-end