--- a/EquidistantGrid.jl	Thu Jan 17 13:53:01 2019 +0100
+++ b/EquidistantGrid.jl	Thu Jan 17 15:19:33 2019 +0100
@@ -3,27 +3,24 @@
 # by the exterior product of the vectors obtained by projecting (x̄₂-x̄₁) onto
 # the coordinate directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2)
 # the domain is defined as (-1,1)x(0,2).
-struct EquidistantGrid <: AbstractGrid
-    numberOfPointsPerDim::Tuple # First coordinate direction stored first, then
-                                # second, then third.
-    limits::NTuple{2,Tuple} # Stores the two points which defines the range of
-                            # the e.g (-1,0) and (1,2) for a domain of size
-                            # (-1,1)x(0,2)
+
+struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
+    numberOfPointsPerDim::NTuple{Dim, Int} # First coordinate direction stored first, then
+
+    limit_lower::NTuple{Dim, T}
+    limit_upper::NTuple{Dim, T}
 
     # 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)
+    function EquidistantGrid(nPointsPerDim::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T
+        @assert all(nPointsPerDim.>0)
+        @assert all(limit_upper.-limit_lower .!= 0)
+        return new{Dim,T}(nPointsPerDim, limit_lower, limit_upper)
     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
+
+    # # 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
 
 end
 
@@ -40,11 +37,7 @@
 # @Input: grid - an EquidistantGrid
 # @Return: numberOfPoints - The total number of points
 function numberOfPoints(grid::EquidistantGrid)
-    numberOfPoints = grid.numberOfPointsPerDim[1];
-    for i = 2:length(grid.numberOfPointsPerDim);
-        numberOfPoints = numberOfPoints*grid.numberOfPointsPerDim[i]
-    end
-    return numberOfPoints
+    return prod(grid.numberOfPointsPerDim)
 end
 
 # Computes the grid spacing of an EquidistantGrid, i.e the unsigned distance
@@ -53,10 +46,7 @@
 # @Input: grid - an EquidistantGrid
 # @Return: h̄ - Grid spacing for each coordinate direction stored in a tuple.
 function spacings(grid::EquidistantGrid)
-    a = grid.limits[1]
-    b = grid.limits[2]
-
-    return abs.(b.-a)./(grid.numberOfPointsPerDim.-1)
+    return abs.(grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
 end
 
 # Computes the points of an EquidistantGrid as a vector of tuples. The vector is ordered
@@ -66,21 +56,17 @@
 # @Input: grid - an EquidistantGrid
 # @Return: points - the points of the grid.
 function points(grid::EquidistantGrid)
-    # Compute signed grid spacings
-    dx̄ = Vector{Real}(undef, numberOfDimensions(grid))
-    for i ∈ eachindex(dx̄)
-        dx̄[i] = (grid.limits[2][i]-grid.limits[1][i])/(grid.numberOfPointsPerDim[i]-1)
-    end
-    dx̄ = Tuple(dx̄)
+    dx̄ = (grid.limit_upper.-grid.limit_lower)./(grid.numberOfPointsPerDim.-1)
 
-    points = Vector{NTuple{numberOfDimensions(grid),Real}}(undef, numberOfPoints(grid))
+    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.limits[1] .+ dx̄.*ci
+        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)