--- a/Notes.md	Sun Feb 19 12:14:43 2023 +0100
+++ b/Notes.md	Sun Feb 19 22:07:57 2023 +0100
@@ -148,6 +148,7 @@
 
 ## Identifiers for regions
 The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probabily be included in the grid module. This allows new grid types to come with their own regions.
+We implement this by refactoring RegionIndices to be agnostic to the region types and then moving the actual types to Grids.
 
 ## Regions and tensormappings
 - [ ] Use a trait to indicate if a LazyTensor uses indices with regions.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/grid_refactor.md	Sun Feb 19 22:07:57 2023 +0100
@@ -0,0 +1,104 @@
+# Grids refactor
+
+## Mål
+  * Embedded grids
+  * Vektorvärda grid funktioner, speciellt koordinat-funktionen
+  * Olika griddar i olika riktningar?
+      * Tex periodiskt i en riktning och intervall i en.
+      * Tex ostrukturerat i ett par och strukturerat i en.
+  * Enkelt att evaluera på 0d-grid
+
+## Frågor
+
+
+### Ska `Grid` vara en AbstractArray?
+Efter som alla nät ska agera som en gridfunktion för koordinaterna måste man
+svara på frågan hur vi hanterar generellla gridfunktioner samtidigt.
+
+Några saker att förhålla sig till:
+  - Multiblock nät?
+  - Unstructured?
+  - Triangular structured grids?
+  - Non-simply connected?
+  - CG/DG-nät
+
+
+Om alla nät är någon slags AbstractArray så kan tillexempel ett multiblock nät vara en AbstractArray{<:Grid, 1} och motsvarande gridfunktioner AbstractArray{<:AbstractArray,1}.
+Där man alltså först indexerar för vilket när man vill åt och sedan indexerar nätet. Tex `mg[2][32,12]`.
+
+Ett ostrukturerat nät med till exempel trianglar skulle vi kunna se på samma sätt som ett multiblocknät. Antagligen har de två typerna av nät olika underliggande datastruktur anpassade efter ändamålet.
+
+Hur fungerar tankarna ovan om man har tex tensorprodukten av ett ostrukturerat nät och en ekvidistant nät?
+```julia
+m = Mesh2DTriangle(10)
+e = EqudistantGrid(range(1:10)
+
+e[4] # fourth point
+
+m[3][5] # Fifth node in third triangle
+m[3,5] # Fifth node in third triangle # Funkar bara om alla nät är samma, (stämmer inte i mb-fallet)
+
+g = TensorGrid(m, e)
+
+g[3,4][5] # ??
+g[3,4] # ??
+
+g[3,5,4] # ??
+
+
+
+```
+
+Alla grids kanske inte är AbstractArray? Måste de vara rektangulära? Det blir svårt för strukturerade trianglar och NSC-griddar. Men de ska i allafall vara indexerbara?
+
+Man skulle kunna utesluta MultiblockGrid i tensorgrids
+
+CG-nät och DG-nät blir olika.
+På CG-nät kanske man både vill indexera noder och trianglar beroende på vad man håller på med?
+
+#### Försök till slutsater
+ * Multiblock-nät indexeras i två nivåer tex `g[3][3,4]`
+     * Vi struntar i att implementera multiblock-nät som en del av ett tensorgrid till att börja med.
+ * En grid kan inte alltid vara en AbstractArray eftersom till exempel ett NCS eller strukturerad triangel inte har rätt form.
+ * Om vi har nod-indexerade ostrukturerade nät borde de fungera med TensorGrid.
+ *
+
+### Kan vi introducera 1d griddar och tensorgriddar?
+  * Vanligt intervallgrid
+  * Periodiskt grid
+  * 0-dimensionellt grid
+
+Det skulle kunna lösa frågan med embedded-griddar
+och frågan om mixade grids.
+
+En svårighet kan vara att implemetera indexeringen för tensorgridden. Är det
+lätt att göra den transparent för kompilatorn?
+
+Läs i notes.md om: Grids embedded in higher dimensions
+
+Periodiska gridfunktioner? Eller ska vi bara skita i det och låta användaren
+hantera det själv? Blir det krångligt i högre dimensioner?
+
+
+### Hur ska vi hantera gridfunktioner med flera komponenter?
+Tex koordinat-funktionen på nätet?
+
+Funkar det att ha StaticArray som element type?
+    Kan man köra rakt på med en operator då?
+
+Vi skulle också kunna använda ArraysOfArrays. Skillnaden blir att den elementtypen är Array{T,N}. Vilket betyder att den är by reference?
+    Frågan är om kompilatorn kommer att bry sig... Jag skulle luta mot StaticArrays for safety
+
+Läs i notes.md om: Vector valued grid functions
+
+## Notes from pluto notebook
+- Är det dåligt att använda ndims om antalet index inte matchar?
+   - Tex ostrukturerat grid
+   - Baserat på hur Array fungerar och ndims docs borde den nog returnera
+     antalet index och inte griddens dimension
+   - Å andra sidan kanske man inte behöver veta viken dimension det är? Fast det känns konstigt...
+- Should we keep the `points` function?
+   - Maybe it's better to support `collect` on the grid?
+- How to handle boundary identifiers and boundary grids on `TensorGrid`?
+
+

--- a/src/Grids/Grids.jl	Sun Feb 19 12:14:43 2023 +0100
+++ b/src/Grids/Grids.jl	Sun Feb 19 22:07:57 2023 +0100
@@ -1,6 +1,7 @@
 module Grids
 
 using Sbplib.RegionIndices
+using Sbplib.LazyTensors
 
 # Grid
 export Grid

--- a/src/Grids/equidistant_grid.jl	Sun Feb 19 12:14:43 2023 +0100
+++ b/src/Grids/equidistant_grid.jl	Sun Feb 19 22:07:57 2023 +0100
@@ -1,28 +1,77 @@
+struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1,1}
+    points::R
+end
+
+Base.eltype(g::EquidistantGrid{T}) where T = T
+Base.getindex(g::EquidistantGrid, i) = g.points[i]
+Base.size(g::EquidistantGrid) = size(g.points)
+Base.length(g::EquidistantGrid) = length(g.points)
+Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
+
+# TODO: Make sure collect works!
+
+
+"""
+    spacing(grid::EquidistantGrid)
+
+The spacing between grid points.
+"""
+spacing(g::EquidistantGrid) = step(g.points)
+
 
 """
-    EquidistantGrid{Dim,T<:Real} <: Grid
+    inverse_spacing(grid::EquidistantGrid)
 
-`Dim`-dimensional equidistant grid with coordinates of type `T`.
+The reciprocal of the spacing between grid points.
 """
-struct EquidistantGrid{Dim,T<:Real} <: Grid
-    size::NTuple{Dim, Int}
-    limit_lower::NTuple{Dim, T}
-    limit_upper::NTuple{Dim, T}
+inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
+
+
+boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
+boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
+boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
+
+
+"""
+    refine(g::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.
 
-    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)
+See also: [`coarsen`](@ref)
+"""
+function refine(g::EquidistantGrid, r::Int)
+    new_sz = (length(g) - 1)*r + 1
+    return EquidistantGrid(change_length(g.points, new_sz))
+end
+
+"""
+    coarsen(grid::EquidistantGrid, r::Int)
+
+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)
+"""
+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
 
 
+
+
+
+
+
 """
-    EquidistantGrid(size, limit_lower, limit_upper)
+    equidistant_grid(size::Dims, limit_lower, limit_upper)
 
 Construct an equidistant grid with corners at the coordinates `limit_lower` and
 `limit_upper`.
@@ -34,158 +83,33 @@
 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)
-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)
-
-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)
-
-Base.size(g::EquidistantGrid) = g.size
-
-Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
-
-
+function equidistant_grid(size::Dims, limit_lower, limit_upper)
+    gs = map(size, limit_lower, limit_upper) do s,l,u
+        EquidistantGrid(range(l, u, length=s)) # TBD: Should it use LinRange instead?
+    end
 
-
-
-"""
-    spacing(grid::EquidistantGrid)
-
-The spacing between grid points.
-"""
-spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1)
-
-
-"""
-    inverse_spacing(grid::EquidistantGrid)
-
-The reciprocal of the spacing between grid points.
-"""
-inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid)
-
-
-"""
-    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)
+    return TensorGrid(gs...)
 end
 
 
 """
-    restrict(::EquidistantGrid, dim)
+    equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
 
-Pick out given dimensions from the grid and return a grid for them.
+Constructs a 1D equidistant grid.
 """
-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 equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
+	return equidistant_grid((size,),(limit_lower,),(limit_upper,))
 end
 
 
-"""
-    orthogonal_dims(grid::EquidistantGrid,dim)
-
-Returns the dimensions of grid orthogonal to that of dim.
-"""
-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
-end
-
-
-"""
-    boundary_identifiers(::EquidistantGrid)
-
-Returns a tuple containing the boundary identifiers for the grid, stored as
-	(CartesianBoundary(1,Lower),
-	 CartesianBoundary(1,Upper),
-	 CartesianBoundary(2,Lower),
-	 ...)
-"""
-boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
-
 
 """
-    boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
+    change_length(::AbstractRange, n)
 
-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.
-"""
-function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
-    orth_dims = orthogonal_dims(grid, dim(id))
-    return restrict(grid, orth_dims)
-end
-boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
-
-
-"""
-    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)
+Change the length of a range to `n`, keeping the same start and stop.
 """
-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)
-
-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.
+function change_length(::AbstractRange, n) end
 
-See also: [`refine`](@ref)
-"""
-function coarsen(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
-
-    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::LinRange, n) = LinRange(r[begin], r[end], n)
+change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
+# TODO: Test the above

--- a/src/Grids/grid.jl	Sun Feb 19 12:14:43 2023 +0100
+++ b/src/Grids/grid.jl	Sun Feb 19 22:07:57 2023 +0100
@@ -1,27 +1,76 @@
 """
-     Grid
+     Grid{T,D,RD} <: AbstractArray{T,D}
+
+The top level type for grids.
 
 Should implement
-    Base.ndims(grid::Grid)
-    points(grid::Grid)
+# TBD:
+"""
+#TBD: Does all the kinds of grids we want fit with this interface?
+# Multigrid?
+# Unstructured?
+# Triangular structured grids?
+# Non-simply connected?
+#
+# Maybe it shouldn't be an abstract array after all?
+abstract type Grid{T,D,RD} <: AbstractArray{T,D} end
+
+
+Base.ndims(::Grid{T,D,RD}) where {T,D,RD} = D # nidms borde nog vara antalet index som används för att indexera nätet. Snarare än vilken dimension nätet har (tänk ostrukturerat)
+nrangedims(::Grid{T,D,RD}) where {T,D,RD} = RD
+Base.eltype(::Grid{T,D,RD}) where {T,D,RD} = T # vad ska eltype vara? Inte T väl... en vektor? SVector{T,D}?
+
+function eval_on(::Grid) end # TODO: Should return a LazyArray and index the grid
+function refine(::Grid) end
+function coarsen(::Grid) end # Should this be here? What if it is not possible?
+
+abstract type BoundaryId end
 
 """
-abstract type Grid end
-function points end
+# TODO
+"""
+function boundary_identifiers(::Grid) end
+"""
+# TODO
+"""
+function boundary_grid(::Grid, ::BoundaryId) end
+
+
+# TODO: Make sure that all grids implement all of the above.
 
 """
     dims(grid::Grid)
 
-A range containing the dimensions of `grid`
+Enumerate the dimensions of the grid.
 """
 dims(grid::Grid) = 1:ndims(grid)
 
+
+
+# TBD: New file grid_functions.jl?
+
 """
-    evalOn(grid::Grid, f::Function)
+    getcomponent(gfun, I::Vararg{Int})
 
-Evaluate function `f` on `grid`
+Return one of the components of gfun as a grid function.
 """
-function evalOn(grid::Grid, f::Function)
-    F(x) = f(x...)
-    return F.(points(grid))
+# Should it be lazy? Could it be a view?
+function getcomponent(gfun, I::Vararg{Int}) end
+# function getcomponent(gfun, s::Symbol) end ?
+
+# TBD: New file zero_dim_grid.jl?
+struct ZeroDimGrid{T,S,RD} <: Grid{T,0,RD}
+    p::S
+
+    function ZeroDimGrid(p)
+        T = eltype(p)
+        S = typeof(p)
+        RD = length(p)
+        return new{T,S,RD}(p)
+    end
 end
+
+Base.size(g::ZeroDimGrid) = ()
+Base.getindex(g::ZeroDimGrid) = g.p
+Base.eachindex(g::ZeroDimGrid) = CartesianIndices(())
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/tensor_grid.jl	Sun Feb 19 22:07:57 2023 +0100
@@ -0,0 +1,87 @@
+struct TensorGrid{T,D,RD,GT<:NTuple{N,Grid} where N} <: Grid{T,D,RD}
+    grids::GT
+
+    function TensorGrid(gs...)
+        T = eltype(gs[1]) # All gs should have the same T
+        D = sum(ndims,gs)
+        RD = sum(nrangedims, gs)
+
+        return new{T,D,RD,typeof(gs)}(gs)
+    end
+end
+
+function Base.size(g::TensorGrid)
+    return concatenate_tuples(size.(g.grids)...)
+end
+
+function Base.getindex(g::TensorGrid, I...)
+    szs = ndims.(g.grids)
+
+    Is = split_tuple(I, szs)
+    ps = map((g,I)->SVector(g[I...]), g.grids, Is)
+
+    return vcat(ps...)
+end
+
+IndexStyle(::TensorGrid) = IndexCartesian()
+
+function Base.eachindex(g::TensorGrid)
+    szs = concatenate_tuples(size.(g.grids)...)
+    return CartesianIndices(szs)
+end
+
+
+
+## Pre refactor code:
+"""
+    orthogonal_dims(grid::EquidistantGrid,dim)
+
+Returns the dimensions of grid orthogonal to that of dim.
+"""
+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
+end
+
+"""
+    restrict(::EquidistantGrid, dim)
+
+Pick out given dimensions from the grid and return a grid for them.
+"""
+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)
+end
+
+
+
+"""
+    boundary_identifiers(::EquidistantGrid)
+
+Returns a tuple containing the boundary identifiers for the grid, stored as
+    (CartesianBoundary(1,Lower),
+     CartesianBoundary(1,Upper),
+     CartesianBoundary(2,Lower),
+     ...)
+"""
+boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
+
+
+"""
+    boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
+
+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.
+"""
+function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
+    orth_dims = orthogonal_dims(grid, dim(id))
+    return restrict(grid, orth_dims)
+end
+boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()

--- a/test/Grids/equidistant_grid_test.jl	Sun Feb 19 12:14:43 2023 +0100
+++ b/test/Grids/equidistant_grid_test.jl	Sun Feb 19 22:07:57 2023 +0100
@@ -1,6 +1,7 @@
 using Sbplib.Grids
 using Test
 using Sbplib.RegionIndices
+using Sbplib.LazyTensors
 
 
 @testset "EquidistantGrid" begin
@@ -44,6 +45,28 @@
         end
     end
 
+    @testset "getindex" begin
+        g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11))
+        @test g[1,1] == (-1.0,0.0)
+        @test g[1,3] == (-1.0,7.11)
+        @test g[5,1] == (0.0,0.0)
+        @test g[5,3] == (0.0,7.11)
+
+        @test g[4,2] == (-0.25,7.11/2)
+
+        @test g[CartesianIndex(1,3)] == (-1.0,7.11)
+    end
+
+    @testset "evalOn" begin
+        g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))
+
+        @test evalOn(g, (x,y) -> 0.) isa LazyArray
+        @test evalOn(g, (x,y) -> 0.) == fill(0., (5,3))
+
+        f(x,y) = sin(x)*cos(y)
+        @test evalOn(g, f) == map(p->f(p...), points(g))
+    end
+
     @testset "restrict" begin
         g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))
         @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0)
@@ -69,36 +92,36 @@
     end
 
     @testset "boundary_grid" begin
-            @testset "1D" begin
-                g = EquidistantGrid(5,0.0,2.0)
-                (id_l, id_r) = boundary_identifiers(g)
-                @test boundary_grid(g,id_l) == EquidistantGrid{Float64}()
-                @test boundary_grid(g,id_r) == EquidistantGrid{Float64}()
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}())
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}())
-            end
-            @testset "2D" begin
-                g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
-                (id_w, id_e, id_s, id_n) = boundary_identifiers(g)
-                @test boundary_grid(g,id_w) == restrict(g,2)
-                @test boundary_grid(g,id_e) == restrict(g,2)
-                @test boundary_grid(g,id_s) == restrict(g,1)
-                @test boundary_grid(g,id_n) == restrict(g,1)
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
-            end
-            @testset "3D" begin
-                g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
-                (id_w, id_e,
-                 id_s, id_n,
-                 id_t, id_b) = boundary_identifiers(g)
-                @test boundary_grid(g,id_w) == restrict(g,[2,3])
-                @test boundary_grid(g,id_e) == restrict(g,[2,3])
-                @test boundary_grid(g,id_s) == restrict(g,[1,3])
-                @test boundary_grid(g,id_n) == restrict(g,[1,3])
-                @test boundary_grid(g,id_t) == restrict(g,[1,2])
-                @test boundary_grid(g,id_b) == restrict(g,[1,2])
-                @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
-            end
+        @testset "1D" begin
+            g = EquidistantGrid(5,0.0,2.0)
+            (id_l, id_r) = boundary_identifiers(g)
+            @test boundary_grid(g,id_l) == EquidistantGrid{Float64}()
+            @test boundary_grid(g,id_r) == EquidistantGrid{Float64}()
+            @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}())
+            @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}())
+        end
+        @testset "2D" begin
+            g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
+            (id_w, id_e, id_s, id_n) = boundary_identifiers(g)
+            @test boundary_grid(g,id_w) == restrict(g,2)
+            @test boundary_grid(g,id_e) == restrict(g,2)
+            @test boundary_grid(g,id_s) == restrict(g,1)
+            @test boundary_grid(g,id_n) == restrict(g,1)
+            @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
+        end
+        @testset "3D" begin
+            g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
+            (id_w, id_e,
+             id_s, id_n,
+             id_t, id_b) = boundary_identifiers(g)
+            @test boundary_grid(g,id_w) == restrict(g,[2,3])
+            @test boundary_grid(g,id_e) == restrict(g,[2,3])
+            @test boundary_grid(g,id_s) == restrict(g,[1,3])
+            @test boundary_grid(g,id_n) == restrict(g,[1,3])
+            @test boundary_grid(g,id_t) == restrict(g,[1,2])
+            @test boundary_grid(g,id_b) == restrict(g,[1,2])
+            @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
+        end
     end
 
     @testset "refine" begin