Developer Documentation

IsubdBase(corners, mapping[, adjacency])

The base domain of an implicit LEB subdivision: one reference-space corner triple per base triangle, ordered for bisection (see leb_order), plus the geometry mapping(base_id, ξ) -> physical point (for FEM cells the geometric map, optionally composed with a warp field). The mapping is applied per decoded vertex, so curved cells subdivide into curved sub-triangles.

adjacency[b][e] names the base triangle across edge e of base triangle b (EDGE_S/EDGE_L/EDGE_R) as (base, edge, reversed), or NO_NEIGHBOR for a boundary. Supplying it enables conforming (watertight) refinement — see refine_keys!. It must be compatible: a split edge may only pair with another split edge, and a leg only with legs. _build_substrate builds such a table for FEData; pairings that violate compatibility are dropped to boundaries, which costs conformity along those edges but nothing else.

IsubdMesh(base::IsubdBase)

The reusable output buffers of decode_keys!: vertex positions (in the mapping's image space), per-vertex reference coordinates and owning base triangle, one face per key, and the lookup table backing vertex sharing.

update_keys!(out, keys, base, lod; max_depth=LEB_MAX_DEPTH, hysteresis=0.0)

One split/merge/keep streaming pass over the sorted key buffer (the GPU compute pass of the demo): every key either emits its two children (its excess_levels is positive), its parent (the parent's excess is at most -hysteresis; only child 0 emits it, so the pair collapses to one key), or itself. Refinement moves at most one level per pass — drive it with refine_keys! to reach the steady state.

Split (excess(key) > 0) and merge (excess(parent) ≤ -hysteresis) test disjoint predicates even at hysteresis = 0 — a parent that just split has positive excess, so its children never immediately merge back; the steady state is a true fixed point. A positive hysteresis additionally keeps keys whose parent hovers around excess 0 from toggling under camera jitter, at the price of the merged state lagging the split state by up to that many levels.

A pair merges only when both siblings are present as leaves, which the sorted order makes an adjacent-element check. The demo omits this and relies on its LoD never jumping levels between siblings; under a criterion with sharp spatial variation (an error estimator on rough data), the unguarded merge lets a parent overlap its sibling's still-deeper subtree, or drops a child while the sibling subtree persists — converging to a state that double-covers or holes the domain. (A GPU port checks sibling presence on the concurrent-binary-tree bitfield instead of the sorted buffer.)

refine_keys!(keys, scratch, base, lod; max_depth, hysteresis, max_passes, conforming) -> keys

Run passes until the key buffer is a fixed point of the criterion (or max_passes is hit — one more than max_depth suffices for any monotone criterion). keys is updated in place, scratch is the ping-pong buffer.

With conforming = true (the default whenever base carries an adjacency table) the passes are conforming_update! and the result is watertight; otherwise they are the streaming update_keys!, whose refinement-level boundaries leave gaps bounded by the criterion's tolerance.

decode_keys!(mesh::IsubdMesh, keys, base; groups=nothing) -> mesh

Emit the key buffer as a drawable mesh: positions, their reference coordinates, and one face per key, with each bisection's handedness flip undone so the winding stays consistent.

Vertices are shared within a groups class and duplicated across classes. FEData passes the owning cell, which shares everything inside a cell — where the drawn field is continuous — while keeping element boundaries duplicated, so discontinuous (L2/DG) fields keep their jumps exactly as in the static tessellation. A conforming mesh has about half as many vertices as triangles, against three per triangle unshared, and every vertex costs a geometry evaluation here and a field evaluation downstream. With groups=nothing each base triangle forms its own class, sharing within its own subtree and duplicating along base edges.

Buffers (including the lookup table) are emptied rather than reallocated, so re-decoding a steady mesh does not grow the heap.

decode_topology!(mesh, keys, base; groups=nothing) -> mesh

The connectivity half of decode_keys!: reference coordinates, their owning base triangle, and the faces. This is the part that depends only on the key set, so a consumer whose mesh is unchanged can re-run just decode_positions! and skip the vertex-sharing lookups entirely.

decode_positions!(mesh, base) -> mesh

Map the vertices laid out by decode_topology! through base.mapping. Cheap to repeat: it touches one point per vertex and no lookup table, which is what an unchanged mesh under a changing solution needs.

UniformLoD(target)

Refine everything to exactly target bisection levels. Deterministic; useful for testing and as a static-refinement fallback.

DeviationLoD(f, tol[, cache]; samples = DEVIATION_SAMPLES)

Split until the triangle's linear interpolation approximates f(base_id, ξ) to within tol. The deviation is sampled over the whole triangle — by default along every edge and across the interior (DEVIATION_SAMPLES) — against the barycentric interpolation of the corner values, and excess_levels is log2(deviation / tol): the deviation of a smooth function under linear interpolation is O(h²) and bisection halves an edge every second level, so each level buys a factor 2. f may return points (geometry error: pass the base's mapping) or scalars (solution error: pass the colour evaluation); tol is absolute, in the units of norm of f's values.

samples is the set of barycentric points the deviation is measured at, each a weight triple over the key's corners. A sampled deviation is a lower bound: finitely many points can miss the peak of a high-order field between them (four samples resolve the quadratic deviation profile exactly, but from cubic order on the true maximum can fall between the sampled points), so pass a denser set to tighten the bound — for diagnostics, or when a high-order field under-refines. Each triple must be nonnegative and sum to 1.

The deviation of a key is a pure function of f and the sample set — it contains neither the tolerance nor any refinement state — so it can be memoized across refinement calls, plots and tolerance changes for as long as f does not change. Pass a Dict{UInt64,Float64} as cache to do so; the caller owns the dict and is responsible for emptying it when f's underlying data changes (the adaptive plots key this to the substrate's solution epoch), and must not share it between criteria with different sample sets. Without a cache every query samples f afresh — orders of magnitude more expensive than a cache hit.

Sampling the interior is what makes the criterion bound what is actually drawn — for curved geometry the deviation peaks in the middle of a face, and an edge-only criterion happily leaves it there. Keeping the edge samples as well bounds the width of any gap or colour seam at a refinement-level boundary, which matters when conformity is disabled.

The criterion is not monotone in depth (a triangle's deviation can vanish while a descendant's does not — e.g. a bilinear field is linear along a quad's outer edges but curved along the fan diagonals), which makes the refined state mildly path-dependent: a state merged down from finer keys may stay finer than one refined up from the roots, because the passes never discard detail whose deviation still exceeds the tolerance. The finer of the two states is the more accurate one.

deviation(lod::DeviationLoD, base, key) -> Float64

The sampled deviation of key's linear interpolation from lod.f, memoized in lod.cache when one is attached. This is the expensive half of excess_levels; the tolerance comparison on top of it is free.

CombinedLoD(lods...)

Split when any member criterion wants to: the excess is the member maximum.

CachedLoD(inner)

Memoize a criterion per key. For a fixed solution the excess is a pure function of the key, but the passes ask for the same keys again and again — every refinement round re-tests the surviving leaves, and the conforming closure additionally asks about parents and neighbours. With an FE evaluation behind every sample point that repetition dominates; refine_keys! therefore wraps its criterion in this for the duration of the call.

excess_levels(lod, base, key) -> Float64

The extension point of a level-of-detail criterion: how many bisection levels below key does the target resolution lie? Positive means the key is too coarse (split), non-positive for the key's parent means the parent is already fine enough (its children merge). For crack-free results the value must be a symmetric function of the key's split edge, varying by less than one level between adjacent triangles (see the file header).

PolyBasis(ip) -> PolyBasis or nothing

The change of basis from an interpolation's nodal values to monomial coefficients, or nothing when the interpolation is not a nodal polynomial space this can represent. Built once per interpolation type.

PolyField(basis, ncells, T) -> PolyField

Per-cell monomial coefficients of one field. refresh! recomputes them from the current nodal values; evaluate then costs N multiply-adds.

key_neighbour(base, key, edge) -> (key, edge, reversed) | nothing

Which triangle at key's own level lies across edge (EDGE_S/EDGE_L/EDGE_R), which of its edges that is, and whether the two traverse the shared segment in opposite directions — or nothing when the edge is a boundary.

Answered purely from the key and the base adjacency table, by recursing to the root:

child 0 = (v1, m, v2)   split edge = parent's left leg
child 1 = (v2, m, v3)   split edge = parent's right leg
both     legs           = halves of the parent's split edge, or the new
                          interior edge shared by the two children

Well-definedness rests on the compatibility invariant — split edges pair only with split edges, legs only with legs — which this recursion preserves at every level, and which the base table is built to satisfy. A GPU port replaces the base lookup with the same recursion over a concurrent binary tree; nothing here needs the leaf set to be materialized.

diamond_partner(base, key) -> UInt64 | nothing

The triangle at key's own level across its split edge — the only leaf that has to split together with key to keep the mesh conforming.

force_split!(leaves, base, key; max_depth)

Split key and everything that has to split with it to keep the mesh conforming: the leaf across its split edge, recursively forced down when it is coarser. Terminates because each forced neighbour is strictly coarser than the triangle that asked for it.

conforming_update!(leaves, base, lod; max_depth, hysteresis) -> changed::Bool

One conforming split/merge pass over the leaf set: every leaf whose excess_levels is positive is split through force_split!, and every diamond of four leaves whose two parents both want to coarsen is merged back. Both operations move whole diamonds, which is what keeps every drawn edge a full edge of the leaf on the other side — no T-vertices, and therefore no slivers or colour seams between refinement levels, on curved geometry as much as on flat.

leb_order(corners::NTuple{3})

Rotate a corner triple so its longest edge connects the first and third corner — the edge update_keys! bisects. Cyclic, so orientation is preserved.

ReferenceTessellation{refdim,T}

Surface triangulation of a reference shape: coords are vertices in reference space, triangles index into coords. Coordinates may be shared between triangles; the per-cell vertex duplication that makes discontinuous (L2) fields render correctly is applied by FEData, not here.

edges are the wireframe segments drawn by meshplot: polylines along the finite element cell's edges (from Ferrite.reference_edges), also indexing into coords. They are deliberately separate from the triangles — the triangulation contains interior diagonals (e.g. the quadrilateral's center fan) that are not cell edges and must not show up in the wireframe. A tessellation without edges renders no wireframe.

reference_tessellation(::Type{<:Ferrite.AbstractRefShape}) -> ReferenceTessellation

The extension point for custom cell types: return the surface triangulation of the reference shape. Implementing this one method makes FEData and all representations work for cells with that reference shape. For 3D shapes, facet_based_tessellation builds a valid tessellation from Ferrite.reference_faces.

facet_based_tessellation(::Type{<:Ferrite.AbstractRefShape{3}}) -> ReferenceTessellation

Build the surface tessellation of a 3D reference shape from its faces: each face's 2D tessellation (triangle or quadrilateral, chosen by vertex count) is mapped into the element via Ferrite.facet_to_element_transformation. The wireframe edges come from Ferrite.reference_edges, with their own (duplicated) endpoint vertices appended after the face vertices. This is the default building block for reference_tessellation of volumetric shapes.

subdivide(tess::ReferenceTessellation, n::Int) -> ReferenceTessellation

Subdivide a tessellation in reference space, n times: every triangle into 4 (orientation preserving) and every edge segment into 2. Midpoints are deduplicated by coordinate, so subdivided triangles share vertices with their neighbours and with edge segments running along the same reference line. New vertices are appended, so indices into the input tessellation stay valid.

This is what resolves curved geometry: the subdivided reference vertices are mapped through the cell's geometric interpolation when the tessellation is instantiated by FEData, so surfaces and wireframe edges of high-order (or nonlinearly deformed) cells bend instead of being drawn as flat facets and straight chords.

QPTessellation{refdim}

Voronoi partition of a reference shape induced by a quadrature rule. coords are vertices in reference space, triangles index into them, and vertex_qp[v] is the quadrature point whose region vertex v belongs to.

Vertices are not shared between regions: every region carries its own copy, so assigning each vertex its quadrature point's value renders the region flat.

qp_voronoi_tessellation(::Type{<:Ferrite.AbstractRefShape}, qr::Ferrite.QuadratureRule) -> QPTessellation

Partition a reference shape into the Voronoi regions of qr's quadrature points and triangulate them. For 3D shapes the partition is intersected with the boundary faces, which is what the surface renderer draws.

Number of triangles a cell tessellates into.

Total number of tessellation vertices, i.e. vertices of the rendered triangulation. These are not the vertices of the finite element cells: cells do not share them (they are duplicated per cell, so discontinuities render), and a tessellated cell generally carries more of them than it has corners.

transfer_solution(ds::FEData, u::Vector; field_name=:u) -> Matrix{Float64}

Evaluate the field at every tessellation vertex of every visible cell from the owning element's dofs (preserving inter-element discontinuities). Vertices of invisible cells or cells outside the field's subdomain stay NaN.

transfer_scalar_celldata(ds::FEData, values::AbstractVector) -> Vector{Float64}

Expand one scalar per cell to the tessellation vertices.

interpolate_gradient_field(dh::DofHandler, u::AbstractVector, field_name::Symbol; copy_fields::Vector{Symbol})

Compute the piecewise discontinuous gradient field for field_name. Returns the flux dof handler and the corresponding flux dof values. If the additional keyword argument copy_fields is provided with a non empty Vector{Symbol}, the corresponding fields of dh will be copied into the returned flux dof handler and flux dof value vector.

_tensorsjl_gradient_accessor(v::Tensors.Vec, field_dim_idx::Int, spatial_dim_idx::Int)

This is a helper to access the correct value in Tensors.jl entities, because the gradient index is the outermost one.