Forms

GridFunctionScale(vₕ, op::LazyOp)

op scaled pointwise by the grid function or function vₕ, as a symbolic node.

IdentityOperator(Wₕ::AbstractSpaceType)

The identity on Wₕ, as a symbolic node.

LazyOp{D} <: OperatorType

A node of the symbolic operator tree over a D-dimensional space. Records an operation without performing it, so that a form can be written as an expression and assembled later.

OperatorAdd(left::LazyOp, right::LazyOp)

The sum of two symbolic nodes over the same space.

OperatorScale(α, op::LazyOp)

op scaled by the number α, as a symbolic node.

OperatorType

Supertype of everything the form layer treats as an operator.

ZeroOperator(Wₕ::AbstractSpaceType)

The zero operator on Wₕ, as a symbolic node. Absorbs multiplication by a scalar.

is_symbolic(op) -> Bool

Whether op still contains a symbolic placeholder, such as a trial or test function, and so cannot be evaluated until one is substituted.

The base cases are here; src/form/stencil_eval.jl adds the methods for the concrete AST nodes, once every node type exists.

AbsoluteColumn

A stencil entry's trial slot, naming a column of the trial space directly rather than an offset from the point being evaluated.

Every other node's stencil says "this many points from here, on the mesh being walked", which is what lets shift_stencil compose operators by relabelling. An interpolation cannot say that: the trial degrees of freedom it reaches live on a different mesh, and which ones depends on where the point falls (locate_cell). So it names them outright, and the bilinear consumers resolve the two kinds of entry by dispatch.

IndexedTestFunction{D} <: LazyOp{D}

An AST node representing the symbolic test function for a specific component of a composite test space. Carries a runtime component_idx. Used by a coupled form to route stencil contributions to the correct block.

IndexedTrialFunction{D} <: LazyOp{D}

An AST node representing the symbolic trial function for a specific component of a composite trial space. Carries a runtime component_idx identifying which leaf scalar space (1-based, depth-first order) it belongs to. Used by a coupled form to route stencil contributions to the correct block.

PointDependentStencil <: StencilShiftTrait

The operator's stencil depends on where it is evaluated in a way relabelling cannot express, so a neighbour's contribution has to be obtained by evaluating the operator again at the neighbour's own point.

SourceConstant{D, T} <: LazyOp{D}

An AST node representing a source term that is the same number everywhere on the mesh.

SourceFunction reaches this value the general way, through f(point(m, I)): a real cost when f is x -> l, discarding the point it just computed, at every grid point of every assembly. SourceConstant skips point entirely; measured behind a function barrier, assembling a constant source is 1.6–2.6× faster than through SourceFunction, the ratio growing with ndofs rather than staying fixed, so this is a per-point saving rather than one-off overhead. source_number is what builds one from a literal Number.

SourceFunction{D,F} <: LazyOp{D}

An AST node representing a source term defined by a continuous function.

SourceVector{D,VType} <: LazyOp{D}

An AST node representing a source term defined by a discrete vector of values.

Note the division of labour with GridFunctionScale, which also carries values per grid point. A SourceFunction holds a function of position, f(x), evaluated at the point. A Function inside a GridFunctionScale is something else entirely: a zero-argument thunk returning the vector or number to scale by, called as f() both here and in resolve_ast. It defers building that vector until the form is resolved.

So (x -> x[1]) * D₋ₓ(u) does not do what it reads as: the thunk call fails, because the function wants a point. A function of position belongs in a SourceFunction, or should be restricted to the grid with Rₕ first and passed as the vector it becomes.

TestFunction{D, N} <: LazyOp{D}

An AST node representing the symbolic test function $v$ in a form over a D-dimensional space with N components.

TranslationInvariantStencil <: StencilShiftTrait

The operator's stencil has the same shape at every point, so a neighbour's contribution is its own stencil with the offsets relabelled (shift_stencil).

TrialFunction{D, N} <: LazyOp{D}

An AST node representing the symbolic trial function $u$ in a bilinear form over a D-dimensional space with N components.

_in_grid(space, I::CartesianIndex) -> Bool

Whether I names a real point of space's mesh.

The check ShiftNode's own local_stencil makes for a PointDependentStencil inner operator, in place of trusting _clamped_shift's clamp; see the note there for why that trust does not extend to this one caller.

shift_offset(offset::NTuple{D, Int}, dim::Int, delta::Int) -> NTuple{D, Int}

Shifts a Cartesian offset tuple by delta in dimension dim.

shift_stencil(inner::Tuple, ::Val{Dim}, delta)

Shifts all coordinates in a stencil tuple by delta in dimension Dim.

map over a Tuple unrolls and stays type-stable at compile time in Julia — measured against a @generated version this once was (gpena/Bramble.jl#63): identical zero allocations and identical inferred return type, so the code generation bought nothing here.

shifted_inner_stencil(inner_op, inner, space, I, markers, ::Val{Dim}, delta)

The stencil inner_op contributes delta points away in direction Dim, given inner, its stencil already evaluated at I.

The one place the "shift by relabelling" assumption is made, so the one place a node that cannot be relabelled has to be handled: TranslationInvariantStencil relabels inner's offsets and never touches inner_op again, PointDependentStencil discards inner and evaluates inner_op at the shifted point instead. Both produce a tuple of the same static length, since it is the same operator either way, so the callers' concatenate_stencils sees exactly the shape it always did.

source_function(f, ::Val{D}) -> SourceFunction{D, typeof(f)}

Constructs a SourceFunction wrapping function f.

sum_stencil_values(stencil::Tuple)

The sum of a stencil's coefficients, ignoring its offsets entirely.

Required by _contracted_left_stencil (form/operators/inner.jl) for a source-only subtree's own local_stencil: not the offsets, which mean nothing for a value that contributes no matrix structure, only their total. false rather than 0 or zero(T) is the empty-stencil answer: RegionRestriction can legitimately produce () for a point outside its region, and there is no T to call zero on when there are no entries to read one from; false promotes to whatever numeric type the other entries (or, empty, the caller's own multiplication) turn out to have — exactly sum(f, itr; init = false)'s own behavior, which is what this calls. This used to be its own @generated unrolled fold "like every other stencil-algebra primitive" in this file; measured against sum directly (gpena/Bramble.jl#63), identical zero allocations and identical inferred type, so the @generated version bought nothing that sum was not already providing.

test_function(::Val{D}) -> TestFunction{D, nothing}
test_function(space::AbstractSpaceType) -> TestFunction{dim(space), leaf_count(space)}

Constructs a TestFunction of dimension D or for a specific grid space.

trial_function(::Val{D}) -> TrialFunction{D, nothing}
trial_function(space::AbstractSpaceType) -> TrialFunction{dim(space), leaf_count(space)}

Constructs a TrialFunction of dimension D or for a specific grid space.

zero_offset(::Val{D}) -> NTuple{D, Int}

Zero-initialized offset tuple of dimension D.

TappedNode{D, Dim}

The nodes whose stencil is ordered taps from {+1, 0, -1} along Dim with per-node weights: the one-sided and extended differences, the two averages, and the jump. ShiftNode is not one of them – it relabels its child's whole stencil rather than combining taps.

UnaryWrapper{D}

The AST nodes that wrap exactly one operand in an inner_op field.

Thirteen node types, and the reason they are worth naming together: a query about a term usually has the same answer for a wrapper as for the operand inside it, so each such query used to be registered against all thirteen by hand – one line apiece, per query (gpena/Bramble.jl#52). A query that forgot one inherited a fallback instead, and every fallback in this family is a plausible wrong answer rather than an error: test_component_or_nothing answering nothing sends a term to every block.

Written as a union of concrete types rather than reached through a generic child accessor, deliberately. Dispatch resolves it at compile time and each method still reads op.inner_op directly, so nothing on the assembly path pays for the generality.

Products and sums are not members. BilinearProduct and LinearProduct hold two operands with different roles – trial on the left, test on the right – so a query about the test component reads right_op alone, and OperatorAdd has to reconcile both sides. Those stay written out, which is the point: what is left explicit is what genuinely differs.

_is_source_only(op::LazyOp) -> Bool

Whether a LazyOp subtree is source-only: built entirely from sources (SourceFunction/SourceVector) and the plain operators that wrap them, never bottoming out in a TrialFunction/IndexedTrialFunction leaf.

innerₕ's l::Function/l::Number/l::VectorElement overloads (operators/inner.jl) never need this: those three types are never anything but a source, so wrapping them in a LinearProduct is unconditional. The question only exists for an argument that already arrived as a LazyOp: πₕ(uₕ) (interpolate_at) or D₋ₓ(πₕ(uₕ)) are sources too, just already wrapped, and the generic innerₕ(::LazyOp, ::LazyOp) used to build a BilinearProduct regardless, which is the wrong AST shape for a LinearForm's assembly walk: a BilinearProduct's stencil carries a pair of offsets (trial and test), where _scatter_term! (form/linear.jl) expects one.

A missing case defaults to false (the fallback ::LazyOp method below): conservative, since that is exactly the behavior every node had before this predicate existed (always BilinearProduct) for anything not explicitly listed as source-only.

_ast_equal(a, b) -> Bool

Whether two LazyOp subtrees are the same expression: same concrete node type, and every field equal – recursively for a field that is itself a LazyOp, by === otherwise.

=== rather than == for a leaf field (a grid function, a closure, a component index) is deliberate: two arrays that hold equal values right now are not the same operator if one is later mutated in place (Rₕ!(cₕ, ...)) and the other is not, and two independently built closures are never "the same" scaling function even if they happen to compute the same thing. Missing an equal-but-distinct pair only forgoes an optimization; treating two different subtrees as equal would silently change what the assembled form computes, which is the one thing this pass may never do.

simplify_ast(op) -> LazyOp

Rewrite a resolved AST into a form that routes to fewer, cheaper mesh sweeps, and exposes a few patterns the router could not assemble at all, without changing what it computes. See the module comment at the top of this file for the rules.

Every node type not named there is a leaf as far as this pass is concerned and is returned unchanged; form(Wₕ, Vₕ, f)/form(Wₕ, f) call this immediately after resolve_ast.

component(op::LazyOp, i::Int) -> LazyOp

The i-th component of the symbolic operator op: the same expression with its trial and test leaves replaced by their indexed forms.

Reached through the functor, so op(i) is component(op, i). The index distributes, so (v + D₋ₓ(v))(1) and v(1) + D₋ₓ(v(1)) are the same tree.

Block{TrialLeaf, TestLeaf}

One leaf-space pair's rectangle within a composite system matrix: the concrete trial and test leaf spaces a term couples, and the row/column offset each contributes to that rectangle's position in the assembled matrix.

Matrix rows are indexed by the test function (see bilinear.jl's file header), so row_offset always comes from test_leaf's offset in leaf_spaces_offsets, and col_offset from trial_leaf's. That asymmetry used to be carried by convention across six call sites, each unpacking a bare (tc, sc) tuple into first/last calls in the right order – one of which got it backwards (gpena/Bramble.jl#48). Naming it here means a caller reads blk.row_offset/blk.col_offset off the type instead of re-deriving which positional element means which.

_collect_region_labels(op) -> NTuple{N, Symbol}

Every marker label a RegionRestriction anywhere in op names: from restrict_to calls written directly, or from the markers = (...) keyword on innerₕ/inner₊ and friends. Flattened into one tuple; a term naming several restrictions (nested, or one on each side of a product) reports all of them, since every one has to exist on every leaf the term reaches for assembly to mean what it says.

Recurses the same way trial_component_or_nothing/test_component_or_nothing do, so a marker nested behind any operator those already see through is found here too.

_validate_term_markers(term, mesh_markers, context::String)

Throws if term names, via restrict_to or markers = (...), a label that does not exist in mesh_markers (the mesh a term is about to be scattered against). Checked once, while the sparsity pattern is built (allocate_system_matrix/_pattern_term!), rather than left to RegionRestriction's own local_stencil: that answers false for a missing key the same way it does for "not marked", so a typo'd or leaf-missing label would otherwise assemble to a silent all-zero contribution instead of failing loudly.

block_of(term, nblocks_trial, nblocks_test) -> Union{Nothing, Tuple{Int, Int}}

The (trial, test) block term belongs to, or nothing when it belongs to every diagonal block.

A term naming neither side is the same integrand on each block, which for a matrix means the diagonal: Σᵢ innerₕ(uᵢ, vᵢ) is block diagonal, not full. A term naming both is one block. A term naming one and not the other is refused: innerₕ(u(1), v) is not something written in a variational formulation, and reading it as a whole row or column of blocks would be a guess about what was meant.

blocks(term, trial_leaves, test_leaves) -> Tuple{Vararg{Block}}

Every Block term must be assembled into.

A term naming both sides (via block_of) resolves to the one Block it names. A term naming neither is the same integrand on every diagonal block, so it resolves to one Block per diagonal leaf pair. trial_leaves/test_leaves are leaf_spaces_offsets results.

test_component_or_nothing(op) -> Union{Int, Nothing}

The component op is written against, or nothing when it names none.

Used when assembling composite right-hand sides, where a term built from an indexed test function belongs to one block while a term built from an unindexed one belongs to all of them. Written as a query rather than catching an exception because it is evaluated once per term per assembly in time-stepping loops.

trial_component_or_nothing(op) -> Union{Int, Nothing}

The trial component op is written against, or nothing when it names none.

The trial-side mirror of test_component_or_nothing, and needed for the same reason one step further on: a bilinear form's term belongs to a block, which takes a component from each side. A term naming neither is the same integrand in every diagonal block; a term naming both is one block; and a term naming one but not the other is not something the mathematics can express, so it is an error rather than a guess.

_stencil_margin(op) -> Int

The largest |offset| component anywhere in op's stencil, trial and test sides combined.

Unlike stencil_offsets, which keeps only the test-side (row) reach at a BilinearProduct because that is all its one caller (_colour_strides) needs, this keeps both. It exists for visit_bilinear_stencil's interior/boundary split (form/bilinear_traversal.jl, gpena/Bramble.jl#160), which has to guard every offset an entry can carry – trial and test alike, not just the row's – before it can skip the guard anywhere.

Answered as a scalar rather than the offset set stencil_offsets returns: the traversal only needs a margin to build a rectangular interior box from, and a scalar costs no allocation to compute, where reusing stencil_offsets (built on union-ing Vectors) would mean paying that cost inside visit_bilinear_stencil itself on the ReplaySink path, once per assemble! call rather than once per form.

A composed reach past magnitude 1 is not a hypothetical this ignores: D₋ₓ(D₋ₓ(u)) and Shift(u, dim, 2) both carry offsets wider than the one step a single difference or a unit shift would suggest, which is why this walks the same tree stencil_offsets does rather than assuming every node caps out at 1.

stencil_offsets(op) -> Vector{NTuple{D, Int}}

The grid offsets the operator op reaches, sorted and without repeats.

Read from the AST rather than from an evaluated stencil, and the two agree: a truncated point keeps its offsets and zeroes its coefficients, so the set does not vary over the grid. The one node whose reach is not fixed by its type is ShiftNode, which carries its step as a field; the value is available here because this walks the built tree.

For a BilinearProduct this is the row (test-side) reach only, not the full row/column pattern – see the note on that method.

AverageNode{D, Dim}

Either average node over a D-dimensional space, averaging along Dim.

The pair carries the same parameters and differs only in which neighbour it reaches, so anything reading the direction rather than choosing a stencil is written against this alias.

BackwardAverage{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a backward spatial averaging operator acting in dimension Dim.

ForwardAverage{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a forward spatial averaging operator acting in dimension Dim.

ShiftNode{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a stencil shift operation by shift_amount grid points in dimension Dim.

avg_backward(op::LazyOp{D}, dim::Int) where D

Applies a backward average operator to op in dimension dim.

avg_forward(op::LazyOp{D}, dim::Int) where D

Applies a forward average operator to op in dimension dim.

shift_op(op::LazyOp{D}, dim::Int, amount::Int) where D

Shifts the stencil of op by amount grid points in dimension dim.

vectorial_avg_backward(op::LazyOp{D}) where D

Applies backward spatial averaging component-wise across all dimensions.

vectorial_avg_forward(op::LazyOp{D}) where D

Applies forward spatial averaging component-wise across all dimensions.

DifferenceNode{D, Dim}

Either one-sided difference node over a D-dimensional space, differencing along Dim.

The two carry the same parameters, so anything that reads only the direction off the node is written against this alias and stays symmetric between them by construction.

Not everything can be: inner₊ takes backward differences alone, because the staggered weights it carries are the ones the summation-by-parts identity pairs with a backward difference. Use this alias where the distinction genuinely does not arise.

ExtendedDifferenceNode{D, Dim}

The three difference nodes that are neither one-sided nor a jump, differencing along Dim. Grouped so that everything reading only the direction off a node covers all of them at once, as DifferenceNode does for the one-sided pair.

BackwardDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a backward finite difference operator acting in dimension Dim.

CenteredDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the centered difference along Dim,

\[Dc(u)_i = \frac{u_{i+1} - u_{i-1}}{h_i + h_{i+1}}\]

Truncated at both ends of Dim, having no neighbour on one side.

CrossWeightedDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the cross-weighted centered difference along Dim,

\[D_h(u)_i = \frac{h_i}{h_i + h_{i+1}} D_{-}(u)_{i+1} + \frac{h_{i+1}}{h_i + h_{i+1}} D_{-}(u)_i\]

The same two one-sided differences the centered difference combines, weighted by the opposite spacings. That swap is what makes it second order on a non-uniform grid where Dc is first, and the two coincide when the spacing is constant. Truncated at both ends.

ForwardDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a forward finite difference operator acting in dimension Dim.

StarDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the starred forward difference along Dim,

\[D^{*}_{+}(u)_i = \frac{u_{i+1} - u_i}{(h_i + h_{i+1})/2}\]

The forward difference over the averaged spacing rather than the forward one, which is what makes the discrete integration by parts close. Truncated at the far end of Dim.

grad_backward(op::LazyOp{D}) where D

Constructs a backward gradient operator tuple, yielding D-tuple of BackwardDifference operators.

grad_forward(op::LazyOp{D}) where D

Constructs a forward gradient operator tuple, yielding D-tuple of ForwardDifference operators.

AbstractInnerProduct

Abstract base type for inner product quadrature weights.

BilinearProduct{D,InnerType,LeftType,RightType} <: LazyOp{D}

An AST node representing a bilinear integration term $(u, v)$ in a bilinear form.

InnerH <: AbstractInnerProduct

Quadrature weights for the standard $L^2$ inner product using trapezoidal integration.

InnerPlus{Dim} <: AbstractInnerProduct

Quadrature weights for the modified $L^2_+$ inner product in a specific coordinate dimension Dim.

LinearProduct{D,InnerType,LeftType,RightType} <: LazyOp{D}

An AST node representing a linear integration term $(f, v)$ in a linear form.

inner_plus(left::NTuple{D, LazyOp{D}}, right::NTuple{D, LazyOp{D}}; markers = ()) -> LazyOp{D}

Constructs the sum of directional modified $L^2_+$ inner products across all dimensions.

Each dimension's term is a LinearProduct or a BilinearProduct independently, following _is_source_only on left[dim] exactly as innerₕ does: a gradient tuple of interpolated sources (πₕ(u1), πₕ(u2)) is source-only dimension by dimension.

InterpolationNode{D, S, OpType} <: LazyOp{D}

The symbolic interpolation operator: inner_op lives on src_space, and this node evaluates it at points of whatever mesh the assembly is walking.

Distinct from the source wrapper πₕ(uₕ), which carries a grid function's values. This node carries no values; it carries the map, and its stencil names trial columns.

JumpNode{D, Dim, OpType <: LazyOp{D}} <: LazyOp{D}

AST node representing the jump across interfaces along dimension Dim, $u_{i+1} - u_i$.

Not truncated at the far end: the absent u_{i+1} is taken as zero, yielding -uᵢ there. This matches the space-layer matrix convention, whose boundary row preserves -1.

RegionRestriction{D, RegionType, OpType <: LazyOp{D}} <: LazyOp{D}

AST node representing a spatial restriction of an operator to a specific mesh region or boundary.

Arguments

• region::RegionType: Identifier for the region (e.g. :interior, :boundary, :left, :right, :top, :bottom).
• inner_op::OpType: Underlying operator being restricted.

restrict_to(region, op::LazyOp{D}) -> RegionRestriction

Restrict the operator op to a specific mesh region or boundary identifier.

Examples

# Restrict the trial function to the interior
restrict_to(:interior, U)

# Restrict to a boundary region
restrict_to(:left, U)

ConstraintMarkers

Union representing either unevaluated Dirichlet constraints (DomainMarkers) or time-evaluated constraints (EvaluatedDomainMarkers).

_dirichlet_bc_indices!(A::AbstractMatrix, index_in_marker::BitVector)

Internal helper to apply Dirichlet boundary conditions to matrix A at the indices marked in index_in_marker: each marked row is zeroed and its diagonal set to one.

Costs the boundary cardinality, not ndofs — the marked indices are walked with _each_marked rather than scanned for.

_dirichlet_bc_indices!(A::SparseMatrixCSC, index_in_marker::BitVector)

Apply Dirichlet boundary conditions to a sparse matrix A by directly manipulating its CSC data structure.

A single sweep of the stored values does both halves of the job: entries in a constrained row are zeroed, and the diagonal of such a row is set to one where the sweep meets it, rather than by a second pass afterwards. Explicit zeros are left in place, so the sparsity pattern is unchanged and the matrix can be refilled without reallocating its columns.

_normalize_dirichlet(dirichlet) -> (labels, conditions)

Internal helper turning every form accepted by the dirichlet keyword into the (labels, conditions) pair apply_dirichlet_labels!/apply_dirichlet_conditions! already take – unifying the keyword this way adds no second implementation of applying constraints to keep in step with the first.

Accepts, in order: nothing; a single label Symbol (bilinear only – no values to carry); a Tuple of label Symbols; a single label => f Pair; a Tuple of such Pairs; or constraints already built by dirichlet_constraints, whose own labels are read back out. Anything else throws.

LinearForm{D, TestSpace, AST}

Represents a linear form defined over a test space.

Arguments

• test_space::TestSpace: Space for the test function.
• ast::AST: Resolved expression tree.

The form resolves its expression tree ast once at construction, referencing the underlying storage of any coefficient grid functions (VectorElement). In-place updates via Rₕ!(fₕ, ...) or parent(fₕ) .= ... are automatically seen by subsequent assemblies with zero heap allocations. The expression itself is not retained: downstream routines evaluate the resolved AST directly.

Constant scalar coefficients can be written directly as numbers (e.g. 2.5 * innerₕ(fₕ, v)). Ref is only needed if a dynamic scalar coefficient changes across loop iterations:

α = Ref(1.0)
l = form(Wₕ, v -> α * innerₕ(fₕ, v))
# Inside time loop:
α[] = 2.5
assemble!(b, l) # zero allocations, evaluates with α = 2.5

_colour_strides(offsets) -> NTuple{D, Int}

Per-dimension stride separating grid points that a parallel assembly may write concurrently, for an operator reaching offsets.

Two points of one colour differ by a multiple of the stride in some dimension (at least span + 1 there, where each writes a footprint span wide about itself). Beyond one width apart, the footprints do not overlap: no two points in a colour ever target the same row, enabling lock-free parallel assembly.

An operator reaching only its own point (such as innerₕ(fₕ, v) or any form whose test argument carries no difference) strides by 1 in every dimension, resulting in a single colour.

_scatter_linear_point!(b, sp, term, I, lin_indices, mesh_markers, offset) -> Nothing

Add one grid point's stencil contributions to b.

Shared by the banded and the point-coloured sweep, so the two cannot drift apart.

_sweep_linear_band_colour!(b, sp, term, ax, parity, nbands, rest, lin_indices, mesh_markers, offset) -> Nothing

Scatter one band colour of term into b across threads.

Two points collide only when they write the same entry of b, which needs their difference to lie inside the stencil's reach in every axis at once. Being at least strides[D] apart along the banded axis rules that out on its own, so alternate slabs never race. A term that reaches only its own point cannot collide at all, and then bidx is every band at once.

resolve_form_ast(form::LinearForm)

Return the resolved AST stored inside the linear form.

test_space(form::LinearForm)

Return the test space of the linear form.

AnySegment{D} = Union{NzvalSegment,DiagonalSegment{D}}

Either recorded shape a (term, block) can cache, for a D-dimensional form: a flat NzvalSegment or, where the structure held, a DiagonalSegment. Kept as a two-concrete-member union per D – see DiagonalSegment – rather than leaving D to vary, so Vector{AnySegment{D}} stores unboxed.

BilinearForm{D, TrialSpace, TestSpace, AST}

Represents a bilinear form defined over a trial space and test space.

Arguments

• trial_space::TrialSpace: Space for the trial function.
• test_space::TestSpace: Space for the test function.
• ast::AST: Resolved expression tree.

The form resolves its expression tree ast once at construction, referencing the underlying storage of any coefficient grid functions (VectorElement). In-place updates via Rₕ!(cₕ, ...) or parent(cₕ) .= ... are automatically seen by subsequent assemblies with zero heap allocations. The expression itself is not kept: downstream routines evaluate the resolved AST directly.

Constant scalar coefficients can be written directly as numbers (e.g. 2.0 * innerₕ(D₋ₓ(u), D₋ₓ(v))). Ref is only needed if a dynamic scalar coefficient changes across loop iterations:

β = Ref(1.0)
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(β * D₋ₓ(u), D₋ₓ(v)))
# Inside time loop:
β[] = 3.0
assemble!(A, a) # zero allocations, evaluates with β = 3.0

DiagonalSegment{D}

One term's recorded nzval positions for one block, on a D-dimensional structured grid where the term's own _stencil_margin let its interior peel away from a boundary shell (gpena/Bramble.jl#160): interior entries are base[k] + stride[k] * n for the n-th point interior's own iteration order visits (n zero-based, _interior_rank), rather than one stored Int per entry – positions never carries the interior's O(N * P) share at all. boundary is an ordinary NzvalSegment covering only the shell, indexed exactly as before.

Built by _record_segment! only when every interior point produces the same number of entries P and the same per-tap stride holds across the whole interior – checked once, not assumed, because a form summing terms of different margins can make a column's true nzval footprint vary inside what this one term calls its own interior (see _stencil_margin). Any point where the check fails falls back to a plain NzvalSegment, the general shape ReplaySink already handles.

Parametrized by D alone – interior's ranges-tuple type is pinned to NTuple{D,UnitRange{Int}} (what _interior_range always produces), never left as an independent free parameter – so that for one BilinearForm's fixed dimension, AnySegment{D} is a two-member union of concrete types. A CartesianIndices{D} alone, or a DiagonalSegment with D left to vary, is not concrete (its ranges type is still a UnionAll) and stores boxed: exactly what made an early version of this allocate 80-400 B on every replay, @test_allocs-checked paths included.

See also: DiagonalReplaySink, AnySegment.

NzvalSegment = Tuple{Vector{Int},Vector{Int}}

One term's recorded nzval positions for one block, as (point_ptr, positions).

point_ptr[lin_idx]:point_ptr[lin_idx + 1] - 1 is the slice of positions holding grid point lin_idx's own entries, in the order a scatter walk visits them. Addressed per point rather than by a shared running counter, so a replay stays correct whatever order the grid is visited in.

See also: RecordSink, ReplaySink.

resolve_form_ast(form::BilinearForm)

Return the resolved AST stored inside the bilinear form.

test_space(form::BilinearForm)

Return the test space of the bilinear form.

trial_space(form::BilinearForm)

Return the trial space of the bilinear form.

DiagonalReplaySink(A::SparseMatrixCSC, interior::CartesianIndices, base::Vector{Int}, stride::Vector{Int}, P::Int)

Add a term's values to A's interior core using DiagonalSegment's per-tap stride instead of a stored position per entry.

Paired with an ordinary ReplaySink for the boundary shell through the two-sink form of visit_bilinear_stencil: this sink is only ever handed the interior region, so _sink_point! addresses a point by its rank n in interior's own iteration order (zero-based, via _interior_rank) rather than by lin_idx, matching the order _record_segment! validated the stride against. The k-th tap of that point (1-based, k in 1:P) then lands at base[k] + stride[k] * n, recovered from the running slot the shared walk already threads (slot = n * P + (k - 1)), so no per-entry lookup runs at all.

See also: visit_bilinear_stencil, _replay_segment!.

PatternSink(I_vec::Vector{Int}, J_vec::Vector{Int})

Collect the (row, col) coordinates a term can reach, for building a sparsity pattern.

Appends each coordinate to I_vec and J_vec, which allocate_system_matrix then hands to sparse!. The only sink that de-duplicates (_sink_dedups): a coordinate named twice by one point's stencil is one entry of the pattern, and the weights it carries are not read here at all.

See also: visit_bilinear_stencil, RecordSink.

RecordSink(A::SparseMatrixCSC, term, point_ptr::Vector{Int}, positions::Vector{Int})

Add a term's values to A and record where each entry landed, building the replay cache.

For each entry it searches A for the (row, col)'s slot in nzval, adds the weight there, and appends the slot to positions. _sink_point! opens each grid point's own slice of that list in point_ptr, so a later replay can address a point directly instead of relying on the walk order.

The search is the expensive half of assembly, which is why it is done once and replayed by ReplaySink afterwards.

Throws

• ArgumentError: A (row, col) the pattern does not contain. This pass builds the cache, so a pattern that cannot hold the term is reported rather than skipped, as add_to_sparse! now does on the threaded path too.

See also: visit_bilinear_stencil, NzvalSegment.

ReplaySink(A::SparseMatrixCSC, point_ptr::Vector{Int}, positions::Vector{Int})

Add a term's values to A using slots recorded earlier by RecordSink.

The same walk and the same fresh stencil evaluation, because weights may be live: a coefficient grid function updated in place through Rₕ! is seen by the next assembly. Only the slot lookup is skipped, taken from positions rather than searched for, which is what the cache buys. row and col are ignored for that reason, and _sink_needs_coordinates tells _step_entry! as much, so it skips computing them at all rather than computing and discarding them.

Immutable: the walk's position advances as a loop-local in visit_bilinear_stencil, handed back through slot, rather than as a field of the sink. Carrying it as a mutable field measured about 20% slower on the cheapest replays.

See also: visit_bilinear_stencil.

_entry_target(lin_indices, I::CartesianIndex, off_u, off_v,
              row_offset::Int, col_offset::Int) -> Tuple{Int,Int}

The matrix position a stencil entry writes to, or (0, 0) when it writes nowhere.

The row comes from the test offset off_v and is dropped when it leaves the grid. The column comes from _trial_column, which answers 0 for a trial offset outside the grid and reads an AbsoluteColumn directly, since an interpolation entry names its source column rather than an offset from I. Both are then shifted into the block by row_offset and col_offset.

(0, 0) is a sentinel rather than nothing so the return type stays concrete on the assembly hot path.

Returns

• Tuple{Int,Int}: The (row, col) to write, or (0, 0) to skip the entry.

_sink_dedups(sink) -> Bool

Whether sink wants repeated (off_u, off_v) pairs within one point's stencil collapsed to a single entry.

The sparsity pattern wants each (row, col) once however many stencil entries name it, and a value sink wants every one of them, because repeated entries accumulate. Getting this backwards gives a wrong matrix rather than an error, so it is answered by dispatch on the sink type: the default is false, and the branch folds away at compile time, leaving the value sweeps with no de-duplication scan at all.

See also: visit_bilinear_stencil, PatternSink.

_sink_entry!(sink, row::Int, col::Int, weight, slot::Int) -> Nothing

Act on one stencil entry landing at (row, col) with coefficient weight.

The one method each sink has to supply. Called by visit_bilinear_stencil only for entries that land inside the matrix, so a sink never has to guard the index itself. slot counts accepted entries from this point's base (_sink_point!) and matters only to ReplaySink; the others ignore it.

_sink_needs_coordinates(sink) -> Bool

Whether sink reads the (row, col) an entry computes, rather than discarding it.

The default is true. ReplaySink answers false: its _sink_entry! ignores row and col entirely, since the nzval position is already recorded in sink.positions. false lets _step_entry! skip the linear-index lookups and offset additions that would only be discarded, while still running the bounds checks that decide whether the entry survives at all – dropped entries must match the record pass exactly, or slot drifts out of step with positions.

See also: _step_entry!, _trial_inbounds.

_sink_point!(sink, lin_idx::Int, I::CartesianIndex) -> Int

Announce grid point lin_idx (at cartesian index I) to sink, and answer the base slot for its entries.

The default answers 0. RecordSink uses the call to open that point's slice of the position list; ReplaySink answers the start of that slice, which the traversal then adds the entry ordinal to. Addressing each point from its own base is what lets a replay stay correct regardless of the order grid points are visited in, without any sink having to carry a mutable cursor. I is passed alongside lin_idx for DiagonalReplaySink, which addresses a point by its rank in LinearIndices(interior) rather than by lin_idx – every other sink ignores it.

See also: visit_bilinear_stencil, _sink_entry!.

_step_entry!(sink, lin_indices, I, off_u, off_v, weight, row_offset::Int, col_offset::Int, slot::Int) -> Bool

One entry: guard it, and hand it to the sink if it lands inside. Answers whether it did, so the caller can advance the slot.

Fused rather than "compute the target, then act on it" (_entry_target, which the tests use and _scatter_point! shares) because returning a (row, col) sentinel tuple has to be merged from three return points: measured against the hand-written loop it cost 5 extra phi nodes, 4 integer adds and 3 comparisons per entry, with identical loads, stores and calls.

_step_entry_unguarded!(sink, lin_indices, I, off_u, off_v, weight, row_offset::Int, col_offset::Int, slot::Int) -> Nothing

_step_entry!, without the guard: the interior-core counterpart, called only where _stencil_margin already guarantees every offset lands inside the grid. Always accepts its entry, so it has nothing to answer back and the caller advances slot unconditionally rather than being told to.

_trial_column(lin_indices, I::CartesianIndex, off_u) -> Int

Which column of the trial block a stencil entry's trial slot names, or 0 for none.

An ordinary offset is bounds-checked against the grid and answers 0 on a boundary, so the entry is dropped. An interpolation entry carries an AbsoluteColumn instead, naming a source column outright, because the trial degrees of freedom it reaches live on a different mesh and which ones depends on where the point falls.

See also: _entry_target.

_trial_column_unguarded(lin_indices, I::CartesianIndex, off_u) -> Int

_trial_column, without the bounds check.

Only called from the interior core of visit_bilinear_stencil, where _stencil_margin has already guaranteed I + CartesianIndex(off_u) lands inside lin_indices for every entry the term's stencil can produce. An AbsoluteColumn names a source column outright either way, matching _trial_column.

_trial_inbounds(lin_indices, I::CartesianIndex, off_u) -> Bool

Whether off_u's trial index survives, without computing the linear index _trial_column would return – for sinks that only need to know whether the entry lands, not where (_sink_needs_coordinates). An AbsoluteColumn always survives, matching _trial_column: it names a source column directly and is never 0.

add_to_sparse!(A::SparseMatrixCSC, row::Int, col::Int, val::Number, term) -> Nothing

Add val to A[row, col], which the preallocated sparsity pattern is required to contain.

Throws

• ArgumentError: (row, col) is not a stored entry of A.

This used to return quietly on a missing entry, which let a matrix whose pattern cannot hold the form assemble to a plausible wrong answer instead of failing. That is how a term naming both components of a composite space once vanished without a word: the pattern held the diagonal blocks only, so every off-diagonal contribution was discarded. The regression test for that is test/form/bilinear.jl, "Composite blocks".

term is carried only to name the offending node in the message, and read only on the branch that throws.

See also: allocate_system_matrix and RecordSink, which raises the same way on the serial recording pass.

visit_bilinear_stencil(sink, term, sp, row_offset::Int, col_offset::Int) -> sink
visit_bilinear_stencil(interior_sink, boundary_sink, term, sp, row_offset::Int, col_offset::Int) -> boundary_sink

Walk every grid point of sp, evaluate term's local stencil there, and hand each entry that lands inside the matrix to a sink.

The walk that five separate sweeps used to re-derive: mesh, markers, linear indices, the loop over grid points, the stencil evaluation, and the decision of where each entry lands (_entry_target). A sink supplies only what to do with an entry, so the pattern passes and the value passes share one traversal instead of agreeing by coincidence.

The fifth argument handed to local_stencil is the leaf's linear index. A term routed to the wrong leaf reads the wrong SourceVector without complaint, which is why the index is derived here rather than by each caller.

Splits into an interior core and a boundary shell when _stencil_margin and the grid size allow it (see the comment above _peelable), so most points skip the bounds guard entirely; every point still gets exactly one visit either way, so a sink sees the same set of entries regardless of which path ran, in a possibly different order. RecordSink/ReplaySink are unaffected by that: they address each point by its own linear index, not by visit order (see their docstrings).

The two-sink form lets the interior and the boundary shell be handled by different sinks – only DiagonalReplaySink needs this, pairing itself (interior) with an ordinary ReplaySink (boundary shell), so the one-sink form below is the thin, common case.

Arguments

• sink (or interior_sink/boundary_sink): What to do per entry. See PatternSink, RecordSink, ReplaySink and DiagonalReplaySink, and the contract in _sink_entry!, _sink_point! and _sink_dedups.
• term: The AST node whose stencil is evaluated at each point.
• sp: The test leaf whose grid is walked and whose markers the stencil sees.
• row_offset, col_offset: The block's origin in the assembled matrix, 0 for a scalar space.

Returns

• The boundary-shell sink (sink itself, in the one-sink form), so a caller can read what it collected.

See also: allocate_system_matrix, add_to_sparse!.

_sweep_band_colour!(A, sp, term, ax, parity, nbands, rest, lin_indices, mesh_markers, row_offset, col_offset) -> Nothing

Scatter one band colour of term into A across threads.

Each thread takes one slab of the last axis and walks it whole. Two grid points can only reach the same matrix entry when they are closer than strides[D] along that axis – their stencil footprints cannot meet otherwise – so slabs of at least that width, taken every other one, never write the same entry concurrently, whatever the remaining axes do. A term that reaches only its own point cannot collide at all, and then bidx is every band at once.