--- a/Notes.md	Fri Feb 12 16:16:45 2021 +0100
+++ b/Notes.md	Sat Feb 13 16:07:46 2021 +0100
@@ -298,3 +298,16 @@
 component(gf,:,2) # Andra kolumnen av en matris
 @ourview gf[:,:][2]
 ```
+
+## Grids embedded in higher dimensions
+
+For grids generated by asking for boundary grids for a regular grid, it would
+make sense if these grids knew they were embedded in a higher dimension. They
+would return coordinates in the full room. This would make sense when
+drawing points for example, or when evaluating functions on the boundary.
+
+Implementation of this is an issue that requires some thought. Adding an extra
+"Embedded" type for each grid would make it easy to understand each type but
+contribute to "type bloat". On the other hand adapting existing types to
+handle embeddedness would complicate the now very simple grid types. Are there
+other ways of doing the implentation?

--- a/src/Grids/EquidistantGrid.jl	Fri Feb 12 16:16:45 2021 +0100
+++ b/src/Grids/EquidistantGrid.jl	Sat Feb 13 16:07:46 2021 +0100
@@ -1,11 +1,19 @@
 """
-    EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}
+    EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T})
+	EquidistantGrid{T}()
+
+`EquidistantGrid` is a grid with equidistant grid spacing per coordinat direction.
 
-EquidistantGrid is a grid with equidistant grid spacing per coordinat direction.
-The domain is defined through the two points P1 = x̄₁, P2 = x̄₂ by the exterior
-product of the vectors obtained by projecting (x̄₂-x̄₁) onto the coordinate
-directions. E.g for a 2D grid with x̄₁=(-1,0) and x̄₂=(1,2) the domain is defined
-as (-1,1)x(0,2). The side lengths of the grid are not allowed to be negative
+`EquidistantGrid(size, limit_lower, limit_upper)` construct the grid with the
+domain defined by the two points P1, and P2 given by `limit_lower` and
+`limit_upper`. The length of the domain sides are given by the components of
+(P2-P1). E.g for a 2D grid with P1=(-1,0) and P2=(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 `size`.
+
+`EquidistantGrid{T}()` constructs a 0-dimensional grid.
+
 """
 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
     size::NTuple{Dim, Int}
@@ -22,6 +30,9 @@
         end
         return new{Dim,T}(size, limit_lower, limit_upper)
     end
+
+	# Specialized constructor for 0-dimensional grid
+	EquidistantGrid{T}() where T = new{0,T}((),(),())
 end
 export EquidistantGrid
 
@@ -35,9 +46,9 @@
 	return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
 end
 
-function Base.eachindex(grid::EquidistantGrid)
-    CartesianIndices(grid.size)
-end
+Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
+
+Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
 
 Base.size(g::EquidistantGrid) = g.size
 
@@ -104,3 +115,23 @@
 """
 boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),dimension(g)))...)...,)
 export boundary_identifiers
+
+
+"""
+    boundary_grid(grid::EquidistantGrid,id::CartesianBoundary)
+	boundary_grid(::EquidistantGrid{1},::CartesianBoundary{1})
+
+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)
+	dims = collect(1:dimension(grid))
+	orth_dims = dims[dims .!= dim(id)]
+	if orth_dims == dims
+		throw(DomainError("boundary identifier not matching grid"))
+	end
+    return restrict(grid,orth_dims)
+end
+export boundary_grid
+boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()

--- a/src/SbpOperators/volumeops/quadratures/quadrature.jl	Fri Feb 12 16:16:45 2021 +0100
+++ b/src/SbpOperators/volumeops/quadratures/quadrature.jl	Sat Feb 13 16:07:46 2021 +0100
@@ -1,36 +1,29 @@
 """
-    Quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils)
+    quadrature(grid::EquidistantGrid, closure_stencils, inner_stencil)
+    quadrature(grid::EquidistantGrid, closure_stencils)
 
 Creates the quadrature operator `H` as a `TensorMapping`
 
-The quadrature approximates the integral operator on the grid using
+`H` approximiates the integral operator on `grid` the using the stencil
 `inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions.
+for the points in the closure regions. If `inner_stencil` is omitted a central
+interior stencil with weight 1 is used.
 
 On a one-dimensional `grid`, `H` is a `VolumeOperator`. On a multi-dimensional
 `grid`, `H` is the outer product of the 1-dimensional quadrature operators in
 each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details.
+`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
+`H` is a 0-dimensional `IdentityMapping`.
 """
-function Quadrature(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
+function quadrature(grid::EquidistantGrid, closure_stencils, inner_stencil = CenteredStencil(one(eltype(grid))))
     h = spacing(grid)
     H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
-    for i ∈ 2:Dim
+    for i ∈ 2:dimension(grid)
         Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
         H = H∘Hᵢ
     end
     return H
 end
-export Quadrature
-
-"""
-    DiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
+export quadrature
 
-Creates the quadrature operator with the inner stencil 1/h and 1-element sized
-closure stencils (i.e the operator is diagonal)
-"""
-function DiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T}
-    inner_stencil = Stencil(one(T), center=1)
-    return Quadrature(grid, inner_stencil, closure_stencils)
-end
-export DiagonalQuadrature
+quadrature(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()

--- a/test/testGrids.jl	Fri Feb 12 16:16:45 2021 +0100
+++ b/test/testGrids.jl	Sat Feb 13 16:07:46 2021 +0100
@@ -13,9 +13,12 @@
     @test_throws DomainError EquidistantGrid(1,1.0,-1.0)
     @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,))
 
-    # size
-    @test size(EquidistantGrid(4,0.0,1.0)) == (4,)
-    @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+    @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
 
     # dimension
     @test dimension(EquidistantGrid(4,0.0,1.0)) == 1
@@ -63,6 +66,39 @@
         @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
 end
 
 end

--- a/test/testSbpOperators.jl	Fri Feb 12 16:16:45 2021 +0100
+++ b/test/testSbpOperators.jl	Sat Feb 13 16:07:46 2021 +0100
@@ -393,24 +393,31 @@
     end
 end
 
-@testset "DiagonalQuadrature" begin
+@testset "Quadrature diagonal" begin
     Lx = π/2.
     Ly = Float64(π)
+    Lz = 1.
     g_1D = EquidistantGrid(77, 0.0, Lx)
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
     integral(H,v) = sum(H*v)
-    @testset "Constructors" begin
+    @testset "quadrature" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "0D" begin
+            H = quadrature(EquidistantGrid{Float64}(),op.quadratureClosure)
+            @test H == IdentityMapping{Float64}()
+            @test H isa TensorMapping{T,0,0} where T
+        end
         @testset "1D" begin
-            H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+            H = quadrature(g_1D,op.quadratureClosure)
             inner_stencil = CenteredStencil(1.)
-            @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure)
+            @test H == quadrature(g_1D,op.quadratureClosure,inner_stencil)
             @test H isa TensorMapping{T,1,1} where T
         end
-        @testset "1D" begin
-            H = DiagonalQuadrature(g_2D,op.quadratureClosure)
-            H_x = DiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
-            H_y = DiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
+        @testset "2D" begin
+            H = quadrature(g_2D,op.quadratureClosure)
+            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
+            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
             @test H == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
@@ -419,12 +426,12 @@
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+            H = quadrature(g_1D,op.quadratureClosure)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+            H = quadrature(g_2D,op.quadratureClosure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -441,7 +448,7 @@
 
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -450,7 +457,7 @@
 
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -464,13 +471,13 @@
             u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-4
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-8
             end
@@ -524,14 +531,14 @@
             u = evalOn(g_1D,x->x^3-x^2+1)
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
@@ -542,14 +549,14 @@
             u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15