--- a/Notes.md	Fri Feb 10 08:36:56 2023 +0100
+++ b/Notes.md	Sat Feb 18 11:37:35 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	Sat Feb 18 11:37:35 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/equidistant_grid.jl	Fri Feb 10 08:36:56 2023 +0100
+++ b/src/Grids/equidistant_grid.jl	Sat Feb 18 11:37:35 2023 +0100
@@ -1,26 +1,77 @@
-"""
-    EquidistantGrid{Dim,T<:Real} <: Grid
+struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1,1}
+    points::R
+end
 
-`Dim`-dimensional equidistant grid with coordinates of type `T`.
+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!
+
+
 """
-struct EquidistantGrid{Dim,T<:Real} <: Grid
-    size::NTuple{Dim, Int}
-    limit_lower::NTuple{Dim, T}
-    limit_upper::NTuple{Dim, T}
+    spacing(grid::EquidistantGrid)
+
+The spacing between grid points.
+"""
+spacing(g::EquidistantGrid) = step(g.points)
+
 
-    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
+"""
+    inverse_spacing(grid::EquidistantGrid)
+
+The reciprocal of the spacing between grid points.
+"""
+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.
+
+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
 
 """
-    EquidistantGrid(size, limit_lower, limit_upper)
+    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
+
+
+
+
+
+
+
+"""
+    equidistant_grid(size::Dims, limit_lower, limit_upper)
 
 Construct an equidistant grid with corners at the coordinates `limit_lower` and
 `limit_upper`.
@@ -32,168 +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 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)...]
-
-# Review: 
-# Is it not strange that evalOn(::Grid) is non-lazy while evalOn(::EquidistantGrid) is?
-# Also: Change name to evalon or eval_on!!!!!!
-function evalOn(grid::EquidistantGrid, f::Function)
-    F(I...) = f(grid[I...]...)
+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
 
-    return LazyFunctionArray(F, size(grid))
-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)
-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)
+    return TensorGrid(gs...)
 end
 
 
 """
-    orthogonal_dims(grid::EquidistantGrid,dim)
+    equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
 
-Returns the dimensions of grid orthogonal to that of dim.
+Constructs a 1D equidistant grid.
 """
-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 equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
+	return equidistant_grid((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)
+    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	Fri Feb 10 08:36:56 2023 +0100
+++ b/src/Grids/grid.jl	Sat Feb 18 11:37:35 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	Sat Feb 18 11:37:35 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}()