--- a/src/Grids/AbstractGrid.jl	Fri Jul 15 09:41:58 2022 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-"""
-     AbstractGrid
-
-Should implement
-    dim(grid::AbstractGrid)
-    points(grid::AbstractGrid)
-
-"""
-abstract type AbstractGrid end
-export AbstractGrid
-function dim end # TODO: Rename to Base.ndims instead? That's the name used for arrays.
-function points end
-export dim, points
-
-"""
-    evalOn(g::AbstractGrid, f::Function)
-
-Evaluate function f on the grid g
-"""
-function evalOn(g::AbstractGrid, f::Function)
-    F(x) = f(x...)
-    return F.(points(g))
-end
-export evalOn
-
-
-"""
-    dims(g::AbstractGrid)
-
-A range containing the dimensions of the grid
-"""
-dims(grid::AbstractGrid) = 1:dim(grid)
-export dims
\ No newline at end of file

--- a/src/Grids/EquidistantGrid.jl	Fri Jul 15 09:41:58 2022 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,205 +0,0 @@
-export EquidistantGrid
-export spacing
-export inverse_spacing
-export restrict
-export orthogonal_dims
-export boundary_identifiers
-export boundary_grid
-export boundary_size
-export refine
-export coarsen
-
-"""
-    EquidistantGrid{Dim,T<:Real} <: AbstractGrid
-
-`Dim`-dimensional equidistant grid with coordinates of type `T`.
-"""
-struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
-    size::NTuple{Dim, Int}
-    limit_lower::NTuple{Dim, T}
-    limit_upper::NTuple{Dim, T}
-
-    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
-
-"""
-    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)
-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
-
-"""
-    dim(grid::EquidistantGrid)
-
-The dimension of the grid.
-"""
-dim(::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)
-
-"""
-    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)
-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
-
-"""
-    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}()),dim(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}()
-
-"""
-    boundary_size(grid::EquidistantGrid, id::CartesianBoundary)
-
-Returns the size of the boundary of `grid` specified by `id`.
-"""
-function boundary_size(grid::EquidistantGrid, id::CartesianBoundary)
-	orth_dims = orthogonal_dims(grid, dim(id))
-    return  grid.size[orth_dims]
-end
-boundary_size(::EquidistantGrid{1,T},::CartesianBoundary{1}) where 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)
-"""
-function refine(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
-    new_sz = (sz .- 1).*r .+ 1
-    return EquidistantGrid{dim(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.
-
-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{dim(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end

--- a/src/Grids/Grids.jl	Fri Jul 15 09:41:58 2022 +0200
+++ b/src/Grids/Grids.jl	Fri Jul 15 09:54:15 2022 +0200
@@ -2,18 +2,33 @@
 
 using Sbplib.RegionIndices
 
-export BoundaryIdentifier, CartesianBoundary
+# Grid
+export Grid
+export dim
+export dims
+export points
+export evalOn
 
-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()
+# BoundaryIdentifier
+export BoundaryIdentifier
+export CartesianBoundary
+export dim
+export region
 
-export dim, region
+# EquidistantGrid
+export EquidistantGrid
+export spacing
+export inverse_spacing
+export restrict
+export orthogonal_dims
+export boundary_identifiers
+export boundary_grid
+export boundary_size
+export refine
+export coarsen
 
-include("AbstractGrid.jl")
-include("EquidistantGrid.jl")
-
-# TODO: Rename AbstractGrid to Grid and move definition here.
+include("grid.jl")
+include("boundary_identifier.jl")
+include("equidistant_grid.jl")
 
 end # module

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/equidistant_grid.jl	Fri Jul 15 09:54:15 2022 +0200
@@ -0,0 +1,206 @@
+
+"""
+    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}
+
+    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
+
+
+"""
+    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)
+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
+
+
+"""
+    dim(grid::EquidistantGrid)
+
+The dimension of the grid.
+"""
+dim(::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)
+
+
+"""
+    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)
+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
+
+
+"""
+    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}()),dim(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}()
+
+
+"""
+    boundary_size(grid::EquidistantGrid, id::CartesianBoundary)
+
+Returns the size of the boundary of `grid` specified by `id`.
+"""
+function boundary_size(grid::EquidistantGrid, id::CartesianBoundary)
+	orth_dims = orthogonal_dims(grid, dim(id))
+    return  grid.size[orth_dims]
+end
+boundary_size(::EquidistantGrid{1,T},::CartesianBoundary{1}) where 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)
+"""
+function refine(grid::EquidistantGrid, r::Int)
+    sz = size(grid)
+    new_sz = (sz .- 1).*r .+ 1
+    return EquidistantGrid{dim(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.
+
+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{dim(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/grid.jl	Fri Jul 15 09:54:15 2022 +0200
@@ -0,0 +1,28 @@
+"""
+     Grid
+
+Should implement
+    dim(grid::Grid)
+    points(grid::Grid)
+
+"""
+abstract type Grid end
+function dim end # TODO: Rename to Base.ndims instead? That's the name used for arrays.
+function points end
+
+"""
+    dims(g::Grid)
+
+A range containing the dimensions of the grid
+"""
+dims(grid::Grid) = 1:dim(grid)
+
+"""
+    evalOn(g::Grid, f::Function)
+
+Evaluate function f on the grid g
+"""
+function evalOn(g::Grid, f::Function)
+    F(x) = f(x...)
+    return F.(points(g))
+end
\ No newline at end of file

--- a/test/Grids/EquidistantGrid_test.jl	Fri Jul 15 09:41:58 2022 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,123 +0,0 @@
-using Sbplib.Grids
-using Test
-using Sbplib.RegionIndices
-
-
-@testset "EquidistantGrid" begin
-    @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid
-    @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid
-    # constuctor
-    @test_throws DomainError EquidistantGrid(0,0.0,1.0)
-    @test_throws DomainError EquidistantGrid(1,1.0,1.0)
-    @test_throws DomainError EquidistantGrid(1,1.0,-1.0)
-    @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,))
-
-    @testset "Base" begin
-        @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64
-        @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int
-        @test size(EquidistantGrid(4,0.0,1.0)) == (4,)
-        @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
-    end
-
-    # dim
-    @test dim(EquidistantGrid(4,0.0,1.0)) == 1
-    @test dim(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == 2
-
-    # spacing
-    @test [spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(1. /3,)...] atol=5e-13
-    @test [spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(0.5, 1.)...] atol=5e-13
-
-    # inverse_spacing
-    @test [inverse_spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(3.,)...] atol=5e-13
-    @test [inverse_spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(2, 1.)...] atol=5e-13
-
-    # points
-    g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11))
-    gp = points(g);
-    p = [(-1.,0.)      (-1.,7.11/2)   (-1.,7.11);
-         (-0.75,0.)    (-0.75,7.11/2) (-0.75,7.11);
-         (-0.5,0.)     (-0.5,7.11/2)  (-0.5,7.11);
-         (-0.25,0.)    (-0.25,7.11/2) (-0.25,7.11);
-         (0.,0.)       (0.,7.11/2)    (0.,7.11)]
-    for i ∈ eachindex(gp)
-        @test [gp[i]...] ≈ [p[i]...] atol=5e-13
-    end
-
-    # restrict
-    g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))
-    @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0)
-    @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0)
-
-    g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
-    @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0)
-    @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0)
-    @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0)
-    @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0))
-    @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
-    @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0))
-    @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0))
-
-    @testset "boundary_identifiers" begin
-        g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
-        bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(),
-                CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(),
-                CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}())
-        @test boundary_identifiers(g) == bids
-        @inferred boundary_identifiers(g)
-    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
-    end
-
-    @testset "refine" begin
-        @test refine(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}()
-        @test refine(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}()
-
-        g = EquidistantGrid((10,5),(0.,1.),(2.,3.))
-        @test refine(g, 1) == g
-        @test refine(g, 2) == EquidistantGrid((19,9),(0.,1.),(2.,3.))
-        @test refine(g, 3) == EquidistantGrid((28,13),(0.,1.),(2.,3.))
-    end
-
-    @testset "coarsen" begin
-        @test coarsen(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}()
-        @test coarsen(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}()
-
-        g = EquidistantGrid((7,13),(0.,1.),(2.,3.))
-        @test coarsen(g, 1) == g
-        @test coarsen(g, 2) == EquidistantGrid((4,7),(0.,1.),(2.,3.))
-        @test coarsen(g, 3) == EquidistantGrid((3,5),(0.,1.),(2.,3.))
-
-        @test_throws DomainError(4, "Size minus 1 must be divisible by the ratio.") coarsen(g, 4) == EquidistantGrid((3,5),(0.,1.),(2.,3.))
-    end
-end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Grids/equidistant_grid_test.jl	Fri Jul 15 09:54:15 2022 +0200
@@ -0,0 +1,123 @@
+using Sbplib.Grids
+using Test
+using Sbplib.RegionIndices
+
+
+@testset "EquidistantGrid" begin
+    @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid
+    @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid
+    # constuctor
+    @test_throws DomainError EquidistantGrid(0,0.0,1.0)
+    @test_throws DomainError EquidistantGrid(1,1.0,1.0)
+    @test_throws DomainError EquidistantGrid(1,1.0,-1.0)
+    @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,))
+
+    @testset "Base" begin
+        @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64
+        @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int
+        @test size(EquidistantGrid(4,0.0,1.0)) == (4,)
+        @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+    end
+
+    # dim
+    @test dim(EquidistantGrid(4,0.0,1.0)) == 1
+    @test dim(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == 2
+
+    # spacing
+    @test [spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(1. /3,)...] atol=5e-13
+    @test [spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(0.5, 1.)...] atol=5e-13
+
+    # inverse_spacing
+    @test [inverse_spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(3.,)...] atol=5e-13
+    @test [inverse_spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(2, 1.)...] atol=5e-13
+
+    # points
+    g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11))
+    gp = points(g);
+    p = [(-1.,0.)      (-1.,7.11/2)   (-1.,7.11);
+         (-0.75,0.)    (-0.75,7.11/2) (-0.75,7.11);
+         (-0.5,0.)     (-0.5,7.11/2)  (-0.5,7.11);
+         (-0.25,0.)    (-0.25,7.11/2) (-0.25,7.11);
+         (0.,0.)       (0.,7.11/2)    (0.,7.11)]
+    for i ∈ eachindex(gp)
+        @test [gp[i]...] ≈ [p[i]...] atol=5e-13
+    end
+
+    # restrict
+    g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))
+    @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0)
+    @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0)
+
+    g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
+    @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0)
+    @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0)
+    @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0)
+    @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0))
+    @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0))
+    @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0))
+    @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0))
+
+    @testset "boundary_identifiers" begin
+        g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0))
+        bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(),
+                CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(),
+                CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}())
+        @test boundary_identifiers(g) == bids
+        @inferred boundary_identifiers(g)
+    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
+    end
+
+    @testset "refine" begin
+        @test refine(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}()
+        @test refine(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}()
+
+        g = EquidistantGrid((10,5),(0.,1.),(2.,3.))
+        @test refine(g, 1) == g
+        @test refine(g, 2) == EquidistantGrid((19,9),(0.,1.),(2.,3.))
+        @test refine(g, 3) == EquidistantGrid((28,13),(0.,1.),(2.,3.))
+    end
+
+    @testset "coarsen" begin
+        @test coarsen(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}()
+        @test coarsen(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}()
+
+        g = EquidistantGrid((7,13),(0.,1.),(2.,3.))
+        @test coarsen(g, 1) == g
+        @test coarsen(g, 2) == EquidistantGrid((4,7),(0.,1.),(2.,3.))
+        @test coarsen(g, 3) == EquidistantGrid((3,5),(0.,1.),(2.,3.))
+
+        @test_throws DomainError(4, "Size minus 1 must be divisible by the ratio.") coarsen(g, 4) == EquidistantGrid((3,5),(0.,1.),(2.,3.))
+    end
+end