--- a/docs/src/grids_and_grid_functions.md	Wed Jun 26 12:42:28 2024 +0200
+++ b/docs/src/grids_and_grid_functions.md	Sun Jun 30 10:50:44 2024 +0200
@@ -7,6 +7,24 @@
 ## Interface for grids
 All grids are expected to work as a grid function for the coordinate function, and thus implements Julia's Indexing- and Iteration-interfaces. Notably they are *not* abstract arrays because that inteface is too restrictive for the types of grids we wish to implement.
 
+
+## Plotting
+Plotting of grids and grid functions is supported through a package extension with Makie.jl.
+
+For grids we have:
+* `plot(::Grid{<:Any,2})` (same as `lines`)
+* `lines(::Grid{<:Any,2})`
+* `scatter(::Grid{<:Any,2})`
+
+For 1D grid functions we have:
+* `plot(::Grid{<:Any,1}, ::AbstractVector)` (same as `lines`)
+* `lines(::Grid{<:Any,1}, ::AbstractVector)`
+* `scatter(::Grid{<:Any,1}, ::AbstractVector)`
+
+For 2D grid functions we have:
+* `plot(::Grid{<:Any,2}, ::AbstractArray{<:Any,2})` (constructs a 2d mesh)
+* `surface(::Grid{<:Any,2}, ::AbstractArray{<:Any,2})`
+
 ## To write about
 <!-- # TODO: -->
 * Grid functions

--- a/ext/SbplibMakieExt.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/ext/SbplibMakieExt.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -8,9 +8,6 @@
 function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
     ps = map(Tuple, g)[:]
     values = gf[:]
-
-    N = length(ps)
-
     faces = Vector{NTuple{3,Int}}()
 
     n,m = size(g)

--- a/src/Grids/Grids.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/Grids/Grids.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -38,11 +38,13 @@
 export boundary_indices
 export boundary_identifiers
 export boundary_grid
+export min_spacing
 export coarsen
 export refine
 export eval_on
 export componentview
 export ArrayComponentView
+export normal
 
 export BoundaryIdentifier
 export TensorGridBoundary

--- a/src/Grids/equidistant_grid.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -47,6 +47,8 @@
 """
 inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
 
+min_spacing(g::EquidistantGrid) = spacing(g)
+
 
 boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
 boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])

--- a/src/Grids/grid.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/Grids/grid.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -82,6 +82,14 @@
 # TODO: Implement `setindex!`?
 # TODO: Implement a more general ComponentView that can handle non-AbstractArrays.
 
+
+"""
+    min_spacing(g::Grid)
+
+The smallest distance between any pair of grid points in `g`.
+"""
+function min_spacing end
+
 """
     refine(g::Grid, r)
 

--- a/src/Grids/mapped_grid.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -81,3 +81,27 @@
     end
 end
 
+"""
+    normal(g::MappedGrid, boundary)
+
+The outward pointing normal as a grid function on the boundary
+"""
+function normal(g::MappedGrid, boundary)
+    b_indices = boundary_indices(g, boundary)
+    σ =_boundary_sign(component_type(g), boundary)
+    return map(jacobian(g)[b_indices...]) do ∂x∂ξ
+        ∂ξ∂x = inv(∂x∂ξ)
+        k = grid_id(boundary)
+        σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
+    end
+end
+
+function _boundary_sign(T, boundary)
+    if boundary_id(boundary) == Upper()
+        return one(T)
+    elseif boundary_id(boundary) == Lower()
+        return -one(T)
+    else
+        throw(ArgumentError("The boundary identifier must be either `Lower()` or `Upper()`"))
+    end
+end

--- a/src/Grids/tensor_grid.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/Grids/tensor_grid.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -50,6 +50,12 @@
 Base.size(g::TensorGrid, d) = size(g)[d]
 
 
+function min_spacing(g::TensorGrid)
+    relevant_grids = filter(g->!isa(g,ZeroDimGrid),g.grids)
+    d = min_spacing.(relevant_grids)
+    return minimum(d)
+end
+
 refine(g::TensorGrid, r::Int) = mapreduce(g->refine(g,r), TensorGrid, g.grids)
 coarsen(g::TensorGrid, r::Int) = mapreduce(g->coarsen(g,r), TensorGrid, g.grids)
 

--- a/src/RegionIndices/RegionIndices.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/src/RegionIndices/RegionIndices.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -13,7 +13,7 @@
     Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
     Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
     Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i)
-    Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed?
+    Index(t::Tuple{T, DataType}) where T<:Integer = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed?
 end
 
 export Index

--- a/test/Grids/equidistant_grid_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -57,6 +57,11 @@
         @test inverse_spacing(EquidistantGrid(0:0.1:10)) == 10
     end
 
+    @testset "min_spacing" begin
+        @test min_spacing(EquidistantGrid(0:10)) == 1
+        @test min_spacing(EquidistantGrid(0:0.1:10)) == 0.1
+    end
+
     @testset "boundary_identifiers" begin
         g = EquidistantGrid(0:0.1:10)
         @test boundary_identifiers(g) == (Lower(), Upper())

--- a/test/Grids/mapped_grid_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/Grids/mapped_grid_test.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -183,4 +183,15 @@
     lg = equidistant_grid((0,0), (1,1), 10, 11)
     @test logicalgrid(mg) == lg
     @test collect(mg) == map(x̄, lg)
+
+
+    @testset "normal" begin
+        @test normal(mg, CartesianBoundary{1,Lower}()) == fill(@SVector[-1,0], 11)
+        @test normal(mg, CartesianBoundary{1,Upper}()) == fill(@SVector[1,0], 11)
+        @test normal(mg, CartesianBoundary{2,Lower}()) == fill(@SVector[0,-1], 10)
+        @test normal(mg, CartesianBoundary{2,Upper}()) ≈ map(boundary_grid(mg,CartesianBoundary{2,Upper}())|>logicalgrid) do ξ̄
+            α = 1-2ξ̄[1]
+            @SVector[α,1]/√(α^2 + 1)
+        end
+    end
 end

--- a/test/Grids/tensor_grid_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/Grids/tensor_grid_test.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -138,6 +138,15 @@
         @test axes(g) == (1:11,1:6)
     end
 
+    @testset "min_spacing" begin
+        g₁ = EquidistantGrid(range(0,1,length=11))
+        g₂ = EquidistantGrid(range(2,3,length=6))
+        g₃ = ZeroDimGrid(@SVector[1,2])
+
+        @test min_spacing(TensorGrid(g₁, g₂)) == 1/10
+        @test min_spacing(TensorGrid(g₂, g₃)) == 1/5
+    end
+
     @testset "refine" begin
         g1(n) = EquidistantGrid(range(0,1,length=n))
         g2(n) = EquidistantGrid(range(2,3,length=n))

--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -4,13 +4,13 @@
 
 using Tullio
 
-struct DummyMapping{T,R,D} <: LazyTensor{T,R,D} end
+struct TransposableDummyMapping{T,R,D} <: LazyTensor{T,R,D} end
 
-LazyTensors.apply(m::DummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply
-LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose
+LazyTensors.apply(m::TransposableDummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply
+LazyTensors.apply_transpose(m::TransposableDummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose
 
-LazyTensors.range_size(m::DummyMapping) = :range_size
-LazyTensors.domain_size(m::DummyMapping) = :domain_size
+LazyTensors.range_size(m::TransposableDummyMapping) = :range_size
+LazyTensors.domain_size(m::TransposableDummyMapping) = :domain_size
 
 
 struct SizeDoublingMapping{T,R,D} <: LazyTensor{T,R,D}
@@ -24,7 +24,7 @@
 
 
 @testset "Mapping transpose" begin
-    m = DummyMapping{Float64,2,3}()
+    m = TransposableDummyMapping{Float64,2,3}()
     @test m' isa LazyTensor{Float64, 3,2}
     @test m'' == m
     @test apply(m',zeros(Float64,(0,0)), 0, 0, 0) == :apply_transpose

--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -13,20 +13,9 @@
 using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
 using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
 
-"""
-    monomial(x,k)
-
-Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
-Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
-"""
-function monomial(x,k)
-    if k < 0
-        return zero(x)
-    end
-    x^k/factorial(k)
-end
 
 @testset "undivided_skewed04" begin
+    monomial(x,k) = k < 0 ? zero(x) : x^k/factorial(k)
     g = equidistant_grid(0., 11., 20)
     D,Dᵀ = undivided_skewed04(g, 1)
 

--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Wed Jun 26 12:42:28 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Sun Jun 30 10:50:44 2024 +0200
@@ -7,19 +7,6 @@
 
 using Sbplib.SbpOperators: closure_size, Stencil, VolumeOperator
 
-"""
-    monomial(x,k)
-
-Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
-Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
-"""
-function monomial(x,k)
-    if k < 0
-        return zero(x)
-    end
-    x^k/factorial(k)
-end
-
 @testset "first_derivative" begin
     @testset "Constructors" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
@@ -39,6 +26,8 @@
     @testset "Accuracy conditions" begin
         N = 20
         g = equidistant_grid(0//1, 2//1, N)
+
+        monomial(x,k) = k < 0 ? zero(x) : x^k/factorial(k)
         @testset for order ∈ [2,4]
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
             D₁ = first_derivative(g, stencil_set)