--- a/Notes.md	Tue Feb 21 21:01:46 2023 +0100
+++ b/Notes.md	Tue Feb 21 21:07:39 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	Tue Feb 21 21:07:39 2023 +0100
@@ -0,0 +1,107 @@
+# 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?
+
+
+Om griddarna inte ska vara AbstractArray finns det många andra ställen som blir konstiga om de är AbstractArray. TensorApplication?! LazyArrays?! Är alla saker vi jobbar med egentligen mer generella object? Finns det något sätt att uttrycka koden så att man kan välja?
+
+#### 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	Tue Feb 21 21:01:46 2023 +0100
+++ b/src/Grids/Grids.jl	Tue Feb 21 21:07:39 2023 +0100
@@ -1,6 +1,7 @@
 module Grids
 
 using Sbplib.RegionIndices
+using Sbplib.LazyTensors
 
 # Grid
 export Grid
@@ -27,5 +28,8 @@
 include("grid.jl")
 include("boundary_identifier.jl")
 include("equidistant_grid.jl")
+include("zero_dim_grid.jl")
+
+abstract type BoundaryIdentifier end
 
 end # module

--- a/src/Grids/boundary_identifier.jl	Tue Feb 21 21:01:46 2023 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-
-abstract type BoundaryIdentifier end
-
-struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end
-dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim
-region(::CartesianBoundary{Dim, R}) where {Dim, R} = R()
\ No newline at end of file

--- a/src/Grids/equidistant_grid.jl	Tue Feb 21 21:01:46 2023 +0100
+++ b/src/Grids/equidistant_grid.jl	Tue Feb 21 21:07:39 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	Tue Feb 21 21:01:46 2023 +0100
+++ b/src/Grids/grid.jl	Tue Feb 21 21:07:39 2023 +0100
@@ -1,27 +1,57 @@
 """
-     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 Grid end
-function points end
+# TODO
+"""
+function boundary_identifiers end
+"""
+# TODO
+"""
+function boundary_grid 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)
 
-"""
-    evalOn(grid::Grid, f::Function)
+
 
-Evaluate function `f` on `grid`
+# TBD: New file grid_functions.jl?
+
 """
-function evalOn(grid::Grid, f::Function)
-    F(x) = f(x...)
-    return F.(points(grid))
-end
+    getcomponent(gfun, I::Vararg{Int})
+
+Return one of the components of gfun as a grid function.
+"""
+# Should it be lazy? Could it be a view?
+function getcomponent(gfun, I::Vararg{Int}) end
+# function getcomponent(gfun, s::Symbol) end ?

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/tensor_grid.jl	Tue Feb 21 21:07:39 2023 +0100
@@ -0,0 +1,62 @@
+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
+
+
+struct TensorBoundary{N, BID<:BoundaryIdentifier} <: BoundaryIdentifier end
+grid_id(::TensorBoundary{N, BID}) where {N, BID} = N
+boundary_id(::TensorBoundary{N, BID}) where {N, BID} = BID()
+
+
+"""
+    boundary_identifiers(::TensorGrid)
+
+Returns a tuple containing the boundary identifiers for the grid.
+"""
+function boundary_identifiers(g::TensorGrid)
+    n = length(g.grids)
+    per_grid = map(eachindex(g.grids)) do i
+        return map(bid -> TensorBoundary{i, bid}(), boundary_identifiers(g.grids[i]))
+    end
+    return concatenate_tuples(per_grid...)
+end
+
+
+"""
+    boundary_grid(grid::TensorGrid, id::TensorBoundary)
+
+The grid for the boundary specified by `id`.
+"""
+function boundary_grid(g::TensorGrid, bid::TensorBoundary)
+    local_boundary_grid = boundary_grid(g.grids[grid_id(bid)], boundary_id(bid))
+    new_grids = Base.setindex(g.grids, local_boundary_grid, grid_id(bid))
+    return TensorGrid(new_grids...)
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/zero_dim_grid.jl	Tue Feb 21 21:07:39 2023 +0100
@@ -0,0 +1,14 @@
+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(())

--- a/test/Grids/equidistant_grid_test.jl	Tue Feb 21 21:01:46 2023 +0100
+++ b/test/Grids/equidistant_grid_test.jl	Tue Feb 21 21:07:39 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