diff -r 6dfa8001eed2 -r f8be66557e0f src/PlanarFE.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/PlanarFE.jl	Wed Dec 22 11:13:38 2021 +0200
@@ -0,0 +1,937 @@
+"""
+Simple planar (2D) finite element discretisation on a box.
+"""
+module PlanarFE
+
+import ..Loops: Box
+import ..LinSolve: linsolve
+using ..Metaprogramming
+
+export FEModel,
+       FEUniform,
+       UniformP0,
+       PlanarP0,
+       LineP0,
+       UniformP1,
+       PlanarP1,
+       LineP1,
+       UniformP2,
+       PlanarP2,
+       LineP2,
+       simplices,
+       nodes,
+       interpolate,
+       differential,
+       minimise
+
+const Location{N} = NTuple{N, Float64}
+abstract type FEModel{N} end
+abstract type FEUniform{N} <: FEModel{N} end
+
+const Index = Int
+const CaseIndex = Int
+
+const NodeId{N} = NTuple{N, Index}
+
+struct Node{N}
+    id :: NodeId{N}
+    x :: Location{N}
+end
+
+@inline function Base.getindex(a :: Array{T,N}, n :: Node{N}) where {T,N}
+    return getindex(a, n.id...)
+end
+
+@inline function Base.setindex!(a :: Array{T,N}, x :: T, n :: Node{N}) where {T,N}
+    return setindex!(a, x, n.id...)
+end
+
+##############################################
+# Simplices
+##############################################
+
+# (D-1)-dimensional simplex in ℝ^N. D is the node count.
+struct Simplex{N,D,Id}
+    id :: Id
+    nodes :: NTuple{D, Node{N}}
+    @inline Simplex{Id}(id :: Id, nodes :: Vararg{Node{N}, D}) where {N,D,Id} = new{N,D,Id}(id, nodes)
+    @inline Simplex{N,D,Id}(id :: Id, nodes :: Vararg{Node{N}, D}) where {N,D,Id} = new{N,D,Id}(id, nodes)
+end
+
+# Uniformly indexed simplices
+const UCubeId{N} = NTuple{N, Int}
+const USimplexId{N} = Tuple{UCubeId{N}, CaseIndex}
+const USimplex{N,D} = Simplex{N,D,USimplexId{N}}
+const PlanarSimplex = USimplex{2,3}
+const RealInterval = USimplex{1,2}
+
+@inline function Base.getindex(a :: Array{T,Nplus1}, s :: USimplex{N,D}) where {T,Nplus1,D,N}
+    idx = (s.id[1]..., s.id[2]) :: NTuple{Nplus1,Int}
+    return getindex(a, idx...)
+end
+
+@inline function Base.setindex!(a :: Array{T,Nplus1}, x :: T, s :: USimplex{N,D}) where {T,Nplus1,D,N}
+    idx = (s.id[1]..., s.id[2]) :: NTuple{Nplus1,Int}
+    return setindex!(a, x, idx...)
+end
+
+# We do this by invididual cases to allow typing return values
+# (N+1 is not calculable at compile time)
+@inline function nodes(s :: Simplex{N,D,Id}) where {N,D,Id}
+    return s.nodes
+end
+
+@inline @generated function dot(x :: Location{N}, y :: Location{N}) where N
+    return foldl_exprs( :( + ), :( x[$i]*y[$i] ) for i=1:N)
+end
+
+@inline function orthog((x, y) :: Location{2})
+    return (y, -x)
+end
+
+@inline function in_simplex(t :: RealInterval, x :: Location{1})
+    (x₀, x₁) = nodes(t)
+    return ((x .- x₀)*(x₁ .- x₀) ≥ 0 &&
+            (x .- x₁)*(x₀ .- x₁) ≥ 0)
+end
+
+@inline function in_simplex(t :: PlanarSimplex, x :: Location{2})
+    (x₀, x₁, x₁) = nodes(t)
+    return (dot(x .- x₀, orthog(x₁ .- x₀)) ≥ 0 &&
+            dot(x .- x₁, orthog(x₂ .- x₁)) ≥ 0 &&
+            dot(x .- x₂, orthog(x₀ .- x₂)) ≥ 0)
+end
+
+
+##############################################
+# Simple planar grid and triangulation
+##############################################
+
+@inline function cellwidth(box :: Box{N, Float64}, gridsize :: Dims{N}) where N
+    return (((box[i][2]-box[i][1])/(gridsize[i]-1) for i=1:N)...,)
+end
+
+# Cannot calculate N+1 compile-time
+struct UniformGrid{N,Nplus1}
+    box :: Box{N, Float64}
+    size :: Dims{N}
+    cellwidth :: NTuple{N, Float64}
+
+    function UniformGrid{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+        # TODO: Nested tuple access is inefficient
+        @assert(N+1 == Nplus1)
+        @assert(all(box[i][2]>box[i][1] for i=1:N))
+        @assert(all(gridsize .≥ 2))
+        return new{N,Nplus1}(box, gridsize, cellwidth(box, gridsize))
+    end
+end
+
+const LineGrid = UniformGrid{1,2}
+const PlanarGrid = UniformGrid{2,3}
+
+# Number of simplices in a cube of dimension N.
+n_simplex_cases( :: LineGrid) = 1 :: CaseIndex
+n_simplex_cases( :: PlanarGrid) = 2 :: CaseIndex
+max_cube_id(grid :: UniformGrid{N, Nplus1}) where {N, Nplus1} = grid.size .- 1
+max_simplex_id(grid :: UniformGrid{N, Nplus1}) where {N, Nplus1} =
+    (max_cube_id(grid)..., n_simplex_cases(grid))
+
+∈(x :: Location{N}, b :: Box{N}) where N = all(b[i][1] ≤ x[i] ≤ b[i][2] for i=1:N)
+∉(x :: Location{N}, b :: Box{N}) where N = !(x ∈ b)
+∈(x :: Location{N}, grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} = x ∈ grid.box
+∉(x :: Location{N}, grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} = x ∉ grid.box
+
+@inline function width(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1}
+    return ((grid.box[i][2]-grid.box[i][1] for i=1:N)...,)
+end
+
+@inline function shift(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1}
+    return ((grid.box[i][1] for i=1:N)...,)
+end
+
+function construct_node(grid :: UniformGrid{N,Nplus1}, idx :: NodeId{N}) :: Node{N} where {N,Nplus1}
+    loc = grid.cellwidth.*(idx.-1).+shift(grid)
+    return Node{N}(idx, loc)
+end
+
+@inline function construct_simplex(grid :: LineGrid, id :: USimplexId{1}) :: RealInterval
+    # TODO: optimise case away
+    (idx, _) = id
+    return RealInterval(construct_node(grid, idx), construct_node(grid, idx.+(1,)))
+end
+
+
+# In 2D we alternate ⧄ and ⧅. The “base” 90° corner is always listed first.
+# It has importance for P0 elements. Otherwise arbitrary we proceed clockwise.
+# Julia is real GARBAGE. It does tons of memory alloc trying to access this in
+# any form, as Array or Tuple. Therefore things are in code below.
+# offset₂ = [[
+#     [(1, 0), (0, 0), (1, 1)], # ◪
+#     [(0, 1), (1, 1), (0, 0)], # ◩
+# ],[
+#     [(0, 0), (1, 0), (0, 1)], # ⬕
+#     [(1, 1), (1, 0), (0, 1)], # ⬔
+# ]] :: Array{Array{Array{Tuple{Int,Int},1},1},1}
+
+
+@inline function sqform₂(idx :: UCubeId{2})
+    return 1+(idx[1]+idx[2])%2
+end
+
+# 2D simplices are identified by their “leading” node, i.e., the lower-left corner of the grid square
+# that contains them, and the “case” number 1 or 2.
+function construct_simplex(grid :: PlanarGrid, simplexid :: USimplexId{2}, form :: Index) :: PlanarSimplex
+    (idx, case) = simplexid
+    @assert(all(idx .< grid.size))
+    @assert(form==1 || form==2)
+    @assert(case==1 || case==2)
+    (off1 :: NodeId{2}, off2 :: NodeId{2}, off3 :: NodeId{2}) = (form==1 ?
+            (case==1 ?
+                ((1, 0), (0, 0), (1, 1)) # ◪
+            :
+                ((0, 1), (1, 1), (0, 0)) # ◩
+        ) : (case==1 ?
+                ((0, 0), (1, 0), (0, 1)) # ⬕
+            :
+                ((1, 1), (1, 0), (0, 1)) # ⬔
+        ))
+
+    return PlanarSimplex(simplexid,
+                         construct_node(grid, idx .+ off1),
+                         construct_node(grid, idx .+ off2),
+                         construct_node(grid, idx .+ off3))
+end
+
+function construct_simplex(grid :: PlanarGrid, simplexid :: USimplexId{2}) :: PlanarSimplex
+    (idx, _) = simplexid
+    return construct_simplex(grid, simplexid, sqform₂(idx))
+end
+
+
+@inline function get_cubeid(grid :: UniformGrid{N,Nplus1}, x :: Location{N}) :: UCubeId{N} where {N,Nplus1}
+    ((
+        # Without the :: Int typeassert this generates lots of memory allocs, although a typeof()
+        # check here confirms that the type is Int. Julia is just as garbage as everything.
+        max(min(ceil(Index, (x[i]-grid.box[i][1])/grid.cellwidth[i]), grid.size[i] - 1), 1) :: Int
+    for i=1:N)...,)
+end
+
+
+@inline function get_offset(grid :: UniformGrid{N,Nplus1}, cubeid :: UCubeId{N}, x :: Location{N}) :: Location{N} where {N,Nplus1}
+    ((
+        (x[i] - (cubeid[i] - 1) * grid.cellwidth[i])
+    for i=1:N)...,)
+end
+
+@inline function simplex_at(grid :: LineGrid, x :: Location{1}) :: RealInterval
+    if x ∉ grid.box
+        throw("Coordinate $x out of domain $(grid.box)")
+    else
+        return construct_simplex(grid, get_cubeid(grid, x))
+    end
+end
+
+@inline function simplex_at(grid :: PlanarGrid, x :: Location{2}) :: PlanarSimplex
+    if x ∉ grid.box
+        throw("Coordinate $x out of domain $(grid.box)")
+    else
+        # Lower left node is is cube id. Therefore maximum cubeid is grid size - 1
+        # in each dimension.
+        cubeid = get_cubeid(grid, x)
+        off = get_offset(grid, cubeid, x)
+        form = sqform₂(cubeid)
+        if form==0
+            # ⧄ → ◪ or ◩
+            case = (off[1]>off[2] ? 1 : 2 ) :: CaseIndex
+        else
+            # ⧅ → ⬕ or ⬔
+            case = (off[1]+off[2]<1 ? 1 : 2) :: CaseIndex
+        end
+        return construct_simplex(grid, (cubeid, case), form)
+    end
+end
+
+##############################################
+# Edges for P2 elements
+##############################################
+
+const Edge{N,Id} = Simplex{N,2,Id}
+const UEdgeId{N} = USimplexId{N}
+const UEdge{N} = Edge{N,UEdgeId{N}}
+
+@inline midpoint(e :: Edge{N,Id}) where {N,Id} = (e.nodes[1].x .+ e.nodes[2].x) ./ 2.0
+
+# For our purposes, an edge is something useful for putting an interpolation midpoint on.
+# This in 1D there's one “edge”, which is just the original simplex.
+n_cube_edges( :: LineGrid) = 1
+n_cube_edges( :: PlanarGrid) = 3
+max_edge_id(grid :: UniformGrid{N, Nplus1}) where {N, Nplus1} = (grid.size..., n_cube_edges(grid))
+
+@inline function edge_id(node₁ :: NodeId{1}, node₂ :: NodeId{1}) :: UEdgeId{1}
+    return (min(node₁[1], node₂[1]),)
+end
+
+# Edges are identified by the lower left node of the “least-numbered” cube that
+# contains them, along with 1=horizontal, 2=vertical, 3=diagonal.
+# Both of the alternating cases of diagonal edges have the same number.
+# Since the least-numbered cube identifies edges, only both diagonal the left and
+# bottom border of the cube are counted for that cube, so there are 3 possible
+# edges counted towards a cube.
+function edge_id(idx₁ :: NodeId{2}, idx₂ :: NodeId{2}) :: UEdgeId{2}
+    # For performance this does not verify that idx₁ and idx₂ are in the
+    # same simplex, it only uses their ordering!
+    if idx₂[1] > idx₁[1] # one of ↗︎, →, ↘︎
+        if idx₂[2] > idx₁[2]
+            return (idx₁, 3) # diagonal ↗︎ in cube(idx₁)
+        elseif idx₂[2] == idx₁[2]
+            return (idx₁, 1) # lower side → in cube(idx₁)
+        else # idx₂[2] < idx₁[2]
+            return ((idx₁[1], idx₂[2]), 3) # ↘︎ in cube(prev(idx₁,idx₂))
+        end
+    elseif idx₂[1] == idx₁[1] # one of ↑, ↓
+        if idx₂[2] > idx₁[2]
+            return (idx₁, 2) # ↑ in cube(idx₁)
+        elseif idx₂[2] == idx₁[2]
+            throw("Edge cannot end in starting idx.")
+        else # idx₂[2] < idx₁[2]
+            return (idx₂, 2) # ↓ cube(idx₂)
+        end
+    else # idx₂[1] < idx₁[1] # one of ↙︎, ←, ↖︎
+        if idx₂[2] > idx₁[2]
+            return ((idx₁[1], idx₂[2]), 3) # diagonal ↖︎ in cube(prev(idx₁,idx₂))
+        elseif idx₂[2] == idx₁[2]
+            return (idx₂, 1) # lower side ← in cube(idx₂)
+        else # idx₂[2] < idx₁[2]
+            return (idx₂, 3) # ↙︎ in cube(idx₂)
+        end
+    end
+end
+
+function construct_edge(n₁ :: Node{N}, n₂ :: Node{N}) where N
+    return UEdge{N}(edge_id(n₁.id, n₂.id), n₁, n₂)
+end
+
+##############################################
+# Iteration over simplices and nodes
+##############################################
+
+struct Iter{G,W}
+    grid :: G
+end
+
+# @inline function nodes(s :: PlanarSimplex{N}) where N
+#     return Iter{PlanarSimplex{N}, :nodes}(s)
+# end
+
+# @inline function Base.iterate(sn :: Iter{PlanarSimplex{N}, :nodes}) where N
+#     return (sn.p.base, Index(0))
+# end
+
+# @inline function Base.iterate(sn :: Iter{PlanarSimplex{N}, :nodes}, idx :: Index) where N
+#     if idx==length(sn.p.tips)
+#         return nothing
+#     else
+#         idx=idx+1
+#         return (sn.p.tips[idx], idx)
+#     end
+# end
+
+@inline function nodes(p :: M) where {N, M <: FEUniform{N}}
+    return nodes(p.grid)
+end
+
+@inline function nodes(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1}
+    return Iter{UniformGrid{N,Nplus1}, Node{N}}(grid)
+end
+
+@inline function first_nodeid( :: UniformGrid{N,Nplus1}) :: NodeId{N} where {N,Nplus1}
+    return (( 1 :: Int for i=1:N)...,)
+end
+
+# Need to hard-code this for each N to avoid Julia munching memory.
+# Neither @generated nor plain code works.
+@inline function next_nodeid(grid :: LineGrid, idx :: NodeId{1}) :: Union{Nothing, NodeId{1}}
+    if idx[1] < grid.size[1]
+        return (idx[1]+1,)
+    else
+        return nothing
+    end
+end
+
+@inline function next_nodeid(grid :: PlanarGrid, idx :: NodeId{2}) :: Union{Nothing, NodeId{2}}
+    if idx[1] < grid.size[1]
+        return (idx[1]+1, idx[2])
+    elseif idx[2] < grid.size[2]
+        return (1 :: Int, idx[2]+1)
+    else
+        return nothing
+    end
+end
+
+function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Node{N}}) ::
+         Union{Nothing, Tuple{Node{N}, NodeId{N}}} where {N,Nplus1}
+    idx = first_nodeid(iter.grid)
+    return construct_node(iter.grid, idx), idx
+end
+
+function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Node{N}}, idx :: NodeId{N}) ::
+         Union{Nothing, Tuple{Node{N}, NodeId{N}}} where {N,Nplus1}
+    nidx = next_nodeid(iter.grid, idx)
+    if isnothing(nidx)
+        return nothing
+    else
+        return construct_node(iter.grid, nidx), nidx
+    end
+end
+
+# @generated function Base.iterate(iter :: Iter{N, Node{N}}, idx :: NodeId{N}) where N
+#     # This generation reduces memory allocations
+#     do_coord(i) = :(
+#         if idx[$i] < iter.grid.size[$i]
+#             nidx = $(lift_exprs(( (:( 1 :: Int ) for _=1:i-1)..., :( idx[$i]+1 ), ( :( idx[$j] ) for j=i+1:N)...,))) :: NTuple{N, Int}
+#             return construct_node(iter.grid, nidx), nidx
+#         end
+#     )
+
+#     return quote
+#         $(sequence_exprs( do_coord(i) for i=1:N ))
+#         return nothing
+#     end
+# end
+
+# function Base.iterate(iter :: Iter{N, Node{N}}, idx :: NodeId{N}) where N
+#     for i=1:N
+#         if idx[i]<iter.grid.size[i]
+#             nidx = (( 1 :: Int for _=1:i-1)..., idx[i]+1, idx[i+1:N]...)
+#             return construct_node(iter.grid, nidx), nidx
+#         end
+#     end
+
+#     return nothing
+# end
+
+@inline function simplices(p :: M) where {N, M <: FEUniform{N}}
+    return simplices(p.grid)
+end
+
+@inline function simplices(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1}
+    return Iter{UniformGrid{N,Nplus1}, USimplex{N,D}}(grid)
+end
+
+@inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Simplex{N,Nplus1}}) ::
+                 Union{Nothing, Tuple{USimplex{N,Nplus1}, USimplexId{N}}} where {N,Nplus1}
+    cubeid = ((Index(1) for i=1:N)...,)
+    simplexid = (cubeid, Index(1))
+    return (construct_simplex(iter.grid, simplexid), simplexid)
+end
+
+@inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Simplex{N,Nplus1}},
+                              (cubeid, case) :: USimplexId{N}) ::
+                 Union{Nothing, Tuple{USimplex{N,Nplus1}, USimplexId{N}}} where {N,Nplus1}
+    if case < n_simplex_cases(iter.grid)
+        simplexid = (cubeid, case+1)
+        return (construct_simplex(iter.grid, simplexid), simplexid)
+    else
+        for i=1:N
+            if cubeid[i]+1<iter.grid.size[i]
+                new_cubeid=((Index(1) for _=1:i-1)..., cubeid[i]+1, cubeid[i+1:N]...)
+                simplexid = (new_cubeid, Index(1))
+                return (construct_simplex(iter.grid, simplexid), simplexid)
+            end
+        end
+    end
+
+    return nothing
+end
+
+# TODO: wrong for far-side boundary? Need to do by cumbe numbers
+@inline function edges(s :: RealInterval) :: Tuple{UEdge{1}}
+    return (construct_edge(s.nodes...),)
+end
+
+# TODO: Edge -> Simplex
+@inline function edges(s :: PlanarSimplex) :: NTuple{3,UEdge{2}}
+    return (construct_edge(s.nodes[1], s.nodes[2]),
+            construct_edge(s.nodes[1], s.nodes[3]),
+            construct_edge(s.nodes[2], s.nodes[3]))
+end
+
+@inline function edges(p :: M) where {N, M <: FEUniform{N}}
+    return edges(p.grid)
+end
+
+@inline function edges(grid :: LineGrid)
+    return Iter{LineGrid, UEdge{1}}(grid)
+end
+
+@inline function edges(grid :: PlanarGrid)
+    return Iter{PlanarGrid, UEdge{2}}(grid)
+end
+
+@inline function UEdge(grid :: LineGrid, edgeid :: UEdgeId{1}) :: UEdge{1}
+    (cubeid, case) = edgeid
+    @assert(case==1)
+    return UEdge(edgeid, construct_node(grid, cubeid), construct_node(grid, cubeid.+(1,)))
+end
+
+@inline function next_edge_case(grid :: LineGrid, (cubeid, case) :: UEdgeId{1}) :: CaseIndex
+    return -1
+end
+
+@inline function first_edge_case(grid :: LineGrid, cubeid :: UCubeId{1}) :: CaseIndex
+    return cubeid[1] < grid.size[1] ? 1 : -1
+end
+
+@inline function UEdge(grid :: PlanarGrid, edgeid :: UEdgeId{2}) :: UEdge{2}
+    (cubeid, case) = edgeid
+    e(idx1, idx2) = UEdge{2}(edgeid, construct_node(grid, idx1), construct_node(grid, idx2))
+    if case==1
+        return e(cubeid, cubeid.+(1,0))
+    elseif case==2
+        return e(cubeid, cubeid.+(0,1))
+    elseif case==3
+        if sqform₂(cubeid)==1
+            return e(cubeid, cubeid.+(1,1))
+        else
+            return e(cubeid.+(1,0), cubeid.+(0,1))
+        end
+    else
+        throw("Invalid edge")
+    end
+end
+
+@inline function first_edge_case(grid :: PlanarGrid, cubeid :: UCubeId{2}) :: CaseIndex
+    if  cubeid[1]<grid.size[1]
+        return 1
+    elseif cubeid[2]<grid.size[2]
+        return 2
+    else
+        return -1
+    end
+end
+
+
+@inline function next_edge_case(grid :: PlanarGrid, (cubeid, case) :: UEdgeId{2}) :: CaseIndex
+    if case<3 && cubeid[1]<grid.size[1] && cubeid[2]<grid.size[2]
+        return case + 1
+    else
+        return -1
+    end
+end
+
+
+@inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, UEdge{N}}) ::
+                 Union{Nothing, Tuple{UEdge{N}, UEdgeId{N}}} where {N,Nplus1}
+    cubeid = first_nodeid(iter.grid)
+    case = first_edge_case(iter.grid, cubeid)
+    @assert(case ≥ 1)
+    edgeid = (cubeid, case)
+    return (UEdge(iter.grid, edgeid), edgeid)
+end
+
+@inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, UEdge{N}}, (cubeid, case) :: UEdgeId{N}) ::
+                 Union{Nothing, Tuple{UEdge{N}, UEdgeId{N}}} where {N,Nplus1}
+    case = next_edge_case(iter.grid, (cubeid, case))
+    if case ≥ 1
+        edgeid = (cubeid, case)
+        return (UEdge(iter.grid, edgeid), edgeid)
+    else
+        while true
+            cubeid = next_nodeid(iter.grid, cubeid)
+            if isnothing(cubeid)
+                return nothing
+            end
+            case = first_edge_case(iter.grid, cubeid) :: Index
+            if case ≥ 1
+                edgeid = (cubeid, case) :: UEdgeId{N}
+                return (UEdge(iter.grid, edgeid), edgeid)
+            end
+        end
+    end
+
+    return nothing
+end
+
+
+##############################################
+# Evaluation, differentiation, minimisation
+##############################################
+
+function interpolate(p :: M, x :: Location{N}) where {N, M <: FEModel{N}}
+    return interpolate(p, simplex_at(p.grid, x), x)
+end
+
+# Apparent Julia bug / incomplete implementation:
+# https://stackoverflow.com/questions/46063872/parametric-functors-in-julia
+# How to do individually.
+# function (p :: M)(x :: Location{N}) where {N, M <: FEModel{N}}
+#     return interpolate(p, x)
+# end
+
+function differential(p :: M, x :: Location{N}) where {N, M <: FEModel{N}}
+    differential(p, simplex_at(p.grid, x), x)
+end
+
+function minimise(p :: M) where {N, M <: FEModel{N}}
+    minloc = nothing
+    minv = nothing
+    for simplex ∈ simplices(p)
+        thisminloc, thisminv = minimise(p, simplex)
+        if isnothing(minv) || minv < thisminv
+            minv=thisminv
+            minloc=thisminloc
+        end
+    end
+    return minv, minloc
+end
+
+##############################################
+# P0 elements
+##############################################
+
+struct UniformP0{N,Nplus1} <: FEUniform{N}
+    simplex_values :: Array{Float64, Nplus1}
+    grid :: UniformGrid{N,Nplus1}
+
+    function UniformP0{N,Nplus1}(simplex_values :: Array{Float64, Nplus1},
+                                 box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+        grid = UniformGrid{N,Nplus1}(box, gridsize)
+        return new{N,Nplus1}(simplex_values, grid)
+    end
+
+    function UniformP0{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+        grid = UniformGrid{N,Nplus1}(box, gridsize)
+        a = Array{Float64, Nplus1}(undef, max_simplex_id(grid))
+        return new{N,Nplus1}(simplex_values, grid)
+    end
+
+end
+
+const LineP0 = UniformP0{1,1}
+const PlanarP0 = UniformP0{2,3}
+
+# f is not specified a type to allow callable objects and functions
+function UniformP0{N,Nplus1}(f, box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+    p = UniformP0{N,Nplus1}(box, gridsize)
+    discretise!(f, p)
+    return p
+end
+
+function (p :: UniformP0{N,Nplus1})(x :: Location{N}) where {N,Nplus1}
+    return interpolate(p, x)
+end
+
+@inline function interpolate(p :: UniformP0{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1}
+    return p.simplex_values[s]
+end
+
+function differential(p :: UniformP0{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1}
+    return ((0.0 for i=1:N)...,)
+end
+
+@inline function minimise(p :: UniformP0{N,Nplus1}, s :: Simplex{N,Nplus1}) where {N,Nplus1}
+    return s.nodes[1].x, p.simplex_values[s]
+end
+
+function discretise!(f, p :: UniformP0{N,Nplus1}) where {N,Nplus1}
+    for s ∈ simplices(p)
+        p.simplex_values[s]=f(s.nodes[1].x)
+    end
+end
+
+##############################################
+# P1 elements
+##############################################
+
+struct UniformP1{N,Nplus1} <: FEUniform{N}
+    node_values :: Array{Float64, N}
+    grid :: UniformGrid{N,Nplus1}
+
+    function UniformP1{N,Nplus1}(node_values :: Array{Float64, N},
+                                 box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+        grid = UniformGrid{N,Nplus1}(box, gridsize)
+        return new{N,Nplus1}(node_values, grid)
+    end
+
+    function UniformP1{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+        return UniformP1{N,Nplus1}(Array{Float64, N}(undef, gridsize), box, gridsize)
+    end
+end
+
+const LineP1 = UniformP1{1,1}
+const PlanarP1 = UniformP1{2,3}
+
+# f is not specified a type to allow callable objects and functions
+function UniformP1{N,Nplus1}(f, box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+    p = UniformP1{N,Nplus1}(box, gridsize)
+    discretise!(f, p)
+    return p
+end
+
+const P1NodalData{N} = Tuple{Location{N}, Float64}
+
+# 1D model functions
+
+@inline function p1_simplicialf(((x₀,), v₀) :: P1NodalData{1},
+                                ((x₁,), v₁) :: P1NodalData{1}) :: NTuple{2,Float64}
+    a₁ = (v₁-v₀)/(x₁-x₀)
+    a₀ = v₀ - a₁*x₀
+    return a₀, a₁
+end
+
+function p1_modelf(p :: UniformP1{1}, s :: RealInterval)
+    n₀, n₁ = nodes(s)
+    return p1_simplicialf((n₀.x, p.node_values[n₀]),
+                          (n₁.x, p.node_values[n₁]))
+end
+
+
+# 2D model functions
+
+@inline function p1_simplicialf((x₀, v₀) :: P1NodalData{2},
+                                (x₁, v₁) :: P1NodalData{2},
+                                (x₂, v₂) :: P1NodalData{2}) :: NTuple{3,Float64}
+    d = x₁ .- x₀
+    q = x₂ .- x₀
+    r = q[1]*d[2] - d[1]*q[2]
+    a₁ = ((v₂-v₀)*d[2] - (v₁-v₀)*q[2])/r
+    a₂ = ((v₁-v₀)*q[1] - (v₂-v₀)*d[1])/r
+    a₀ = v₀-a₁*x₀[1]-a₂*x₀[2]
+    return a₀, a₁, a₂
+end
+
+@inline function p1_modelf(p :: UniformP1{N,Nplus1}, s :: Simplex{N,Nplus1}) where {N,Nplus1}
+    return  p1_simplicialf(((n.x, p.node_values[n]) for n ∈ nodes(s))...)
+end
+
+# @generated function p1_modelf(p :: UniformP1{N}, s :: PlanarSimplex{N}) where N
+#     return quote
+#         ns = nodes(s)
+#         p1_simplicialf($(lift_exprs( :((ns[$i].x, p.node_values[ns[$i].id...])) for i=1:3 ))...)
+#     end
+# end
+
+#@inline function g(a, (i,))
+#    return a[i]
+#end
+
+# @inline function g(a, (i, j))
+#     return a[i,j]
+# end
+
+# function p1_modelf(p :: UniformP1{2}, s :: PlanarSimplex{2})
+#     n₀, n₁, n₂ = nodes(s)
+#     a, b = n₀.id[1], n₀.id[2]
+#     v₀ = p.node_values[a, b]
+#     a, b = n₁.id[1], n₁.id[2]
+#     v₁ = p.node_values[a, b]
+#     a, b = n₂.id[1], n₂.id[2]
+#     v₂ = p.node_values[a, b]
+#     return p1_simplicialf(n₀.x, v₀,
+#                           n₁.x, v₁,
+#                           n₂.x, v₂)
+# end
+
+
+# Common
+
+function (p :: UniformP1{N,Nplus1})(x :: Location{N}) where {N,Nplus1}
+    return interpolate(p, x)
+end
+
+function interpolate(p :: UniformP1{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1}
+    a = p1_modelf(p, s)
+    return +((a.*(1, x...))...)
+end
+
+@inline function differential(p :: UniformP1{N,Nplus1}, s :: Simplex{N,Nplus1}, :: Location{N}) :: Gradient{N} where {N,Nplus1}
+    a = p1_modelf(p, s)
+    return a[2:end]
+end
+
+@inline function minimise(p :: UniformP1{N,Nplus1}, t :: Simplex{N,Nplus1}) where {N,Nplus1}
+    bestnode = s.nodes[1]
+    bestvalue = p.node_values[bestnode]
+    # The minima is found at a node, so only need to go through them
+    for node ∈ s.nodes[2:end]
+        value = p.node_values[node]
+        if value < bestvalue
+            bestvalue = value
+            bestnode = node
+        end
+    end
+    return bestnode.x, bestvalue
+end
+
+function discretise!(f, p :: UniformP1{N,Nplus1}) where {N,Nplus1}
+    for n ∈ nodes(p)
+        p.node_values[n]=f(n.x)
+    end
+end
+
+##############################################
+# P2 elements
+##############################################
+
+const Gradient{N} = NTuple{N, Float64}
+const P2NodalData{N} = Tuple{Location{N}, Float64, Gradient{N}}
+const LocationAndGradient{N} = Tuple{Location{N}, Gradient{N}}
+
+struct UniformP2{N,Nplus1} <: FEUniform{N}
+    node_values :: Array{Float64, N}
+    edge_values :: Array{Float64, Nplus1}
+    grid :: UniformGrid{N,Nplus1}
+
+    function UniformP2{N,Nplus1}(
+                (node_values, edge_values) :: Tuple{Array{Float64, N}, Array{Float64, Nplus1}},
+                box :: Box{N, Float64}, gridsize :: Dims{N}) where {N, Nplus1}
+        return new{N,Nplus1}(node_values, edge_values, UniformGrid{N,Nplus1}(box, gridsize))
+    end
+
+    function UniformP2{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+        grid = UniformGrid{N,Nplus1}(box, gridsize)
+        na = Array{Float64, N}(undef, gridsize)
+        ne = Array{Float64, Nplus1}(undef, max_edge_id(grid))
+        return new{N,Nplus1}(na, ne, grid)
+    end
+end
+
+const LineP2 = UniformP2{1,2}
+const PlanarP2 = UniformP2{2,3}
+
+# f is not specified a type to allow callable objects and functions
+function UniformP2{N,Nplus1}(f, box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1}
+    p = UniformP2{N,Nplus1}(box, gridsize)
+    discretise!(f, p)
+    return p
+end
+
+# 1D model function
+
+@inline function p2powers(x :: Location{1}) :: NTuple{3, Float64}
+    return (1, x[1], x[1]*x[1])
+end
+
+@inline function p2powers_diff(x :: Location{1}) :: NTuple{3, Float64}
+    return ((0, 1, 2*x[1]),)
+end
+
+@inline function p2_simplicialf(xv :: Vararg{P1NodalData{1},3}) :: NTuple{3, Float64}
+    Ab = (((p2powers(xv[i][1])..., xv[i][2]) for i = 1:3)...,)
+    return linsolve(Ab)
+end
+
+# 2D model function
+
+@inline function p2powers(x :: Location{2}) #:: NTuple{6, Float64}
+    return (1.0, x[1], x[2], x[1]*x[1], x[1]*x[2], x[2]*x[2])
+end
+
+function p2powers_diff(x :: Location{2}) :: NTuple{6, Float64}
+    return ((0.0, 1.0, 0.0, 2*x[1], 2*x[2], 0.0),
+            (0.0, 0.0, 1.0, 0.0, 2*x[1], 2*x[2]))
+end
+
+function p2_simplicialf(xv :: Vararg{P1NodalData{2},6}) :: NTuple{6, Float64}
+    Ab = (((p2powers(xv[i][1])..., xv[i][2]) for i = 1:6)...,)
+    return linsolve(Ab)
+end
+
+# Common
+
+@noinline function p2_modelf(p :: UniformP2{N,Nplus1}, s :: Simplex{N,Nplus1}) where {N,Nplus1}
+    a = (((n.x, p.node_values[n]) for n ∈ nodes(s))...,)
+    b = (((midpoint(e), p.edge_values[e]) for e ∈ edges(s))...,)
+    return  p2_simplicialf(a...,
+                           b...)
+end
+
+@inline function interpolate(p :: UniformP2{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1}
+    a = p2_modelf(p, s)
+    p = p2powers(x)
+    return +((a.*p)...)
+end
+
+function (p :: UniformP2{N,Nplus1})(x :: Location{N}) where {N,Nplus1}
+    return interpolate(p, x)
+end
+
+@inline function differential(p :: UniformP2{N,Nplus1}, s :: Simplex{N,Nplus1},
+                              x :: Location{N}) :: Gradient{N} where {N,Nplus1}
+    a = p2_modelf(p, s)
+    d = p2powers_diff(x)
+    return (+((a.*di)...) for di ∈ d)
+end
+
+@inline function minimise_p2_1d(a, (x₀, x₁))
+    (_, a₁, a₁₁) = a
+    # We do this in cases, first trying for an interior solution, then edges.
+    # For interior solution, first check determinant; no point trying if non-positive
+    if a₁₁ > 0
+        # An interior solution x[1] has to satisfy
+        #   2a₁₁*x[1] + a₁ =0
+        # This gives
+        x = (-a₁/(2*a₁₁),)
+        if in_simplex(t, x)
+            return x, +(a.*p2powers(x))
+        end
+    end
+
+    (x₀, x₁) = nodes(t)
+    v₀ = +((a.*p2powers(x₀))...)
+    v₁ = +((a.*p2powers(x₁))...)
+    return v₀ > v₁ ? (x₀, v₀) : (x₁, v₁)
+end
+
+function minimise(p :: LineP2, t :: RealInterval)
+    minimise_p2_1d(p2_modelf(p, s), nodes(t))
+end
+
+@inline function minimise_p2_edge((a₀, a₁, a₂, a₁₁, a₁₂, a₂₂), x₀, x₁)
+    # Minimise the 2D model on the edge {x₀ + t(x₁ - x₀) | t ∈ [0, 1] }
+    d = x₁ .- x₀
+    b₀ = a₀ + a₁*x₀[1] + a₂*x₀[2] + a₁₁*x₀[1]*x₀[1] + a₁₂*x₀[1]*x₀[2] + a₂₂*x₀[2]*x₀[2]
+    b₁ = a₁*d[1] + a₂*d[2] + 2*a₁₁*d[1]*x₀[1] + a₁₂*(d[1]*x₀[2] + d[2]*x₀[1]) + 2*a₂₂*d[2]*x₀[2]
+    b₁₁ = a₁₁*d[1]*d[1] + a₁₂*d[1]*d[2] + a₂₂*d[2]*d[2]
+    t, v = minimise_p2_1d((b₀, b₁, b₁₁), (0, 1))
+    return @.(x₀ + t*d), v
+end
+
+@inline function minimise(p :: PlanarP2, t :: PlanarSimplex)
+    a = (_, a₁, a₂, a₁₁, a₁₂, a₂₂) =  p2_modelf(p, s)
+    # We do this in cases, first trying for an interior solution, then edges.
+
+    # For interior solution, first check determinant; no point trying if non-positive
+    r = 2*(a₁₁*a₂₂-a₁₂*a₁₂)
+    if r > 0
+        # An interior solution (x[1], x[2]) has to satisfy
+        #   2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0
+        # This gives
+        x = ((a₂₂*a₁-a₁₂*a₂)/r, (a₁₂*a₁-a₁₁*a₂)/r)
+        if in_simplex(t, x)
+            return x, +(a.*p2powers(x))
+        end
+    end
+
+    (x₀, x₁, x₁) = nodes(t)
+
+    x₀₁, v₀₁ = minimise_p2_edge(a, x₀, x₁)
+    x₁₂, v₁₂ = minimise_p2_edge(a, x₁, x₂)
+    x₂₀, v₂₀ = minimise_p2_edge(a, x₂, x₀)
+    if v₀₁ > v₁₂
+        return v₀₁ > v₂₀ ? (x₀₁, v₀₁) : (x₂₀, v₂₀)
+    else
+        return v₁₂ > v₂₀ ? (x₁₂, v₁₂) : (x₂₀, v₂₀)
+    end
+end
+
+function discretise!(f, p :: UniformP2{N,Nplus1}) where {N,Nplus1}
+    for n ∈ nodes(p)
+        p.node_values[n] = f(n.x)
+    end
+    for e ∈ edges(p)
+        p.edge_values[e] = f(midpoint(e))
+    end
+end
+
+end # module
\ No newline at end of file