--- a/src/Grids/equidistant_grid.jl	Sat May 20 14:26:36 2023 +0200
+++ b/src/Grids/equidistant_grid.jl	Sun May 21 21:55:14 2023 +0200
@@ -1,83 +1,40 @@
-
 """
-    EquidistantGrid{Dim,T<:Real} <: Grid
-
-`Dim`-dimensional equidistant grid with coordinates of type `T`.
-"""
-struct EquidistantGrid{Dim,T<:Real} <: Grid
-    size::NTuple{Dim, Int}
-    limit_lower::NTuple{Dim, T}
-    limit_upper::NTuple{Dim, T}
+    EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1}
 
-    function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T}
-        if any(size .<= 0)
-            throw(DomainError("all components of size must be postive"))
-        end
-        if any(limit_upper.-limit_lower .<= 0)
-            throw(DomainError("all side lengths must be postive"))
-        end
-        return new{Dim,T}(size, limit_lower, limit_upper)
-    end
-end
+A one-dimensional equidistant grid. Most users are expected to use
+[`equidistant_grid`](@ref) for constructing equidistant grids.
+
+See also: [`equidistant_grid`](@ref)
 
 
+## Note
+The type of range used for the points can likely impact performance.
 """
-    EquidistantGrid(size, limit_lower, limit_upper)
-
-Construct an equidistant grid with corners at the coordinates `limit_lower` and
-`limit_upper`.
-
-The length of the domain sides are given by the components of
-`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and `limit_upper=(1,2)` the domain is defined
-as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative.
-
-The number of equidistantly spaced points in each coordinate direction are given
-by the tuple `size`.
-"""
-function EquidistantGrid(size, limit_lower, limit_upper)
-    return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper)
+struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1}
+    points::R
 end
 
-
-"""
-    EquidistantGrid{T}()
-
-Constructs a 0-dimensional grid.
-"""
-EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid
-
-
-"""
-    EquidistantGrid(size::Int, limit_lower::T, limit_upper::T)
+# Indexing interface
+Base.getindex(g::EquidistantGrid, i) = g.points[i]
+Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
+Base.firstindex(g::EquidistantGrid) = firstindex(g.points)
+Base.lastindex(g::EquidistantGrid) = lastindex(g.points)
 
-Convenience constructor for 1D grids.
-"""
-function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T
-	return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
-end
-
-Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
-
-Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
+# Iteration interface
+Base.iterate(g::EquidistantGrid) = iterate(g.points)
+Base.iterate(g::EquidistantGrid, state) = iterate(g.points, state)
 
-Base.size(g::EquidistantGrid) = g.size
+Base.IteratorSize(::Type{<:EquidistantGrid}) = Base.HasShape{1}()
+Base.length(g::EquidistantGrid) = length(g.points)
+Base.size(g::EquidistantGrid) = size(g.points)
 
-function Base.getindex(g::EquidistantGrid, I::Vararg{Int})
-    h = spacing(g)
-    return g.limit_lower .+ (I.-1).*h
-end
-
-Base.getindex(g::EquidistantGrid, I::CartesianIndex) = g[Tuple(I)...]
-# TBD: Can this method be removed if `EquidistantGrid` is an AbstractArray?
-
-Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
 
 """
     spacing(grid::EquidistantGrid)
 
 The spacing between grid points.
 """
-spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1)
+spacing(g::EquidistantGrid) = step(g.points)
 
 
 """
@@ -85,111 +42,98 @@
 
 The reciprocal of the spacing between grid points.
 """
-inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid)
+inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
 
 
-"""
-    points(grid::EquidistantGrid)
-
-The point of the grid as an array of tuples with the same dimension as the grid.
-The points are stored as [(x1,y1), (x1,y2), … (x1,yn);
-						  (x2,y1), (x2,y2), … (x2,yn);
-						  	⋮		 ⋮            ⋮
-						  (xm,y1), (xm,y2), … (xm,yn)]
-"""
-function points(grid::EquidistantGrid)
-    indices = Tuple.(CartesianIndices(grid.size))
-    h = spacing(grid)
-    return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
-end
+boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
+boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
+boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
 
 
 """
-    restrict(::EquidistantGrid, dim)
+    refine(g::EquidistantGrid, r::Int)
 
-Pick out given dimensions from the grid and return a grid for them.
+The grid where `g` is refined by the factor `r`. The factor is applied to the number of
+intervals, i.e., 1 less than the size of  `g`.
+
+See also: [`coarsen`](@ref)
 """
-function restrict(grid::EquidistantGrid, dim)
-    size = grid.size[dim]
-    limit_lower = grid.limit_lower[dim]
-    limit_upper = grid.limit_upper[dim]
-
-    return EquidistantGrid(size, limit_lower, limit_upper)
+function refine(g::EquidistantGrid, r::Int)
+    new_sz = (length(g) - 1)*r + 1
+    return EquidistantGrid(change_length(g.points, new_sz))
 end
 
+"""
+    coarsen(g::EquidistantGrid, r::Int)
 
-"""
-    orthogonal_dims(grid::EquidistantGrid,dim)
+The grid where `g` is coarsened by the factor `r`. The factor is applied to the number of
+intervals, i.e., 1 less than the size of `g`. If the number of
+intervals are not divisible by `r` an error is raised.
 
-Returns the dimensions of grid orthogonal to that of dim.
+See also: [`refine`](@ref)
 """
-function orthogonal_dims(grid::EquidistantGrid, dim)
-    orth_dims = filter(i -> i != dim, dims(grid))
-	if orth_dims == dims(grid)
-		throw(DomainError(string("dimension ",string(dim)," not matching grid")))
-	end
-    return orth_dims
+function coarsen(g::EquidistantGrid, r::Int)
+    if (length(g)-1)%r != 0
+        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
+    end
+
+    new_sz = (length(g) - 1)÷r + 1
+
+    return EquidistantGrid(change_length(g.points, new_sz))
 end
 
 
 """
-    boundary_identifiers(::EquidistantGrid)
+    equidistant_grid(size::Dims, limit_lower, limit_upper)
+
+Construct an equidistant grid with corners at the coordinates `limit_lower` and
+`limit_upper`.
+
+The length of the domain sides are given by the components of
+`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and
+`limit_upper=(1,2)` the domain is defined as `(-1,1)x(0,2)`. The side lengths
+of the grid are not allowed to be negative.
 
-Returns a tuple containing the boundary identifiers for the grid, stored as
-	(CartesianBoundary(1,Lower),
-	 CartesianBoundary(1,Upper),
-	 CartesianBoundary(2,Lower),
-	 ...)
+The number of equispaced points in each coordinate direction are given
+by the tuple `size`.
+
+Note: If `limit_lower` and `limit_upper` are integers and `size` would allow a
+completely integer grid, `equidistant_grid` will still return a floating point
+grid. This simlifies the implementation and avoids certain surprise
+behaviours.
 """
-boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
-
+function equidistant_grid(size::Dims, limit_lower, limit_upper)
+    gs = map(equidistant_grid, size, limit_lower, limit_upper)
+    return TensorGrid(gs...)
+end
 
 """
-    boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
+    equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
 
-Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
-orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
-grid is a zero-dimensional grid.
+Constructs a 1D equidistant grid.
 """
-function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
-    orth_dims = orthogonal_dims(grid, dim(id))
-    return restrict(grid, orth_dims)
+function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
+    if any(size .<= 0)
+        throw(DomainError("size must be postive"))
+    end
+
+    if any(limit_upper.-limit_lower .<= 0)
+        throw(DomainError("side length must be postive"))
+    end
+	return EquidistantGrid(range(limit_lower, limit_upper, length=size)) # TBD: Should it use LinRange instead?
 end
-boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
+
+CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
 
 
 """
-    refine(grid::EquidistantGrid, r::Int)
-
-Refines `grid` by a factor `r`. The factor is applied to the number of
-intervals which is 1 less than the size of the grid.
-
-See also: [`coarsen`](@ref)
-"""
-function refine(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
-    new_sz = (sz .- 1).*r .+ 1
-    return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end
-
-
-"""
-    coarsen(grid::EquidistantGrid, r::Int)
+    change_length(r::AbstractRange, n)
 
-Coarsens `grid` by a factor `r`. The factor is applied to the number of
-intervals which is 1 less than the size of the grid. If the number of
-intervals are not divisible by `r` an error is raised.
-
-See also: [`refine`](@ref)
+Change the length of `r` to `n`, keeping the same start and stop.
 """
-function coarsen(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
+function change_length end
 
-    if !all(n -> (n % r == 0), sz.-1)
-        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
-    end
-
-    new_sz = (sz .- 1).÷r .+ 1
-
-    return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end
+change_length(r::UnitRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
+change_length(r::StepRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
+change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
+change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)