comparison: src/PlanarFE.jl
src/PlanarFE.jl
- changeset 37
- f8be66557e0f
equal
deleted
inserted
replaced
| 36:6dfa8001eed2 | 37:f8be66557e0f |
|---|---|
| 1 """ | |
| 2 Simple planar (2D) finite element discretisation on a box. | |
| 3 """ | |
| 4 module PlanarFE | |
| 5 | |
| 6 import ..Loops: Box | |
| 7 import ..LinSolve: linsolve | |
| 8 using ..Metaprogramming | |
| 9 | |
| 10 export FEModel, | |
| 11 FEUniform, | |
| 12 UniformP0, | |
| 13 PlanarP0, | |
| 14 LineP0, | |
| 15 UniformP1, | |
| 16 PlanarP1, | |
| 17 LineP1, | |
| 18 UniformP2, | |
| 19 PlanarP2, | |
| 20 LineP2, | |
| 21 simplices, | |
| 22 nodes, | |
| 23 interpolate, | |
| 24 differential, | |
| 25 minimise | |
| 26 | |
| 27 const Location{N} = NTuple{N, Float64} | |
| 28 abstract type FEModel{N} end | |
| 29 abstract type FEUniform{N} <: FEModel{N} end | |
| 30 | |
| 31 const Index = Int | |
| 32 const CaseIndex = Int | |
| 33 | |
| 34 const NodeId{N} = NTuple{N, Index} | |
| 35 | |
| 36 struct Node{N} | |
| 37 id :: NodeId{N} | |
| 38 x :: Location{N} | |
| 39 end | |
| 40 | |
| 41 @inline function Base.getindex(a :: Array{T,N}, n :: Node{N}) where {T,N} | |
| 42 return getindex(a, n.id...) | |
| 43 end | |
| 44 | |
| 45 @inline function Base.setindex!(a :: Array{T,N}, x :: T, n :: Node{N}) where {T,N} | |
| 46 return setindex!(a, x, n.id...) | |
| 47 end | |
| 48 | |
| 49 ############################################## | |
| 50 # Simplices | |
| 51 ############################################## | |
| 52 | |
| 53 # (D-1)-dimensional simplex in ℝ^N. D is the node count. | |
| 54 struct Simplex{N,D,Id} | |
| 55 id :: Id | |
| 56 nodes :: NTuple{D, Node{N}} | |
| 57 @inline Simplex{Id}(id :: Id, nodes :: Vararg{Node{N}, D}) where {N,D,Id} = new{N,D,Id}(id, nodes) | |
| 58 @inline Simplex{N,D,Id}(id :: Id, nodes :: Vararg{Node{N}, D}) where {N,D,Id} = new{N,D,Id}(id, nodes) | |
| 59 end | |
| 60 | |
| 61 # Uniformly indexed simplices | |
| 62 const UCubeId{N} = NTuple{N, Int} | |
| 63 const USimplexId{N} = Tuple{UCubeId{N}, CaseIndex} | |
| 64 const USimplex{N,D} = Simplex{N,D,USimplexId{N}} | |
| 65 const PlanarSimplex = USimplex{2,3} | |
| 66 const RealInterval = USimplex{1,2} | |
| 67 | |
| 68 @inline function Base.getindex(a :: Array{T,Nplus1}, s :: USimplex{N,D}) where {T,Nplus1,D,N} | |
| 69 idx = (s.id[1]..., s.id[2]) :: NTuple{Nplus1,Int} | |
| 70 return getindex(a, idx...) | |
| 71 end | |
| 72 | |
| 73 @inline function Base.setindex!(a :: Array{T,Nplus1}, x :: T, s :: USimplex{N,D}) where {T,Nplus1,D,N} | |
| 74 idx = (s.id[1]..., s.id[2]) :: NTuple{Nplus1,Int} | |
| 75 return setindex!(a, x, idx...) | |
| 76 end | |
| 77 | |
| 78 # We do this by invididual cases to allow typing return values | |
| 79 # (N+1 is not calculable at compile time) | |
| 80 @inline function nodes(s :: Simplex{N,D,Id}) where {N,D,Id} | |
| 81 return s.nodes | |
| 82 end | |
| 83 | |
| 84 @inline @generated function dot(x :: Location{N}, y :: Location{N}) where N | |
| 85 return foldl_exprs( :( + ), :( x[$i]*y[$i] ) for i=1:N) | |
| 86 end | |
| 87 | |
| 88 @inline function orthog((x, y) :: Location{2}) | |
| 89 return (y, -x) | |
| 90 end | |
| 91 | |
| 92 @inline function in_simplex(t :: RealInterval, x :: Location{1}) | |
| 93 (x₀, x₁) = nodes(t) | |
| 94 return ((x .- x₀)*(x₁ .- x₀) ≥ 0 && | |
| 95 (x .- x₁)*(x₀ .- x₁) ≥ 0) | |
| 96 end | |
| 97 | |
| 98 @inline function in_simplex(t :: PlanarSimplex, x :: Location{2}) | |
| 99 (x₀, x₁, x₁) = nodes(t) | |
| 100 return (dot(x .- x₀, orthog(x₁ .- x₀)) ≥ 0 && | |
| 101 dot(x .- x₁, orthog(x₂ .- x₁)) ≥ 0 && | |
| 102 dot(x .- x₂, orthog(x₀ .- x₂)) ≥ 0) | |
| 103 end | |
| 104 | |
| 105 | |
| 106 ############################################## | |
| 107 # Simple planar grid and triangulation | |
| 108 ############################################## | |
| 109 | |
| 110 @inline function cellwidth(box :: Box{N, Float64}, gridsize :: Dims{N}) where N | |
| 111 return (((box[i][2]-box[i][1])/(gridsize[i]-1) for i=1:N)...,) | |
| 112 end | |
| 113 | |
| 114 # Cannot calculate N+1 compile-time | |
| 115 struct UniformGrid{N,Nplus1} | |
| 116 box :: Box{N, Float64} | |
| 117 size :: Dims{N} | |
| 118 cellwidth :: NTuple{N, Float64} | |
| 119 | |
| 120 function UniformGrid{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 121 # TODO: Nested tuple access is inefficient | |
| 122 @assert(N+1 == Nplus1) | |
| 123 @assert(all(box[i][2]>box[i][1] for i=1:N)) | |
| 124 @assert(all(gridsize .≥ 2)) | |
| 125 return new{N,Nplus1}(box, gridsize, cellwidth(box, gridsize)) | |
| 126 end | |
| 127 end | |
| 128 | |
| 129 const LineGrid = UniformGrid{1,2} | |
| 130 const PlanarGrid = UniformGrid{2,3} | |
| 131 | |
| 132 # Number of simplices in a cube of dimension N. | |
| 133 n_simplex_cases( :: LineGrid) = 1 :: CaseIndex | |
| 134 n_simplex_cases( :: PlanarGrid) = 2 :: CaseIndex | |
| 135 max_cube_id(grid :: UniformGrid{N, Nplus1}) where {N, Nplus1} = grid.size .- 1 | |
| 136 max_simplex_id(grid :: UniformGrid{N, Nplus1}) where {N, Nplus1} = | |
| 137 (max_cube_id(grid)..., n_simplex_cases(grid)) | |
| 138 | |
| 139 ∈(x :: Location{N}, b :: Box{N}) where N = all(b[i][1] ≤ x[i] ≤ b[i][2] for i=1:N) | |
| 140 ∉(x :: Location{N}, b :: Box{N}) where N = !(x ∈ b) | |
| 141 ∈(x :: Location{N}, grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} = x ∈ grid.box | |
| 142 ∉(x :: Location{N}, grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} = x ∉ grid.box | |
| 143 | |
| 144 @inline function width(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} | |
| 145 return ((grid.box[i][2]-grid.box[i][1] for i=1:N)...,) | |
| 146 end | |
| 147 | |
| 148 @inline function shift(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} | |
| 149 return ((grid.box[i][1] for i=1:N)...,) | |
| 150 end | |
| 151 | |
| 152 function construct_node(grid :: UniformGrid{N,Nplus1}, idx :: NodeId{N}) :: Node{N} where {N,Nplus1} | |
| 153 loc = grid.cellwidth.*(idx.-1).+shift(grid) | |
| 154 return Node{N}(idx, loc) | |
| 155 end | |
| 156 | |
| 157 @inline function construct_simplex(grid :: LineGrid, id :: USimplexId{1}) :: RealInterval | |
| 158 # TODO: optimise case away | |
| 159 (idx, _) = id | |
| 160 return RealInterval(construct_node(grid, idx), construct_node(grid, idx.+(1,))) | |
| 161 end | |
| 162 | |
| 163 | |
| 164 # In 2D we alternate ⧄ and ⧅. The “base” 90° corner is always listed first. | |
| 165 # It has importance for P0 elements. Otherwise arbitrary we proceed clockwise. | |
| 166 # Julia is real GARBAGE. It does tons of memory alloc trying to access this in | |
| 167 # any form, as Array or Tuple. Therefore things are in code below. | |
| 168 # offset₂ = [[ | |
| 169 # [(1, 0), (0, 0), (1, 1)], # ◪ | |
| 170 # [(0, 1), (1, 1), (0, 0)], # ◩ | |
| 171 # ],[ | |
| 172 # [(0, 0), (1, 0), (0, 1)], # ⬕ | |
| 173 # [(1, 1), (1, 0), (0, 1)], # ⬔ | |
| 174 # ]] :: Array{Array{Array{Tuple{Int,Int},1},1},1} | |
| 175 | |
| 176 | |
| 177 @inline function sqform₂(idx :: UCubeId{2}) | |
| 178 return 1+(idx[1]+idx[2])%2 | |
| 179 end | |
| 180 | |
| 181 # 2D simplices are identified by their “leading” node, i.e., the lower-left corner of the grid square | |
| 182 # that contains them, and the “case” number 1 or 2. | |
| 183 function construct_simplex(grid :: PlanarGrid, simplexid :: USimplexId{2}, form :: Index) :: PlanarSimplex | |
| 184 (idx, case) = simplexid | |
| 185 @assert(all(idx .< grid.size)) | |
| 186 @assert(form==1 || form==2) | |
| 187 @assert(case==1 || case==2) | |
| 188 (off1 :: NodeId{2}, off2 :: NodeId{2}, off3 :: NodeId{2}) = (form==1 ? | |
| 189 (case==1 ? | |
| 190 ((1, 0), (0, 0), (1, 1)) # ◪ | |
| 191 : | |
| 192 ((0, 1), (1, 1), (0, 0)) # ◩ | |
| 193 ) : (case==1 ? | |
| 194 ((0, 0), (1, 0), (0, 1)) # ⬕ | |
| 195 : | |
| 196 ((1, 1), (1, 0), (0, 1)) # ⬔ | |
| 197 )) | |
| 198 | |
| 199 return PlanarSimplex(simplexid, | |
| 200 construct_node(grid, idx .+ off1), | |
| 201 construct_node(grid, idx .+ off2), | |
| 202 construct_node(grid, idx .+ off3)) | |
| 203 end | |
| 204 | |
| 205 function construct_simplex(grid :: PlanarGrid, simplexid :: USimplexId{2}) :: PlanarSimplex | |
| 206 (idx, _) = simplexid | |
| 207 return construct_simplex(grid, simplexid, sqform₂(idx)) | |
| 208 end | |
| 209 | |
| 210 | |
| 211 @inline function get_cubeid(grid :: UniformGrid{N,Nplus1}, x :: Location{N}) :: UCubeId{N} where {N,Nplus1} | |
| 212 (( | |
| 213 # Without the :: Int typeassert this generates lots of memory allocs, although a typeof() | |
| 214 # check here confirms that the type is Int. Julia is just as garbage as everything. | |
| 215 max(min(ceil(Index, (x[i]-grid.box[i][1])/grid.cellwidth[i]), grid.size[i] - 1), 1) :: Int | |
| 216 for i=1:N)...,) | |
| 217 end | |
| 218 | |
| 219 | |
| 220 @inline function get_offset(grid :: UniformGrid{N,Nplus1}, cubeid :: UCubeId{N}, x :: Location{N}) :: Location{N} where {N,Nplus1} | |
| 221 (( | |
| 222 (x[i] - (cubeid[i] - 1) * grid.cellwidth[i]) | |
| 223 for i=1:N)...,) | |
| 224 end | |
| 225 | |
| 226 @inline function simplex_at(grid :: LineGrid, x :: Location{1}) :: RealInterval | |
| 227 if x ∉ grid.box | |
| 228 throw("Coordinate $x out of domain $(grid.box)") | |
| 229 else | |
| 230 return construct_simplex(grid, get_cubeid(grid, x)) | |
| 231 end | |
| 232 end | |
| 233 | |
| 234 @inline function simplex_at(grid :: PlanarGrid, x :: Location{2}) :: PlanarSimplex | |
| 235 if x ∉ grid.box | |
| 236 throw("Coordinate $x out of domain $(grid.box)") | |
| 237 else | |
| 238 # Lower left node is is cube id. Therefore maximum cubeid is grid size - 1 | |
| 239 # in each dimension. | |
| 240 cubeid = get_cubeid(grid, x) | |
| 241 off = get_offset(grid, cubeid, x) | |
| 242 form = sqform₂(cubeid) | |
| 243 if form==0 | |
| 244 # ⧄ → ◪ or ◩ | |
| 245 case = (off[1]>off[2] ? 1 : 2 ) :: CaseIndex | |
| 246 else | |
| 247 # ⧅ → ⬕ or ⬔ | |
| 248 case = (off[1]+off[2]<1 ? 1 : 2) :: CaseIndex | |
| 249 end | |
| 250 return construct_simplex(grid, (cubeid, case), form) | |
| 251 end | |
| 252 end | |
| 253 | |
| 254 ############################################## | |
| 255 # Edges for P2 elements | |
| 256 ############################################## | |
| 257 | |
| 258 const Edge{N,Id} = Simplex{N,2,Id} | |
| 259 const UEdgeId{N} = USimplexId{N} | |
| 260 const UEdge{N} = Edge{N,UEdgeId{N}} | |
| 261 | |
| 262 @inline midpoint(e :: Edge{N,Id}) where {N,Id} = (e.nodes[1].x .+ e.nodes[2].x) ./ 2.0 | |
| 263 | |
| 264 # For our purposes, an edge is something useful for putting an interpolation midpoint on. | |
| 265 # This in 1D there's one “edge”, which is just the original simplex. | |
| 266 n_cube_edges( :: LineGrid) = 1 | |
| 267 n_cube_edges( :: PlanarGrid) = 3 | |
| 268 max_edge_id(grid :: UniformGrid{N, Nplus1}) where {N, Nplus1} = (grid.size..., n_cube_edges(grid)) | |
| 269 | |
| 270 @inline function edge_id(node₁ :: NodeId{1}, node₂ :: NodeId{1}) :: UEdgeId{1} | |
| 271 return (min(node₁[1], node₂[1]),) | |
| 272 end | |
| 273 | |
| 274 # Edges are identified by the lower left node of the “least-numbered” cube that | |
| 275 # contains them, along with 1=horizontal, 2=vertical, 3=diagonal. | |
| 276 # Both of the alternating cases of diagonal edges have the same number. | |
| 277 # Since the least-numbered cube identifies edges, only both diagonal the left and | |
| 278 # bottom border of the cube are counted for that cube, so there are 3 possible | |
| 279 # edges counted towards a cube. | |
| 280 function edge_id(idx₁ :: NodeId{2}, idx₂ :: NodeId{2}) :: UEdgeId{2} | |
| 281 # For performance this does not verify that idx₁ and idx₂ are in the | |
| 282 # same simplex, it only uses their ordering! | |
| 283 if idx₂[1] > idx₁[1] # one of ↗︎, →, ↘︎ | |
| 284 if idx₂[2] > idx₁[2] | |
| 285 return (idx₁, 3) # diagonal ↗︎ in cube(idx₁) | |
| 286 elseif idx₂[2] == idx₁[2] | |
| 287 return (idx₁, 1) # lower side → in cube(idx₁) | |
| 288 else # idx₂[2] < idx₁[2] | |
| 289 return ((idx₁[1], idx₂[2]), 3) # ↘︎ in cube(prev(idx₁,idx₂)) | |
| 290 end | |
| 291 elseif idx₂[1] == idx₁[1] # one of ↑, ↓ | |
| 292 if idx₂[2] > idx₁[2] | |
| 293 return (idx₁, 2) # ↑ in cube(idx₁) | |
| 294 elseif idx₂[2] == idx₁[2] | |
| 295 throw("Edge cannot end in starting idx.") | |
| 296 else # idx₂[2] < idx₁[2] | |
| 297 return (idx₂, 2) # ↓ cube(idx₂) | |
| 298 end | |
| 299 else # idx₂[1] < idx₁[1] # one of ↙︎, ←, ↖︎ | |
| 300 if idx₂[2] > idx₁[2] | |
| 301 return ((idx₁[1], idx₂[2]), 3) # diagonal ↖︎ in cube(prev(idx₁,idx₂)) | |
| 302 elseif idx₂[2] == idx₁[2] | |
| 303 return (idx₂, 1) # lower side ← in cube(idx₂) | |
| 304 else # idx₂[2] < idx₁[2] | |
| 305 return (idx₂, 3) # ↙︎ in cube(idx₂) | |
| 306 end | |
| 307 end | |
| 308 end | |
| 309 | |
| 310 function construct_edge(n₁ :: Node{N}, n₂ :: Node{N}) where N | |
| 311 return UEdge{N}(edge_id(n₁.id, n₂.id), n₁, n₂) | |
| 312 end | |
| 313 | |
| 314 ############################################## | |
| 315 # Iteration over simplices and nodes | |
| 316 ############################################## | |
| 317 | |
| 318 struct Iter{G,W} | |
| 319 grid :: G | |
| 320 end | |
| 321 | |
| 322 # @inline function nodes(s :: PlanarSimplex{N}) where N | |
| 323 # return Iter{PlanarSimplex{N}, :nodes}(s) | |
| 324 # end | |
| 325 | |
| 326 # @inline function Base.iterate(sn :: Iter{PlanarSimplex{N}, :nodes}) where N | |
| 327 # return (sn.p.base, Index(0)) | |
| 328 # end | |
| 329 | |
| 330 # @inline function Base.iterate(sn :: Iter{PlanarSimplex{N}, :nodes}, idx :: Index) where N | |
| 331 # if idx==length(sn.p.tips) | |
| 332 # return nothing | |
| 333 # else | |
| 334 # idx=idx+1 | |
| 335 # return (sn.p.tips[idx], idx) | |
| 336 # end | |
| 337 # end | |
| 338 | |
| 339 @inline function nodes(p :: M) where {N, M <: FEUniform{N}} | |
| 340 return nodes(p.grid) | |
| 341 end | |
| 342 | |
| 343 @inline function nodes(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} | |
| 344 return Iter{UniformGrid{N,Nplus1}, Node{N}}(grid) | |
| 345 end | |
| 346 | |
| 347 @inline function first_nodeid( :: UniformGrid{N,Nplus1}) :: NodeId{N} where {N,Nplus1} | |
| 348 return (( 1 :: Int for i=1:N)...,) | |
| 349 end | |
| 350 | |
| 351 # Need to hard-code this for each N to avoid Julia munching memory. | |
| 352 # Neither @generated nor plain code works. | |
| 353 @inline function next_nodeid(grid :: LineGrid, idx :: NodeId{1}) :: Union{Nothing, NodeId{1}} | |
| 354 if idx[1] < grid.size[1] | |
| 355 return (idx[1]+1,) | |
| 356 else | |
| 357 return nothing | |
| 358 end | |
| 359 end | |
| 360 | |
| 361 @inline function next_nodeid(grid :: PlanarGrid, idx :: NodeId{2}) :: Union{Nothing, NodeId{2}} | |
| 362 if idx[1] < grid.size[1] | |
| 363 return (idx[1]+1, idx[2]) | |
| 364 elseif idx[2] < grid.size[2] | |
| 365 return (1 :: Int, idx[2]+1) | |
| 366 else | |
| 367 return nothing | |
| 368 end | |
| 369 end | |
| 370 | |
| 371 function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Node{N}}) :: | |
| 372 Union{Nothing, Tuple{Node{N}, NodeId{N}}} where {N,Nplus1} | |
| 373 idx = first_nodeid(iter.grid) | |
| 374 return construct_node(iter.grid, idx), idx | |
| 375 end | |
| 376 | |
| 377 function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Node{N}}, idx :: NodeId{N}) :: | |
| 378 Union{Nothing, Tuple{Node{N}, NodeId{N}}} where {N,Nplus1} | |
| 379 nidx = next_nodeid(iter.grid, idx) | |
| 380 if isnothing(nidx) | |
| 381 return nothing | |
| 382 else | |
| 383 return construct_node(iter.grid, nidx), nidx | |
| 384 end | |
| 385 end | |
| 386 | |
| 387 # @generated function Base.iterate(iter :: Iter{N, Node{N}}, idx :: NodeId{N}) where N | |
| 388 # # This generation reduces memory allocations | |
| 389 # do_coord(i) = :( | |
| 390 # if idx[$i] < iter.grid.size[$i] | |
| 391 # nidx = $(lift_exprs(( (:( 1 :: Int ) for _=1:i-1)..., :( idx[$i]+1 ), ( :( idx[$j] ) for j=i+1:N)...,))) :: NTuple{N, Int} | |
| 392 # return construct_node(iter.grid, nidx), nidx | |
| 393 # end | |
| 394 # ) | |
| 395 | |
| 396 # return quote | |
| 397 # $(sequence_exprs( do_coord(i) for i=1:N )) | |
| 398 # return nothing | |
| 399 # end | |
| 400 # end | |
| 401 | |
| 402 # function Base.iterate(iter :: Iter{N, Node{N}}, idx :: NodeId{N}) where N | |
| 403 # for i=1:N | |
| 404 # if idx[i]<iter.grid.size[i] | |
| 405 # nidx = (( 1 :: Int for _=1:i-1)..., idx[i]+1, idx[i+1:N]...) | |
| 406 # return construct_node(iter.grid, nidx), nidx | |
| 407 # end | |
| 408 # end | |
| 409 | |
| 410 # return nothing | |
| 411 # end | |
| 412 | |
| 413 @inline function simplices(p :: M) where {N, M <: FEUniform{N}} | |
| 414 return simplices(p.grid) | |
| 415 end | |
| 416 | |
| 417 @inline function simplices(grid :: UniformGrid{N,Nplus1}) where {N,Nplus1} | |
| 418 return Iter{UniformGrid{N,Nplus1}, USimplex{N,D}}(grid) | |
| 419 end | |
| 420 | |
| 421 @inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Simplex{N,Nplus1}}) :: | |
| 422 Union{Nothing, Tuple{USimplex{N,Nplus1}, USimplexId{N}}} where {N,Nplus1} | |
| 423 cubeid = ((Index(1) for i=1:N)...,) | |
| 424 simplexid = (cubeid, Index(1)) | |
| 425 return (construct_simplex(iter.grid, simplexid), simplexid) | |
| 426 end | |
| 427 | |
| 428 @inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, Simplex{N,Nplus1}}, | |
| 429 (cubeid, case) :: USimplexId{N}) :: | |
| 430 Union{Nothing, Tuple{USimplex{N,Nplus1}, USimplexId{N}}} where {N,Nplus1} | |
| 431 if case < n_simplex_cases(iter.grid) | |
| 432 simplexid = (cubeid, case+1) | |
| 433 return (construct_simplex(iter.grid, simplexid), simplexid) | |
| 434 else | |
| 435 for i=1:N | |
| 436 if cubeid[i]+1<iter.grid.size[i] | |
| 437 new_cubeid=((Index(1) for _=1:i-1)..., cubeid[i]+1, cubeid[i+1:N]...) | |
| 438 simplexid = (new_cubeid, Index(1)) | |
| 439 return (construct_simplex(iter.grid, simplexid), simplexid) | |
| 440 end | |
| 441 end | |
| 442 end | |
| 443 | |
| 444 return nothing | |
| 445 end | |
| 446 | |
| 447 # TODO: wrong for far-side boundary? Need to do by cumbe numbers | |
| 448 @inline function edges(s :: RealInterval) :: Tuple{UEdge{1}} | |
| 449 return (construct_edge(s.nodes...),) | |
| 450 end | |
| 451 | |
| 452 # TODO: Edge -> Simplex | |
| 453 @inline function edges(s :: PlanarSimplex) :: NTuple{3,UEdge{2}} | |
| 454 return (construct_edge(s.nodes[1], s.nodes[2]), | |
| 455 construct_edge(s.nodes[1], s.nodes[3]), | |
| 456 construct_edge(s.nodes[2], s.nodes[3])) | |
| 457 end | |
| 458 | |
| 459 @inline function edges(p :: M) where {N, M <: FEUniform{N}} | |
| 460 return edges(p.grid) | |
| 461 end | |
| 462 | |
| 463 @inline function edges(grid :: LineGrid) | |
| 464 return Iter{LineGrid, UEdge{1}}(grid) | |
| 465 end | |
| 466 | |
| 467 @inline function edges(grid :: PlanarGrid) | |
| 468 return Iter{PlanarGrid, UEdge{2}}(grid) | |
| 469 end | |
| 470 | |
| 471 @inline function UEdge(grid :: LineGrid, edgeid :: UEdgeId{1}) :: UEdge{1} | |
| 472 (cubeid, case) = edgeid | |
| 473 @assert(case==1) | |
| 474 return UEdge(edgeid, construct_node(grid, cubeid), construct_node(grid, cubeid.+(1,))) | |
| 475 end | |
| 476 | |
| 477 @inline function next_edge_case(grid :: LineGrid, (cubeid, case) :: UEdgeId{1}) :: CaseIndex | |
| 478 return -1 | |
| 479 end | |
| 480 | |
| 481 @inline function first_edge_case(grid :: LineGrid, cubeid :: UCubeId{1}) :: CaseIndex | |
| 482 return cubeid[1] < grid.size[1] ? 1 : -1 | |
| 483 end | |
| 484 | |
| 485 @inline function UEdge(grid :: PlanarGrid, edgeid :: UEdgeId{2}) :: UEdge{2} | |
| 486 (cubeid, case) = edgeid | |
| 487 e(idx1, idx2) = UEdge{2}(edgeid, construct_node(grid, idx1), construct_node(grid, idx2)) | |
| 488 if case==1 | |
| 489 return e(cubeid, cubeid.+(1,0)) | |
| 490 elseif case==2 | |
| 491 return e(cubeid, cubeid.+(0,1)) | |
| 492 elseif case==3 | |
| 493 if sqform₂(cubeid)==1 | |
| 494 return e(cubeid, cubeid.+(1,1)) | |
| 495 else | |
| 496 return e(cubeid.+(1,0), cubeid.+(0,1)) | |
| 497 end | |
| 498 else | |
| 499 throw("Invalid edge") | |
| 500 end | |
| 501 end | |
| 502 | |
| 503 @inline function first_edge_case(grid :: PlanarGrid, cubeid :: UCubeId{2}) :: CaseIndex | |
| 504 if cubeid[1]<grid.size[1] | |
| 505 return 1 | |
| 506 elseif cubeid[2]<grid.size[2] | |
| 507 return 2 | |
| 508 else | |
| 509 return -1 | |
| 510 end | |
| 511 end | |
| 512 | |
| 513 | |
| 514 @inline function next_edge_case(grid :: PlanarGrid, (cubeid, case) :: UEdgeId{2}) :: CaseIndex | |
| 515 if case<3 && cubeid[1]<grid.size[1] && cubeid[2]<grid.size[2] | |
| 516 return case + 1 | |
| 517 else | |
| 518 return -1 | |
| 519 end | |
| 520 end | |
| 521 | |
| 522 | |
| 523 @inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, UEdge{N}}) :: | |
| 524 Union{Nothing, Tuple{UEdge{N}, UEdgeId{N}}} where {N,Nplus1} | |
| 525 cubeid = first_nodeid(iter.grid) | |
| 526 case = first_edge_case(iter.grid, cubeid) | |
| 527 @assert(case ≥ 1) | |
| 528 edgeid = (cubeid, case) | |
| 529 return (UEdge(iter.grid, edgeid), edgeid) | |
| 530 end | |
| 531 | |
| 532 @inline function Base.iterate(iter :: Iter{UniformGrid{N,Nplus1}, UEdge{N}}, (cubeid, case) :: UEdgeId{N}) :: | |
| 533 Union{Nothing, Tuple{UEdge{N}, UEdgeId{N}}} where {N,Nplus1} | |
| 534 case = next_edge_case(iter.grid, (cubeid, case)) | |
| 535 if case ≥ 1 | |
| 536 edgeid = (cubeid, case) | |
| 537 return (UEdge(iter.grid, edgeid), edgeid) | |
| 538 else | |
| 539 while true | |
| 540 cubeid = next_nodeid(iter.grid, cubeid) | |
| 541 if isnothing(cubeid) | |
| 542 return nothing | |
| 543 end | |
| 544 case = first_edge_case(iter.grid, cubeid) :: Index | |
| 545 if case ≥ 1 | |
| 546 edgeid = (cubeid, case) :: UEdgeId{N} | |
| 547 return (UEdge(iter.grid, edgeid), edgeid) | |
| 548 end | |
| 549 end | |
| 550 end | |
| 551 | |
| 552 return nothing | |
| 553 end | |
| 554 | |
| 555 | |
| 556 ############################################## | |
| 557 # Evaluation, differentiation, minimisation | |
| 558 ############################################## | |
| 559 | |
| 560 function interpolate(p :: M, x :: Location{N}) where {N, M <: FEModel{N}} | |
| 561 return interpolate(p, simplex_at(p.grid, x), x) | |
| 562 end | |
| 563 | |
| 564 # Apparent Julia bug / incomplete implementation: | |
| 565 # https://stackoverflow.com/questions/46063872/parametric-functors-in-julia | |
| 566 # How to do individually. | |
| 567 # function (p :: M)(x :: Location{N}) where {N, M <: FEModel{N}} | |
| 568 # return interpolate(p, x) | |
| 569 # end | |
| 570 | |
| 571 function differential(p :: M, x :: Location{N}) where {N, M <: FEModel{N}} | |
| 572 differential(p, simplex_at(p.grid, x), x) | |
| 573 end | |
| 574 | |
| 575 function minimise(p :: M) where {N, M <: FEModel{N}} | |
| 576 minloc = nothing | |
| 577 minv = nothing | |
| 578 for simplex ∈ simplices(p) | |
| 579 thisminloc, thisminv = minimise(p, simplex) | |
| 580 if isnothing(minv) || minv < thisminv | |
| 581 minv=thisminv | |
| 582 minloc=thisminloc | |
| 583 end | |
| 584 end | |
| 585 return minv, minloc | |
| 586 end | |
| 587 | |
| 588 ############################################## | |
| 589 # P0 elements | |
| 590 ############################################## | |
| 591 | |
| 592 struct UniformP0{N,Nplus1} <: FEUniform{N} | |
| 593 simplex_values :: Array{Float64, Nplus1} | |
| 594 grid :: UniformGrid{N,Nplus1} | |
| 595 | |
| 596 function UniformP0{N,Nplus1}(simplex_values :: Array{Float64, Nplus1}, | |
| 597 box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 598 grid = UniformGrid{N,Nplus1}(box, gridsize) | |
| 599 return new{N,Nplus1}(simplex_values, grid) | |
| 600 end | |
| 601 | |
| 602 function UniformP0{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 603 grid = UniformGrid{N,Nplus1}(box, gridsize) | |
| 604 a = Array{Float64, Nplus1}(undef, max_simplex_id(grid)) | |
| 605 return new{N,Nplus1}(simplex_values, grid) | |
| 606 end | |
| 607 | |
| 608 end | |
| 609 | |
| 610 const LineP0 = UniformP0{1,1} | |
| 611 const PlanarP0 = UniformP0{2,3} | |
| 612 | |
| 613 # f is not specified a type to allow callable objects and functions | |
| 614 function UniformP0{N,Nplus1}(f, box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 615 p = UniformP0{N,Nplus1}(box, gridsize) | |
| 616 discretise!(f, p) | |
| 617 return p | |
| 618 end | |
| 619 | |
| 620 function (p :: UniformP0{N,Nplus1})(x :: Location{N}) where {N,Nplus1} | |
| 621 return interpolate(p, x) | |
| 622 end | |
| 623 | |
| 624 @inline function interpolate(p :: UniformP0{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1} | |
| 625 return p.simplex_values[s] | |
| 626 end | |
| 627 | |
| 628 function differential(p :: UniformP0{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1} | |
| 629 return ((0.0 for i=1:N)...,) | |
| 630 end | |
| 631 | |
| 632 @inline function minimise(p :: UniformP0{N,Nplus1}, s :: Simplex{N,Nplus1}) where {N,Nplus1} | |
| 633 return s.nodes[1].x, p.simplex_values[s] | |
| 634 end | |
| 635 | |
| 636 function discretise!(f, p :: UniformP0{N,Nplus1}) where {N,Nplus1} | |
| 637 for s ∈ simplices(p) | |
| 638 p.simplex_values[s]=f(s.nodes[1].x) | |
| 639 end | |
| 640 end | |
| 641 | |
| 642 ############################################## | |
| 643 # P1 elements | |
| 644 ############################################## | |
| 645 | |
| 646 struct UniformP1{N,Nplus1} <: FEUniform{N} | |
| 647 node_values :: Array{Float64, N} | |
| 648 grid :: UniformGrid{N,Nplus1} | |
| 649 | |
| 650 function UniformP1{N,Nplus1}(node_values :: Array{Float64, N}, | |
| 651 box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 652 grid = UniformGrid{N,Nplus1}(box, gridsize) | |
| 653 return new{N,Nplus1}(node_values, grid) | |
| 654 end | |
| 655 | |
| 656 function UniformP1{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 657 return UniformP1{N,Nplus1}(Array{Float64, N}(undef, gridsize), box, gridsize) | |
| 658 end | |
| 659 end | |
| 660 | |
| 661 const LineP1 = UniformP1{1,1} | |
| 662 const PlanarP1 = UniformP1{2,3} | |
| 663 | |
| 664 # f is not specified a type to allow callable objects and functions | |
| 665 function UniformP1{N,Nplus1}(f, box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 666 p = UniformP1{N,Nplus1}(box, gridsize) | |
| 667 discretise!(f, p) | |
| 668 return p | |
| 669 end | |
| 670 | |
| 671 const P1NodalData{N} = Tuple{Location{N}, Float64} | |
| 672 | |
| 673 # 1D model functions | |
| 674 | |
| 675 @inline function p1_simplicialf(((x₀,), v₀) :: P1NodalData{1}, | |
| 676 ((x₁,), v₁) :: P1NodalData{1}) :: NTuple{2,Float64} | |
| 677 a₁ = (v₁-v₀)/(x₁-x₀) | |
| 678 a₀ = v₀ - a₁*x₀ | |
| 679 return a₀, a₁ | |
| 680 end | |
| 681 | |
| 682 function p1_modelf(p :: UniformP1{1}, s :: RealInterval) | |
| 683 n₀, n₁ = nodes(s) | |
| 684 return p1_simplicialf((n₀.x, p.node_values[n₀]), | |
| 685 (n₁.x, p.node_values[n₁])) | |
| 686 end | |
| 687 | |
| 688 | |
| 689 # 2D model functions | |
| 690 | |
| 691 @inline function p1_simplicialf((x₀, v₀) :: P1NodalData{2}, | |
| 692 (x₁, v₁) :: P1NodalData{2}, | |
| 693 (x₂, v₂) :: P1NodalData{2}) :: NTuple{3,Float64} | |
| 694 d = x₁ .- x₀ | |
| 695 q = x₂ .- x₀ | |
| 696 r = q[1]*d[2] - d[1]*q[2] | |
| 697 a₁ = ((v₂-v₀)*d[2] - (v₁-v₀)*q[2])/r | |
| 698 a₂ = ((v₁-v₀)*q[1] - (v₂-v₀)*d[1])/r | |
| 699 a₀ = v₀-a₁*x₀[1]-a₂*x₀[2] | |
| 700 return a₀, a₁, a₂ | |
| 701 end | |
| 702 | |
| 703 @inline function p1_modelf(p :: UniformP1{N,Nplus1}, s :: Simplex{N,Nplus1}) where {N,Nplus1} | |
| 704 return p1_simplicialf(((n.x, p.node_values[n]) for n ∈ nodes(s))...) | |
| 705 end | |
| 706 | |
| 707 # @generated function p1_modelf(p :: UniformP1{N}, s :: PlanarSimplex{N}) where N | |
| 708 # return quote | |
| 709 # ns = nodes(s) | |
| 710 # p1_simplicialf($(lift_exprs( :((ns[$i].x, p.node_values[ns[$i].id...])) for i=1:3 ))...) | |
| 711 # end | |
| 712 # end | |
| 713 | |
| 714 #@inline function g(a, (i,)) | |
| 715 # return a[i] | |
| 716 #end | |
| 717 | |
| 718 # @inline function g(a, (i, j)) | |
| 719 # return a[i,j] | |
| 720 # end | |
| 721 | |
| 722 # function p1_modelf(p :: UniformP1{2}, s :: PlanarSimplex{2}) | |
| 723 # n₀, n₁, n₂ = nodes(s) | |
| 724 # a, b = n₀.id[1], n₀.id[2] | |
| 725 # v₀ = p.node_values[a, b] | |
| 726 # a, b = n₁.id[1], n₁.id[2] | |
| 727 # v₁ = p.node_values[a, b] | |
| 728 # a, b = n₂.id[1], n₂.id[2] | |
| 729 # v₂ = p.node_values[a, b] | |
| 730 # return p1_simplicialf(n₀.x, v₀, | |
| 731 # n₁.x, v₁, | |
| 732 # n₂.x, v₂) | |
| 733 # end | |
| 734 | |
| 735 | |
| 736 # Common | |
| 737 | |
| 738 function (p :: UniformP1{N,Nplus1})(x :: Location{N}) where {N,Nplus1} | |
| 739 return interpolate(p, x) | |
| 740 end | |
| 741 | |
| 742 function interpolate(p :: UniformP1{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1} | |
| 743 a = p1_modelf(p, s) | |
| 744 return +((a.*(1, x...))...) | |
| 745 end | |
| 746 | |
| 747 @inline function differential(p :: UniformP1{N,Nplus1}, s :: Simplex{N,Nplus1}, :: Location{N}) :: Gradient{N} where {N,Nplus1} | |
| 748 a = p1_modelf(p, s) | |
| 749 return a[2:end] | |
| 750 end | |
| 751 | |
| 752 @inline function minimise(p :: UniformP1{N,Nplus1}, t :: Simplex{N,Nplus1}) where {N,Nplus1} | |
| 753 bestnode = s.nodes[1] | |
| 754 bestvalue = p.node_values[bestnode] | |
| 755 # The minima is found at a node, so only need to go through them | |
| 756 for node ∈ s.nodes[2:end] | |
| 757 value = p.node_values[node] | |
| 758 if value < bestvalue | |
| 759 bestvalue = value | |
| 760 bestnode = node | |
| 761 end | |
| 762 end | |
| 763 return bestnode.x, bestvalue | |
| 764 end | |
| 765 | |
| 766 function discretise!(f, p :: UniformP1{N,Nplus1}) where {N,Nplus1} | |
| 767 for n ∈ nodes(p) | |
| 768 p.node_values[n]=f(n.x) | |
| 769 end | |
| 770 end | |
| 771 | |
| 772 ############################################## | |
| 773 # P2 elements | |
| 774 ############################################## | |
| 775 | |
| 776 const Gradient{N} = NTuple{N, Float64} | |
| 777 const P2NodalData{N} = Tuple{Location{N}, Float64, Gradient{N}} | |
| 778 const LocationAndGradient{N} = Tuple{Location{N}, Gradient{N}} | |
| 779 | |
| 780 struct UniformP2{N,Nplus1} <: FEUniform{N} | |
| 781 node_values :: Array{Float64, N} | |
| 782 edge_values :: Array{Float64, Nplus1} | |
| 783 grid :: UniformGrid{N,Nplus1} | |
| 784 | |
| 785 function UniformP2{N,Nplus1}( | |
| 786 (node_values, edge_values) :: Tuple{Array{Float64, N}, Array{Float64, Nplus1}}, | |
| 787 box :: Box{N, Float64}, gridsize :: Dims{N}) where {N, Nplus1} | |
| 788 return new{N,Nplus1}(node_values, edge_values, UniformGrid{N,Nplus1}(box, gridsize)) | |
| 789 end | |
| 790 | |
| 791 function UniformP2{N,Nplus1}(box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 792 grid = UniformGrid{N,Nplus1}(box, gridsize) | |
| 793 na = Array{Float64, N}(undef, gridsize) | |
| 794 ne = Array{Float64, Nplus1}(undef, max_edge_id(grid)) | |
| 795 return new{N,Nplus1}(na, ne, grid) | |
| 796 end | |
| 797 end | |
| 798 | |
| 799 const LineP2 = UniformP2{1,2} | |
| 800 const PlanarP2 = UniformP2{2,3} | |
| 801 | |
| 802 # f is not specified a type to allow callable objects and functions | |
| 803 function UniformP2{N,Nplus1}(f, box :: Box{N, Float64}, gridsize :: Dims{N}) where {N,Nplus1} | |
| 804 p = UniformP2{N,Nplus1}(box, gridsize) | |
| 805 discretise!(f, p) | |
| 806 return p | |
| 807 end | |
| 808 | |
| 809 # 1D model function | |
| 810 | |
| 811 @inline function p2powers(x :: Location{1}) :: NTuple{3, Float64} | |
| 812 return (1, x[1], x[1]*x[1]) | |
| 813 end | |
| 814 | |
| 815 @inline function p2powers_diff(x :: Location{1}) :: NTuple{3, Float64} | |
| 816 return ((0, 1, 2*x[1]),) | |
| 817 end | |
| 818 | |
| 819 @inline function p2_simplicialf(xv :: Vararg{P1NodalData{1},3}) :: NTuple{3, Float64} | |
| 820 Ab = (((p2powers(xv[i][1])..., xv[i][2]) for i = 1:3)...,) | |
| 821 return linsolve(Ab) | |
| 822 end | |
| 823 | |
| 824 # 2D model function | |
| 825 | |
| 826 @inline function p2powers(x :: Location{2}) #:: NTuple{6, Float64} | |
| 827 return (1.0, x[1], x[2], x[1]*x[1], x[1]*x[2], x[2]*x[2]) | |
| 828 end | |
| 829 | |
| 830 function p2powers_diff(x :: Location{2}) :: NTuple{6, Float64} | |
| 831 return ((0.0, 1.0, 0.0, 2*x[1], 2*x[2], 0.0), | |
| 832 (0.0, 0.0, 1.0, 0.0, 2*x[1], 2*x[2])) | |
| 833 end | |
| 834 | |
| 835 function p2_simplicialf(xv :: Vararg{P1NodalData{2},6}) :: NTuple{6, Float64} | |
| 836 Ab = (((p2powers(xv[i][1])..., xv[i][2]) for i = 1:6)...,) | |
| 837 return linsolve(Ab) | |
| 838 end | |
| 839 | |
| 840 # Common | |
| 841 | |
| 842 @noinline function p2_modelf(p :: UniformP2{N,Nplus1}, s :: Simplex{N,Nplus1}) where {N,Nplus1} | |
| 843 a = (((n.x, p.node_values[n]) for n ∈ nodes(s))...,) | |
| 844 b = (((midpoint(e), p.edge_values[e]) for e ∈ edges(s))...,) | |
| 845 return p2_simplicialf(a..., | |
| 846 b...) | |
| 847 end | |
| 848 | |
| 849 @inline function interpolate(p :: UniformP2{N,Nplus1}, s :: Simplex{N,Nplus1}, x :: Location{N}) where {N,Nplus1} | |
| 850 a = p2_modelf(p, s) | |
| 851 p = p2powers(x) | |
| 852 return +((a.*p)...) | |
| 853 end | |
| 854 | |
| 855 function (p :: UniformP2{N,Nplus1})(x :: Location{N}) where {N,Nplus1} | |
| 856 return interpolate(p, x) | |
| 857 end | |
| 858 | |
| 859 @inline function differential(p :: UniformP2{N,Nplus1}, s :: Simplex{N,Nplus1}, | |
| 860 x :: Location{N}) :: Gradient{N} where {N,Nplus1} | |
| 861 a = p2_modelf(p, s) | |
| 862 d = p2powers_diff(x) | |
| 863 return (+((a.*di)...) for di ∈ d) | |
| 864 end | |
| 865 | |
| 866 @inline function minimise_p2_1d(a, (x₀, x₁)) | |
| 867 (_, a₁, a₁₁) = a | |
| 868 # We do this in cases, first trying for an interior solution, then edges. | |
| 869 # For interior solution, first check determinant; no point trying if non-positive | |
| 870 if a₁₁ > 0 | |
| 871 # An interior solution x[1] has to satisfy | |
| 872 # 2a₁₁*x[1] + a₁ =0 | |
| 873 # This gives | |
| 874 x = (-a₁/(2*a₁₁),) | |
| 875 if in_simplex(t, x) | |
| 876 return x, +(a.*p2powers(x)) | |
| 877 end | |
| 878 end | |
| 879 | |
| 880 (x₀, x₁) = nodes(t) | |
| 881 v₀ = +((a.*p2powers(x₀))...) | |
| 882 v₁ = +((a.*p2powers(x₁))...) | |
| 883 return v₀ > v₁ ? (x₀, v₀) : (x₁, v₁) | |
| 884 end | |
| 885 | |
| 886 function minimise(p :: LineP2, t :: RealInterval) | |
| 887 minimise_p2_1d(p2_modelf(p, s), nodes(t)) | |
| 888 end | |
| 889 | |
| 890 @inline function minimise_p2_edge((a₀, a₁, a₂, a₁₁, a₁₂, a₂₂), x₀, x₁) | |
| 891 # Minimise the 2D model on the edge {x₀ + t(x₁ - x₀) | t ∈ [0, 1] } | |
| 892 d = x₁ .- x₀ | |
| 893 b₀ = a₀ + a₁*x₀[1] + a₂*x₀[2] + a₁₁*x₀[1]*x₀[1] + a₁₂*x₀[1]*x₀[2] + a₂₂*x₀[2]*x₀[2] | |
| 894 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] | |
| 895 b₁₁ = a₁₁*d[1]*d[1] + a₁₂*d[1]*d[2] + a₂₂*d[2]*d[2] | |
| 896 t, v = minimise_p2_1d((b₀, b₁, b₁₁), (0, 1)) | |
| 897 return @.(x₀ + t*d), v | |
| 898 end | |
| 899 | |
| 900 @inline function minimise(p :: PlanarP2, t :: PlanarSimplex) | |
| 901 a = (_, a₁, a₂, a₁₁, a₁₂, a₂₂) = p2_modelf(p, s) | |
| 902 # We do this in cases, first trying for an interior solution, then edges. | |
| 903 | |
| 904 # For interior solution, first check determinant; no point trying if non-positive | |
| 905 r = 2*(a₁₁*a₂₂-a₁₂*a₁₂) | |
| 906 if r > 0 | |
| 907 # An interior solution (x[1], x[2]) has to satisfy | |
| 908 # 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 | |
| 909 # This gives | |
| 910 x = ((a₂₂*a₁-a₁₂*a₂)/r, (a₁₂*a₁-a₁₁*a₂)/r) | |
| 911 if in_simplex(t, x) | |
| 912 return x, +(a.*p2powers(x)) | |
| 913 end | |
| 914 end | |
| 915 | |
| 916 (x₀, x₁, x₁) = nodes(t) | |
| 917 | |
| 918 x₀₁, v₀₁ = minimise_p2_edge(a, x₀, x₁) | |
| 919 x₁₂, v₁₂ = minimise_p2_edge(a, x₁, x₂) | |
| 920 x₂₀, v₂₀ = minimise_p2_edge(a, x₂, x₀) | |
| 921 if v₀₁ > v₁₂ | |
| 922 return v₀₁ > v₂₀ ? (x₀₁, v₀₁) : (x₂₀, v₂₀) | |
| 923 else | |
| 924 return v₁₂ > v₂₀ ? (x₁₂, v₁₂) : (x₂₀, v₂₀) | |
| 925 end | |
| 926 end | |
| 927 | |
| 928 function discretise!(f, p :: UniformP2{N,Nplus1}) where {N,Nplus1} | |
| 929 for n ∈ nodes(p) | |
| 930 p.node_values[n] = f(n.x) | |
| 931 end | |
| 932 for e ∈ edges(p) | |
| 933 p.edge_values[e] = f(midpoint(e)) | |
| 934 end | |
| 935 end | |
| 936 | |
| 937 end # module |
