--- a/Notes.md	Thu Nov 23 16:08:47 2023 +0100
+++ b/Notes.md	Fri Dec 01 11:52:26 2023 +0100
@@ -313,3 +313,15 @@
 Could the implementation of LazyOuterProduct be simplified by making it a
 struct containing two or more LazyTensors? (using split_tuple in a similar way
 as TensorGrid)
+
+## Implementation of boundary_indices for more complex grids
+To represent boundaries of for example tet-elements we can use a type `IndexCollection` to index a grid function directly.
+
+```julia
+I = IndexCollection(...)
+v[I]
+```
+
+* This would impact how tensor grid works.
+* To make things homogenous maybe these index collections should be used for the more simple grids too.
+* The function `to_indices` from Base could be useful to implement for `IndexCollection`

--- a/src/Grids/Grids.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/src/Grids/Grids.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -9,32 +9,27 @@
 export Grid
 export coordinate_size
 export component_type
+export grid_id
+export boundary_id
+export boundary_indices
+export boundary_identifiers
+export boundary_grid
+export eval_on
+export coarsen
+export getcomponent
+export refine
+
+export BoundaryIdentifier
+export TensorGridBoundary
+export CartesianBoundary
 
 export TensorGrid
 export ZeroDimGrid
 
-export TensorGridBoundary
-
-export grid_id
-export boundary_id
-
-export eval_on
-export getcomponent
-
-# BoundaryIdentifier
-export BoundaryIdentifier
-
-
-# EquidistantGrid
 export EquidistantGrid
+export inverse_spacing
 export spacing
-export inverse_spacing
-export boundary_identifiers
-export boundary_grid
-export refine
-export coarsen
 export equidistant_grid
-export CartesianBoundary
 
 
 # CurvilinearGrid

--- a/src/Grids/equidistant_grid.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/src/Grids/equidistant_grid.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -51,7 +51,8 @@
 boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
 boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
 boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
-
+boundary_indices(g::EquidistantGrid, id::Lower) = (1,)
+boundary_indices(g::EquidistantGrid, id::Upper) = (length(g),)
 
 """
     refine(g::EquidistantGrid, r::Int)

--- a/src/Grids/grid.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/src/Grids/grid.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -74,6 +74,18 @@
 # TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff?
 
 """
+    boundary_indices(g::Grid, id::BoundaryIdentifier)
+
+A collection of indices corresponding to the boundary with given id. For grids
+with Cartesian indexing these collections will be tuples with elements of type
+``Union{Int,Colon}``.
+
+When implementing this method it is expected that the returned collection can
+be used to index grid functions to obtain grid functions on the boundary grid.
+"""
+function boundary_indices end
+
+"""
     eval_on(g::Grid, f)
 
 Lazy evaluation of `f` on the grid. `f` can either be on the form `f(x,y,...)`
@@ -100,3 +112,5 @@
 
 _ncomponents(::Type{<:Number}) = 1
 _ncomponents(T::Type{<:SVector}) = length(T)
+
+

--- a/src/Grids/tensor_grid.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/src/Grids/tensor_grid.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -75,7 +75,6 @@
     return LazyTensors.concatenate_tuples(per_grid...)
 end
 
-
 """
     boundary_grid(g::TensorGrid, id::TensorGridBoundary)
 
@@ -87,6 +86,16 @@
     return TensorGrid(new_grids...)
 end
 
+function boundary_indices(g::TensorGrid, id::TensorGridBoundary)
+    per_grid_ind = map(g.grids) do g
+        ntuple(i->:, ndims(g))
+    end
+
+    local_b_ind = boundary_indices(g.grids[grid_id(id)], boundary_id(id))
+    b_ind = Base.setindex(per_grid_ind, local_b_ind, grid_id(id))
+
+    return LazyTensors.concatenate_tuples(b_ind...)
+end
 
 function combined_coordinate_vector_type(coordinate_types...)
     combined_coord_length = mapreduce(_ncomponents, +, coordinate_types)

--- a/src/Grids/zero_dim_grid.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/src/Grids/zero_dim_grid.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -25,3 +25,4 @@
 
 boundary_identifiers(g::ZeroDimGrid) = ()
 boundary_grid(g::ZeroDimGrid, ::Any) = throw(ArgumentError("ZeroDimGrid has no boundaries"))
+boundary_indices(g::ZeroDimGrid, ::Any) = throw(ArgumentError("ZeroDimGrid has no boundaries"))

--- a/test/Grids/equidistant_grid_test.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/test/Grids/equidistant_grid_test.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -68,6 +68,17 @@
         @test boundary_grid(g, Upper()) == ZeroDimGrid(1.0)
     end
 
+    @testset "boundary_indices" begin
+        g = EquidistantGrid(0:0.1:1)
+        @test boundary_indices(g, Lower()) == (1,)
+        @test boundary_indices(g, Upper()) == (11,)
+
+        g = EquidistantGrid(2:0.1:10)
+        @test boundary_indices(g, Lower()) == (1,)
+        @test boundary_indices(g, Upper()) == (81,)
+
+    end
+
     @testset "refine" begin
         g = EquidistantGrid(0:0.1:1)
         @test refine(g, 1) == g

--- a/test/Grids/tensor_grid_test.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/test/Grids/tensor_grid_test.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -169,6 +169,21 @@
         @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂)
         @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄)
     end
+
+    @testset "boundary_indices" begin
+        g₁ = EquidistantGrid(range(0,1,length=11))
+        g₂ = EquidistantGrid(range(2,3,length=6))
+        g₄ = ZeroDimGrid(@SVector[1,2])
+
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Lower}()) == (1,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == (11,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Lower}()) == (:,1)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Upper}()) == (:,6)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Lower}()) == (1,)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == (11,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Lower}()) == (1,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Upper}()) == (11,)
+    end
 end
 
 @testset "combined_coordinate_vector_type" begin

--- a/test/Grids/zero_dim_grid_test.jl	Thu Nov 23 16:08:47 2023 +0100
+++ b/test/Grids/zero_dim_grid_test.jl	Fri Dec 01 11:52:26 2023 +0100
@@ -41,4 +41,8 @@
     @testset "boundary_grid" begin
         @test_throws ArgumentError("ZeroDimGrid has no boundaries") boundary_grid(ZeroDimGrid(@SVector[1.0,2.0]), :bid)
     end
+
+    @testset "boundary_indices" begin
+        @test_throws ArgumentError("ZeroDimGrid has no boundaries") boundary_indices(ZeroDimGrid(@SVector[1.0,2.0]), :bid)
+    end
 end