--- /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