Adaptive mesh refinement (AMR)

isancestor(o1, o2, b) -> Bool

Is o1 a strict ancestor of o2, i.e. is o2 a descendant of o1? Walks o2's parent chain up to the root, which is itself a valid ancestor (parent of the root is the root, hence the level-0 stop instead of a levels-remaining counter).

From Burstedde et al. [15];

The octant coordinates are stored as integers of a fixed number b of bits, where the highest (leftmost) bit represents the first vertical level of the octree (counting the root as level zero), the second highest bit the second level of the octree, and so on.

A morton index can thus be constructed by interleaving the integer bits (2D): $m(\text{Oct}) := (y_b,x_b,y_{b-1},x_{b-1},...y_0,x_0)_2$ further we assume the following

Due to the two-complement representation of integers in practically all current hardware, where the highest digit denotes the negated appropriate power of two, bitwise operations as used, for example, in Algorithm 1 yield the correct result even for negative coordinates.

also from Burstedde et al. [15]

children(octant::OctantBWG{dim, N, T}, b::Integer) -> NTuple{N, OctantBWG}

Compute the N = 2^dim children of octant, returned in z-order (x before y before z).

The first child's vertices of octant are utilized. Its 2^dim vertices coincide exactly with the anchors (lower-left corners) of all children, so each child is simply the level-l+1 octant placed at the corresponding vertex of family.

In 2D, with parent anchor ⊙ at (x,y) and edge length H = 2h:

(x,y+H) +───────┬───────+ (x+H,y+H)        │  c3   │  c4   │            family = first child (c1), anchored at (x,y),(x,y+h) ├───────┼───────┤            edge length h. Its four vertices        │  c1   │  c2   │              v1=(x,  y  )  v2=(x+h, y  )  (x,y) ⊙───────┴───────+ (x+H,y)      v3=(x,  y+h)  v4=(x+h, y+h)           (x+h,y)                 are exactly the anchors of c1..c4.

vertices(octant::OctantBWG{dim}, b::Integer)

Computes all vertices of a given octant. Each vertex is encoded within the octree coordinates i.e. by integers.

edges(octant::OctantBWG{dim}, b::Integer)

Computes all edges of a given octant. Each edge is encoded within the octree coordinates i.e. by integers. Further, each edge consists of two three-dimensional integer coordinates.

faces(octant::OctantBWG{dim}, b::Integer)

Computes all faces of a given octant. Each face is encoded within the octree coordinates i.e. by integers. Further, each face consists of either two two-dimensional integer coordinates or four three-dimensional integer coordinates.

corner_neighbor(octant::OctantBWG, c::Integer, b::Integer)

Computes the corner neighbor octant which is only connected by the corner c to octant

edge_neighbor(octant::OctantBWG, e::Integer, b::Integer)

Computes the edge neighbor octant which is only connected by the edge e to octant.

facet_neighbor(octant::OctantBWG{dim, N, T}, f::T, b::T = DEFAULT_MAXLEVEL[dim]) -> OctantBWG{dim, N, T}

Intraoctree face neighbor for a given faceindex f (in p4est, i.e. z order convention) and specified maximum refinement level b. Implements Algorithm 5 of Burstedde et al. [15].

x-------x-------x|       |       ||   3   |   4   ||       |       |x-------x-------x|       |       |o   1   *   2   ||       |       |x-------x-------x

Consider octant 1 at xyz=(0,0), a maximum refinement level of 1 and faceindex 2 (marked as *). Then, the computed face neighbor will be octant 2 with xyz=(1,0). Note that the function is not sensitive in terms of leaving the octree boundaries. For the above example, a query for face index 1 (marked as o) will return an octant outside of the octree with xyz=(-1,0).

possibleneighbors(o::OctantBWG{2}, l, b)

Returns a tuple of possible neighbors, where the first four are corner neighbors that are exclusively connected via a corner. The other four possible neighbors are face neighbors. Always returns the full NTuple (type-stable); callers that want only in-tree neighbours filter with inside.

possibleneighbors(o::OctantBWG{3}, l, b)

Returns a tuple of possible neighbors, where the first eight are corner neighbors that are exclusively connected via a corner. After the first eight corner neighbors, the 6 possible face neighbors follow and after them, the edge neighbors. Always returns the full NTuple (type-stable); callers that want only in-tree neighbours filter with inside.

transform_corner(forest, k, c, oct, inside::Bool)
transform_corner(forest, v::VertexIndex, oct, inside::Bool)

Algorithm 12 but with flipped logic in Burstedde et al. [15] to transform corner into different octree coordinate system Implements flipped logic in the sense of pushing the Octant oct through vertex v and stays within octree coordinate system k.

c is the corner of tree k (in BWG corner numbering) at the shared vertex; oct is placed at that corner of k, inside the root (inside = true) or diagonally outside of it (inside = false). A corner octant is fully determined by the corner index and the level, so no connectivity lookup is needed — in particular the corner must not be re-derived from vertex_vertex_neighbor[k, ...][1], which is ambiguous (and wrong) as soon as more than two trees meet at the vertex.

transform_edge(forest, k, e, k′, e′, oct, inside::Bool)
transform_edge(forest, k′, e′, oct, inside::Bool)
transform_edge(forest, e′::EdgeIndex, oct, inside::Bool)

Algorithm 10 in Burstedde et al. [15] to transform an edge into a different octree coordinate system, but with reversed logic. See transform_edge_remote with logic from paper. Transform the octant oct, which sits at edge e of pivot tree k (in k's coordinates), into the coordinate system of the neighbouring tree k′ at its edge e′. The along-edge coordinate is taken from the pivot's edge axis and mirrored iff trees k and k′ traverse the shared macro edge in opposite directions.

Both the pivot pair (k, e) and the target pair (k′, e′) must be passed explicitly: a macro edge can be shared by more than two trees, so neither the pivot nor the relative orientation can be re-derived from edge_edge_neighbor[..][1] lookups (those pick an arbitrary incident tree). The two shorter forms assume exactly two trees at the macro edge and take edge_edge_neighbor[k′, e′][1] as the pivot.

transform_facet(forest::ForestBWG, k', f', o::OctantBWG) -> OctantBWG
transform_facet(forest::ForestBWG, f'::FacetIndex, o::OctantBWG) -> OctantBWG

Interoctree coordinate transformation of an given octant o that lies outside of the pivot octree k, namely in neighbor octree k'. However, the coordinate of o is given in octree coordinates of k. Thus, this algorithm implements the transformation of the octree coordinates of o into the octree coordinates of k'. Useful in order to check whether or not a possible neighbor exists in a neighboring octree. Implements Algorithm 8 of Burstedde et al. [15].

x-------x-------x|       |       ||   3   |   4   ||       |       |x-------x-------x|       |       ||   1   *   2   ||       |       |x-------x-------x

Consider 4 octrees with a single leaf each and a maximum refinement level of 1 This function transforms octant 1 into the coordinate system of octant 2 by specifying k=1 and f=2. While from the perspective of octree coordinates k=2 octant 1 is at xyz=(-2,0), the returned and transformed octant is located at xyz=(0,0)

transform_corner_remote(forest, k, c, oct, inside::Bool)
transform_corner_remote(forest, v::VertexIndex, oct, inside::Bool)

Algorithm 12 in Burstedde et al. [15] to transform corner into different octree coordinate system. Follows exactly the version of the paper by taking oct and looking from the neighbor octree coordinate system (neighboring to k,v) at oct.

transform_edge_remote(forest, k, e, oct, inside::Bool)
transform_edge_remote(forest, e::EdgeIndex, oct, inside::Bool)

Algorithm 10 in Burstedde et al. [15] to transform edge into different octree coordinate system. This function looks at the octant from the octree coordinate system of the neighbor that can be found at (k,e)

transform_facet_remote(forest::ForestBWG, k::T1, f::T1, o::OctantBWG{dim, N, T2}) -> OctantBWG{dim, N, T2}
transform_facet_remote(forest::ForestBWG, f::FacetIndex, o::OctantBWG{dim, N, T2}) -> OctantBWG{dim, N, T2}

Interoctree coordinate transformation of an given octant o to the face-neighboring of octree k by virtually pushing os coordinate system through ks face f. Implements Algorithm 8 of Burstedde et al. [15].

x-------x-------x|       |       ||   3   |   4   ||       |       |x-------x-------x|       |       ||   1   *   2   ||       |       |x-------x-------x

Consider 4 octrees with a single leaf each and a maximum refinement level of 1 This function transforms octant 1 into the coordinate system of octant 2 by specifying k=2 and f=1. While in the own octree coordinate system octant 1 is at xyz=(0,0), the returned and transformed octant is located at xyz=(-2,0)

refine_octant!(octree::OctreeBWG, pivot_octant::OctantBWG)

Internal, octree-level refinement primitive; the user-facing entry point is refine!(forest, cellids).

Replace the leaf pivot_octant in octree.leaves by its 2^dim children, spliced in z-order into the parent's slot so the Morton order of leaves is preserved. A no-op if pivot_octant is already at the tree's maximum level octree.b.

pivot_octant is located with a searchsortedfirst binary search, which requires octree.leaves to be Morton-sorted (the Base.isless total order) — the standard invariant for a BWG octree. A single call is O(n) because of the in-place insert!; to refine many leaves at once prefer refine_all! or the refine!(forest, cellids) vector method, which rebuild the leaf list in one linear pass instead of n shifts.

coarsen_octant!(octree::OctreeBWG, o::OctantBWG)

Internal, octree-level coarsening primitive; the user-facing entry point is coarsen!(forest, cellids).

Replace the 2^dim-sibling family that o belongs to with their common parent in the tree's Morton-sorted leaves. o is snapped back to the family's first sibling (via its child_id/morton), the parent is written in its slot, and the remaining 2^dim - 1 siblings are deleted — the inverse of refine_octant!. Assumes the whole family is present and at the same level (e.g. after balanceforest!).

_coarsen_all!(forest::ForestBWG)

Internal convenience for tests and development — not part of the public API.

Coarsen every 2^dim-sibling family in forest by one level — the inverse of refine_all!. Each leaf that is a first sibling (child_id == 1) is replaced by its parent via coarsen!.

_balance_leaf!(forest, k, tree, o, perm_face, perm_face_inv, perm_corner, perm_corner_inv, rootfaces, rootedges, rootvertices, facet_neighborhood)

Per-leaf kernel of balanceforest! handling the inter-tree part of the 2:1 balance.

Operates on a single "pivot" leaf o of tree k. In-tree balancing is already taken care of by balancetree; this function only propagates balance across tree boundaries. It walks the possibleneighbors of o, keeps those lying outside the current tree (reachable only through a corner/face/edge connection to another tree), decodes the neighbour type from the possibleneighbors index s_i, maps the pivot's local index into the neighbour tree via the permutation tables, and calls balance_face/balance_corner/balance_edge to refine the neighbour tree where the balance condition requires it.

_touches_tree_boundary(o::OctantBWG{dim}, b) -> Bool

true iff octant o has a face on its tree's boundary (some axis anchor at 0 or at the root extent 2^b). Only such leaves can have out-of-tree neighbours, so balanceforest! uses this to skip the interior leaves (the majority) when balancing across tree interfaces.

inside(oct::OctantBWG{dim}, b) -> Bool

Whether oct lies within its tree's root domain [0, 2^b)^dim (see _maximum_size). A false means the octant has crossed a tree boundary — the signal that an inter-tree transform (transform_facet/transform_corner) is needed to express it in the neighbouring tree's coordinate system.

_maximum_size(b::Integer) -> Int

Integer extent 2^b of the root octant along each axis (computed as 1 << b): the octree coordinate domain of a tree is [0, 2^b)^dim. Used by inside to detect when an octant leaves its tree.

IteratePoint{dim}

A point in the sense of Isaac et al. [16] §2.1: an octant volume or one of its lower-dimensional features (face, 3D edge, corner), encoded topologically by the box anchor (minimum integer corner), level (box extent h = _compute_size(b, level)) and axes (the directions the box extends along). point_dim(c) = count(axes) is the paper's dim(c); equality is field-wise — no physical coordinates, no rounding.

iterate_points(visit, tree::OctreeBWG, sc::IterScratch; mindim = 0, maxdim = dim, skip_conforming = false)

Isaac et al. [16] Algorithm 5.3 (Iterate), serial: drive _iterate_interior! from the closure of each tree root. visit(c::IteratePoint, leaf_supp::LeafSupport) is called once for every point c ∈ PΩ (5.1) — every non-hanging volume / face / edge / corner — with leaf_supp the leaves surrounding it plus their leaf indices (§6.4: "Iterate provides the index"). Use point_dim(c) to dispatch per dimension (volume = dim, face = dim-1, edge = 1, corner = 0), or pass mindim/maxdim for the §5.4 specialization (e.g. mindim = dim - 1 to visit only volumes + faces, or maxdim = dim - 1 to skip the volume callback like a NULL volume callback in p4est_iterate). With skip_conforming = true, face/edge points whose supports are all leaves of equal level (conforming interfaces) are skipped as well — the specialization for callbacks that only act on non-conforming interfaces, like the hanging-node detection; corner points always fire. The descent is allocation-free: callers traversing many trees (creategrid, facetskeleton) allocate one IterScratch and pass it to every tree of equal maximum level b (the buffers are depth-indexed, so they only depend on b); leaf_supp wraps reused buffers — copy what you retain.

Looping the trees of a forest is the serial Alg 5.3. Within a tree this visits PΩ exactly; at shared tree boundaries a feature is currently visited once per incident tree (its per-tree leaf_supp covers only that tree's leaves) — creategrid reconciles the per-tree visits through its boundary tables. Cross-tree coordinated descent (single-visit boundary leaf_supp via the orientation transforms, fully mirroring p4est_iterate) is the documented next step, consistent with _iterate_interface_hanging!'s inter-tree face descent.

_iterate_interior!(visit, c::IteratePoint, depth, sc::IterScratch, leaves, b, mindim, maxdim, skipconf)

Isaac et al. [16] Algorithm 5.2 (Iterate_interior), serial and allocation-free. On entry sc.supp[depth] holds c's support octants (eq 2.11) and sc.S[depth][i] the leaves index range under each. Exactly as in Alg 5.2, the recursion stops when some support octant is itself a leaf (line 7; for a corner point always, lines 16-18) — so hanging points are never visited (PΩ, eq 5.1) — and calls visit(c, leaf_supp) with the LeafSupport built via _child_touches_point (line 14) / _descend_to_corner (line 18); otherwise it descends part(c), slicing each range with split_bounds. leaf_supp wraps reused buffers — copy if retained. mindim/maxdim specialize the callback (§5.4): dims below mindim are not recursed into, and visit fires only for point_dim(c) ∈ mindim:maxdim. skipconf additionally suppresses the visit at conforming (all-supports-leaf) interfaces, see the comment inline.

IterScratch{N, M, OT}
IterScratch(tree::OctreeBWG)

Preallocated per-depth working memory of the recursive descent — the sc argument of iterate_points — so the traversal allocates nothing (Isaac et al. [16] §5.4). In a DFS only one root-to-node path is live, so buffers are indexed by recursion depth and reused across siblings; see the field comments for what each buffer holds. M = N + 1 is the split_bounds tuple length. Sizes depend only on the maximum level b, so one scratch is shared across all trees of a forest (creategrid does exactly that).

_foreach_partc(f, c::IteratePoint, b)

Call f(e, combo) for each e ∈ part(c), the partition of c (eq 2.7 of Isaac et al. [16]): the 3^dim(c) one-level-finer points strictly interior to dom(c) — per extending axis the lower half, the degenerate mid plane, or the upper half (boundary features belong to the parent point's partition, which is what makes the descent visit every entity exactly once). combo is the base-3 slot encoding (one digit per extending axis, ascending), the key into the precomputed _part_mask tables.

_foreach_root_closure(f, ::Val{dim}, b)

Call f(c) for each point in the closure of the tree root (Alg 5.3 line 4 of Isaac et al. [16], single tree): the root volume and all its boundary faces/edges/corners — the 3^dim seeds of the recursive descent. The root's boundary features have no parent split to produce them (unlike interior features, which arise as part of an ancestor), so they get their own descent seeds here. Along each axis the feature is pinned to the low face (coord 0), spans the full root, or is pinned to the high face.

_child_touches_point(ch::OctantBWG, c::IteratePoint, b) -> Bool

Whether point c (a feature of ch's parent) lies in the closure of child octant ch — the child-boundary-intersection set B_∩^j (eq 4.5 / Alg 5.2 line 14 of Isaac et al. [16]): true iff ch's box straddles c's coordinate along every degenerate axis of c.

_descend_to_corner(c::IteratePoint, s::OctantBWG, lo, hi, leaves, b) -> (leaf, index)

Find the leaf under support octant s (leaves[lo:hi]) whose closure contains the 0-point c — the atom supp(c) search (Alg 5.2 line 18 / Prop 2.8 of Isaac et al. [16]) — returning it with its index into leaves (the paper's element index j, §6.4). c is a corner of s, so the descent follows the fixed ci-most path, narrowed with split_bounds; the Morton-first/-last slots resolve in O(1).

split_bounds(leaves, lo, hi, a::OctantBWG, b) -> NTuple{2^dim + 1, Int}

Algorithm 3.3 of Isaac et al. [16]. Given the contiguous, Morton-sorted leaf sub-range leaves[lo:hi] (all strict descendants of a), return boundary indices k such that child i of a occupies leaves[k[i]:k[i+1]-1]. Non-allocating: returns a stack NTuple — no 𝐤 vector and no SubArray views, so the recursive descent that calls it at every internal octant stays allocation-free. Short ranges (the vast majority: the recursion halves the range per level) are split by a single linear ancestor_id sweep; long ranges by 2^dim - 1 binary searches (the child key ancestor_id is monotone along the range).

Shared descent helper of the point iterator (_iterate_interior!, _descend_to_corner) and the interface descent (_iter_interface!).

LeafSupport{OT}

The local leaf support set leaf_supp(c) handed to the iterate_points callback: the leaves whose closure touches the visited point c, each paired with its index into the tree's Morton-sorted leaves (the element index j of Isaac et al. [16] §6.4). Iterating yields the octants; ls.idxs[i] belongs to ls.octs[i]; an index of 0 marks a non-leaf entry — impossible on a 2:1-balanced forest, and turned into an error by creategrid. Both vectors are reused buffers of the traversal — copy them if you retain them past the callback.

LnodesVisitor{dim, N, T, TI}

Per-tree numbering context of the creategrid traversal — the Lnodes_callback of Isaac et al. [16] Alg 6.2, as a callable struct (a concrete top-level type instead of a closure, so the captures don't box). One instance per tree: E, the provisional node data (nodecoords_prov, cnt) and the constraint records (cons2/cons4) are shared across trees; bnd, the geometry and the element offset are the tree's own.

Callback dispatch by point_dim(c): corners create + scatter node ids (_visit_corner!); non-conforming faces and (3D) edges create the hanging vertex and record its constraint (_visit_face!, _visit_edge3d!); volumes need no work (cell connectivity IS the filled E).

_visit_corner!(v::LnodesVisitor, c, ls)

Corner callback (dim(c) == 0): c is a non-hanging mesh vertex — hanging points are never visited by the iterator — so every supporting leaf has c as a corner. Create its node and scatter the id into each leaf's E slot (the leaf indices come with the support set).

_mixed_support(c::IteratePoint, ls::LeafSupport) -> (coarse, fine)

Non-conformity test of a finalized face/edge point (Isaac et al. [16] Fig 5): a leaf support at level(c) (the coarse side) coexisting with finer children. Returns the index of the first coarse leaf in ls (its corners are the constraint masters) and whether finer leaves exist — the interface hangs iff finer leaves coexist with the coarse one. A finalized face/edge always has at least one support that is itself a leaf at level(c) (that is what stopped the recursion in _iterate_interior!), so with all supports equal-level leaves the interface is conforming and there is nothing to do.

_global_numbering(E, alias, nodecoords_prov) -> (final_of_prov, nodecoords)

Serial Global_numbering (Isaac et al. [16] Alg 6.1): sweep the element-node matrix in (element, element-node) lexicographic order — the linear (column-major) order of E — and assign each canonical node its final dense id at first encounter. With the ownership rule owner(c) = min leaf supp(c) (eq 6.2), the first element referencing a node is its owner, so this reproduces the paper's partition-independent numbering without any per-node search. Also gathers the final node coordinates (the canonical provisional node's) in final-id order.

A zero entry in E means some element vertex was never assigned a node — impossible on a 2:1-balanced forest — and raises an error.

_merge_intertree_nodes!(forest::ForestBWG{dim}, bnd, alias)

Identify nodes shared across tree boundaries (creategrid's cross-tree pass). For each tree k, walk its root-vertex, root-face and (3D) root-edge neighbours; a node on a shared boundary is matched to its image in the lower-index neighbour k′ via transform_facet/ transform_corner/transform_edge (handling tree rotations), and aliased onto that owner. Only the lower-index tree owns a shared node (k > k′), giving a single canonical id per geometric node across all incident trees.

Node ids are looked up in bnd, the per-tree boundary node tables: bnd[k] holds (packed coord, provisional id) for every node of tree k lying on its root boundary, sorted by key (filled by the numbering traversal, see creategrid). This is the only node-lookup structure of the whole materializer — O(surface) per tree instead of a global coordinate hash map — and the walk itself visits only boundary leaves. The canonicalization is recorded in alias::Vector{Int} (indexed by provisional id, identity-initialized): alias[p] is the provisional id that p is merged onto, so the per-node canonical lookup in creategrid is an array index.

_bnd_lookup(bnd, coord::NTuple{dim, <:Integer}, b) -> Int

O(log) lookup in a tree's sorted boundary node table bnd ((packed coord, provisional id) pairs, see _merge_intertree_nodes!): the provisional id of the node at integer octree coord coord, or 0 if coord is out of the tree's [0, 2^b] range or is not a node of that tree. Both happen routinely during the cross-tree merge: the transform_facet/transform_edge images can land outside the neighbour's domain, and a hanging node exists as a leaf vertex on the refined side of an interface only.

_build_cells(::Type{CT}, E, node_map, final_of_prov, ::Val{NV}) -> Vector{CT}

Materialize the cell vector from the element-node matrix E (provisional ids, z-order slots): column gid is cell gid's connectivity, remapped through final_of_prov to the final node ids and reordered to Ferrite's vertex order via node_map, wrapped in cell type CT (Quadrilateral/Hexahedron). A top-level function barrier so the cell construction compiles concretely (building cells in the type-unstable creategrid body boxes every cell).

reconstruct_facetsets(forest::ForestBWG{dim}) -> Dict{String, OrderedSet{FacetIndex}}

Transfer the macro-mesh facet sets onto the materialized (refined) grid. For each original FacetIndex (tree, face), emit a FacetIndex for every leaf of that tree lying on the root face, converting between p4est and Ferrite face ordering (𝒱₂_perm/𝒱₃_perm). This keeps named boundaries (e.g. Dirichlet/Neumann sets) valid after refinement.

Staying inside one tree there is no rotation, so a leaf is on the root face iff its anchor lies on that face's axis-aligned plane (leaf.xyz[axis] == 0 for a low face, == 2^b - leafsize for a high face), and the contributing local face index is exactly the root face index. This is an O(#leaves) plane test, replacing a former O(#leaves · 2dim) contains_facet scan over each leaf's faces.

reconstruct_cellsets(forest::ForestBWG) -> Dict{String, OrderedSet{Int}}

Transfer the macro-mesh cell sets onto the materialized (refined) grid: every leaf inherits the set membership of its tree (macro cell), so each macro cell id in a set is replaced by the cell ids of that tree's leaves (contiguous by _element_offsets). This keeps named subdomains (e.g. material regions) valid after refinement.

_treecorners(forest::ForestBWG{dim}, k::Integer) -> NTuple{2^dim, Vec{dim}}

Physical coordinates of macro-tree k's 2^dim corner nodes, in Ferrite's vertex order for the tree's cell. These are the interpolation support points for _interp_treepoint; indexing forest.nodes through forest.cells[k].nodes directly keeps the result concrete and allocation-free.

_interp_treepoint(corners::NTuple{N, Vec{dim}}, b, vertex::NTuple{dim}) -> Vec{dim}

Map an integer octree coordinate vertex of a tree to physical space — the isoparametric $Q_1$ geometry map of the macro element (tree). Two steps:

1. affine-scale the octree coordinate (in [0, 2^b]^dim, see _maximum_size) to the reference cube $\xi \in [-1,1]^{dim}$ via $\xi = \texttt{vertex} \cdot 2/2^b - 1$;
2. interpolate the tree's physical corners with the bi-/trilinear Lagrange shape functions, $x = \sum_{j=1}^{N} N_j(\xi)\, \texttt{corners}[j]$.

corners are the tree's 2^dim physical corner nodes (see _treecorners), passed in explicitly so the per-tree corners are computed once and reused for every node of the tree. This is the single bridge from the integer/topological octree world into physical coordinates.

_visit_face!(v::LnodesVisitor, c, ls)

Face callback (dim(c) == dim-1). Conforming faces need no work. On a non-conforming face — a coarse leaf on one side, its refined neighbour's children on the other — the face's interior vertices hang: in 2D the face midpoint (constrained by the 2 endpoints), in 3D the face center (constrained by the 4 face corners). The 3D face's edge midpoints are not emitted here: each hanging edge is its own iterator point, visited exactly once even when shared by several non-conforming faces, and handled by _visit_edge3d!.

_visit_edge3d!(v::LnodesVisitor, c, ls)

Edge callback (3D, dim(c) == 1): the midpoint of a non-conforming coarse edge hangs, constrained by the edge's endpoints. On a 2:1-balanced forest every hanging vertex is either such an edge midpoint or a face center (_visit_face!) — level jumps are capped at one, so no ¼-points exist — and each is created by exactly one callback: the edge midpoint belongs to exactly one coarse edge (two distinct octant edges cannot share their midpoints), even when that edge borders several non-conforming faces.

_iter_interface!(
    cons2, cons4, E, offR, forest, kL, lvsL, octL, loL, hiL, fL,
    kR, lvsR, octR, loR, hiR, fR, bL, bR
)

Synchronized two-sided descent of a shared tree face, emitting inter-tree hanging-node constraints. octL ∈ tree kL (leaves lvsL[loL:hiL], native frame) and octR ∈ tree kR are images of each other across the shared face — fL/fR are the local face indices toward it — and descend in lock-step at equal levels.

• both sides leaves of equal size → conforming, nothing emitted;
• one side a leaf, the other refined → the leaf is the coarse side and the hanging nodes lie on the refined side's face in the fine tree's frame (genuine fine-leaf vertices), emitted via _emit_interface_face! as (element, slot) references into E (offR is tree kR's element offset). Only the kL-coarse case emits here; the kR-coarse case is emitted when the descent is run from (kR, fR), so each interface is handled once per direction;
• both refined → recurse, matching child i on the kL side to its image on the kR side via transform_facet (the validated cross-tree orientation pattern — no new logic).

The integer/topological analogue, across trees, of the intra-tree face/edge callbacks of the numbering traversal (see creategrid).

_emit_interface_face!(cons2, cons4, E, off, lvs, lo, hi, octR, fR, b)

The kL side of an interface octant pair is a leaf while the kR side octR (leaves lvs[lo:hi], element offset off) is refined: the interior points of the shared face hang. Every one of them is a vertex of octR's face children — leaves under the 2:1 balance, already numbered by tree kR's own traversal — so the constrained ids are read straight off E and the constraints recorded exactly like the intra-tree emitters: the face midpoint (2D) / face center (3D) constrained by the face corners, and in 3D each face-edge midpoint constrained by its edge's endpoints, masters in ascending coordinate order. (On a >1 level jump the ids are read from the corner leaves via _subtree_corner_ref's descent.)

_iterate_interface_hanging!(cons2, cons4, E, offsets, forest::ForestBWG)

Collect all inter-tree hanging-node constraints of forest into cons2/cons4. For each shared tree face (found via the facet–facet neighbourhood), seed _iter_interface! with the two tree roots and let it descend both sides in lock-step. The forest-level counterpart of the intra-tree face/edge callbacks; together they capture every hanging node of the materialized grid. Single-tree forests have no shared faces, so this is a no-op there.

center(pivot_face) -> NTuple{dim, Int}

Integer centroid of a set of octree coordinates (the corners of a face or edge): the elementwise sum divided by the number of points. For a coarse face bordering a refined neighbour this is exactly the hanging node — the face center (4 corners) or edge midpoint (2 corners) — in the same integer/topological frame (no physical coordinates).

contains_facet(mface, sface) -> Bool

Whether the sub-facet sface is geometrically contained in the master facet mface, both given as octree corner coordinates. In 2D the facets are axis-aligned segments and containment is an interval test along the shared axis; in 3D mface is a face and sface is accepted iff its center lies inside mface's bounding box. Used to decide whether a leaf face lies on a coarse/root face (refinement interfaces, boundary-set transfer).

FacetSkeletonVisitor{dim}

Per-tree face callback of facetskeleton's intra-tree traversal. At each (non-hanging) face point the LeafSupport holds the leaves of both sides of the face; the visitor recovers each side's local facet index from the face normal (the point's degenerate axis, sides told apart by the anchor comparison of _side_mask) and pushes the facet pair(s): one pair for a conforming face, one pair per fine child — fine side first, the coarse leaf second — for a non-conforming one. Face points on the tree boundary (anchor 0/2^b along the normal) are skipped: their support is one-sided within the tree, and they belong to either the inter-tree descent or the domain boundary.

_iter_interface_facets!(
    skel, perm, offL, offR, forest, kL, lvsL, octL, loL, hiL, fL,
    kR, lvsR, octR, loR, hiR, fR, bL, bR
)

facetskeleton's inter-tree counterpart of _iter_interface!: the same synchronized two-sided descent of a shared tree face, but emitting the facet pairs of the interface instead of hanging-node constraints. Each shared face is descended once per direction, so to emit exactly once a conforming leaf pair fires only from the direction with the smaller (tree, face) key, and a hanging interface only when the coarse side is kL — the fine subfacets then live in tree kR's frame and are enumerated by _emit_interface_facets!, fine side first.

_emit_interface_facets!(skel, perm, off, lvs, lo, hi, oct, f, b, coarse)

One-sided descent of the refined side of a hanging inter-tree interface: enumerate the leaf subfacets on oct's face f (leaves lvs[lo:hi], element offset off) and pair each with the coarse side's FacetIndex coarse, fine side first. Recursing to the leaves instead of stopping at oct's face children tolerates interfaces with a >1 level jump, like _subtree_corner_ref's descent.

_element_offsets(forest::ForestBWG) -> Vector{Int}

Element offset of each tree into the materialized cell vector: tree k's leaf j (Morton order) is cell offsets[k] + j of the grid returned by creategrid.