Forms#

_batch_linear_band_sweep!(b, sp, term, ax, bidx, nbands, rest, lin_indices, mesh_markers, offset, α) -> Nothing

CpuPolyester's counterpart of the Threads.@threads body in _sweep_linear_band_colour!, filled by BramblePolyesterExt. Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_batch_linear_colour_sweep!(b, sp, term, idxs, lin_indices, mesh_markers, offset, α) -> Nothing

CpuPolyester's counterpart of the Threads.@threads body in _sweep_colour!, filled by BramblePolyesterExt. Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_zero_stored!(A::AbstractMatrix) -> AbstractMatrix

Zero every value a refill will touch, before _assemble_bilinear_core_cached!/ _assemble_bilinear_parallel_cached! add each term's contribution back in.

Part of the matrix-type seam (gpena/Bramble.jl#12): SparseMatrixCSC only has to zero its stored nzval entries (the sparsity pattern itself never changes between refills), while a dense Matrix – the seam's positive control – has no such distinction and zeros the whole backing array. A future backend (tridiagonal, banded, ...) implements whichever of the two shapes its own storage has.

_scatter_add!(A::AbstractMatrix, pos::Integer, val) -> Nothing

Add val at the position _scatter_position named in A's own storage.

SparseMatrixCSC's method writes nzval[pos]; the generic AbstractMatrix fallback writes A[pos] directly (linear indexing into the backing array), the other half of the seam _scatter_position documents. pos::Integer on every method: a DeviceLocality matrix's own mirror (gpena/Bramble.jl#313) searches in its own index type, which need not be Int.

_scatter_position(A::AbstractMatrix, row::Int, col::Int) -> Int

Where entry (row, col) lives in A's own storage, or 0 for a matrix type that can answer "not stored" (a CSC or CSR search, host or mirror, does; the dense fallback always answers a valid position).

The seam a new backend's matrix type implements to plug into assembly (gpena/Bramble.jl#12): DiagonalReplaySink, the positions search of the coordinate walk (_coordinates_to_positions!) and add_to_sparse! reduce to this and _scatter_add! once the raw nzval/ linear-index field access each used to do directly is factored out here. ReplaySink and its kin replay those positions into _scatter_storage(A), the array _scatter_add! writes.

The SparseMatrixCSC method is exactly the position search this file always ran: a linear scan of the column when it holds few entries, a binary search otherwise – both rely on SparseMatrixCSC's own invariant that rowval is sorted within each column. It returns Int whatever the matrix's index type. The column bounds are converted before the search, so an Int32 matrix answers the same Int positions as an Int one. The generic AbstractMatrix fallback needs no search at all: every (row, col) inside the matrix's bounds is "stored" for a dense backend, at LinearIndices(A)[row, col].

_allocate_from_pattern(::Type{MT}, nrows::Int, ncols::Int, I::Vector{Int}, J::Vector{Int}, V::AbstractVector) -> MT

Build the nrows × ncols matrix a form's sparsity pattern describes, from the coordinate triplet (I, J, V) the coordinate walk (_form_coordinates) collected – summing V[k] into any (row, col) that I/J name more than once, matching sparse!'s own combiner.

The one place a fresh system matrix is born (gpena/Bramble.jl#12): every backend's matrix type implements exactly this to be usable with allocate_system_matrix. SparseMatrixCSC's method is sparse! itself, consuming I/J in place. The generic AbstractMatrix fallback – the dense Matrix{Float64} positive control among them – allocates zeros and scatters into it, since a dense matrix has no sparsity pattern to build.

_batch_bilinear_band_replay!(target, sp, term, ax, bidx, nbands, rest, lin_indices, mesh_markers, row_offset, col_offset) -> Nothing

CpuPolyester's counterpart of the Threads.@threads body in _sweep_band_colour! when it replays: for each band b in bidx, every I in CartesianIndices((rest..., _band_range(ax, nbands, b))) gets _replay_point!(target, term, sp, I, lin_indices, mesh_markers, row_offset, col_offset). target is a _ReplayTarget or an action target, crossing @batch as _batch_bilinear_colour_replay! describes. Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_batch_bilinear_colour_replay!(target, sp, term, idxs, lin_indices, mesh_markers, row_offset, col_offset) -> Nothing

CpuPolyester's counterpart of the Threads.@threads body in _sweep_bilinear_colour! when it replays: the same loop over idxs, calling _replay_point!(target, term, sp, I, lin_indices, mesh_markers, row_offset, col_offset) per point. target is a _ReplayTarget or an action target. The loop captures none of them whole: the target's storage vector (nzval, or y and x) crosses as a top-level argument, the rest of target, sp, term and mesh_markers as one _batch_split, rebuilt in each task. Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

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}

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.

ZeroOperator(Wₕ::AbstractSpaceType)

The zero operator on Wₕ, as a symbolic node. Absorbs multiplication by a scalar. Keeps only the dimension of Wₕ: two zero operators of one dimension are the same node.

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

AbsoluteRow

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

The mirror of AbsoluteColumn, and it exists for the mirror reason: a test-side interpolation (innerₕ(u, πₕ(w))) determines which rows a point's contribution scatters into, on a mesh other than the one being walked, so it names them outright. The walked mesh is then the trial function's – whichever side stays native is the side that supplies the quadrature weight and the grid being swept.

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 function wants a point, and resolve_ast throws an ArgumentError naming both fixes. 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 stencil evaluated at the neighbour 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.

_grid_function_value(grid_function, lin_idx::Int)

The value a GridFunctionScale coefficient contributes at linear index lin_idx.

grid_function is whatever GridFunctionScale.grid_function may hold: a zero-argument Function thunk deferred until resolve_ast calls it (see SourceVector's docstring for why a thunk, not a plain vector, is the wire format), a bare Number for a uniform scale, or an array of per-point values. Shared by local_stencil(::GridFunctionScale, ...) (assembly/stencil_eval.jl), which reads the coefficient at the point being visited, and the shifted_inner_stencil override below, which reads it at the point a tap shifts to – previously the same three-way branch, written out twice.

_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. For a source the Bool itself scales the stencil read at the clamped point, so the stencil keeps one length and one type at every point.

entry_offsets(stencil::Tuple)
entry_weights(stencil::Tuple)

A stencil's offsets and its coefficients, as two separate containers: entry_offsets keeps each entry's offsets alone ((off_u, off_v) for a bilinear entry, (off_v,) for a linear one) and entry_weights keeps the coefficients.

Taken apart for Enzyme, which cannot type a stencil entry's mixed Int/Float64 tuple in a function that reads both halves of it and is not inlined into the function being differentiated. All three conditions are needed, measured one at a time: the same loop written inside the differentiated closure compiles, so does one that reads only the offsets, and so does one that reads only the weights – what fails is the combination, with EnzymeNoTypeError inside _visit_guarded_region!. Reading the offsets from a container of Ints and the weights from a container of Float64s removes it, for a coefficient scaling the form and for a VectorElement coefficient alike, in 1D, 2D and 3D.

The stencil itself is unchanged: local_stencil returns what it always did, and the split happens where the entries are consumed (_visit_entries). That is enough, and the narrower change: every local_stencil method, the stencil algebra above and the tests that compare stencils against literal tuples all stay as they are. Enzyme differentiates local_stencil and scale_stencil themselves correctly at any stencil size measured, up to 63 machine words – a failure at a large stencil is really _peelable selecting the guarded walk, not an aggregate Enzyme cannot type.

Both are an ntuple over Val(length(stencil)), which unrolls at any length, so neither call allocates whatever the stencil's size. map would not do: over a Tuple of more than 32 entries Base's map collects into a Vector{Any} and boxes, and a product of two nested operands (innerₕ(D̽ₓ(D₋ₓ(u)), D̽ₓ(D₋ₓ(v)))) is that long at every grid point.

entry_offsets(stencil::Tuple)
entry_weights(stencil::Tuple)

A stencil's offsets and its coefficients, as two separate containers: entry_offsets keeps each entry's offsets alone ((off_u, off_v) for a bilinear entry, (off_v,) for a linear one) and entry_weights keeps the coefficients.

Taken apart for Enzyme, which cannot type a stencil entry's mixed Int/Float64 tuple in a function that reads both halves of it and is not inlined into the function being differentiated. All three conditions are needed, measured one at a time: the same loop written inside the differentiated closure compiles, so does one that reads only the offsets, and so does one that reads only the weights – what fails is the combination, with EnzymeNoTypeError inside _visit_guarded_region!. Reading the offsets from a container of Ints and the weights from a container of Float64s removes it, for a coefficient scaling the form and for a VectorElement coefficient alike, in 1D, 2D and 3D.

The stencil itself is unchanged: local_stencil returns what it always did, and the split happens where the entries are consumed (_visit_entries). That is enough, and the narrower change: every local_stencil method, the stencil algebra above and the tests that compare stencils against literal tuples all stay as they are. Enzyme differentiates local_stencil and scale_stencil themselves correctly at any stencil size measured, up to 63 machine words – a failure at a large stencil is really _peelable selecting the guarded walk, not an aggregate Enzyme cannot type.

Both are an ntuple over Val(length(stencil)), which unrolls at any length, so neither call allocates whatever the stencil's size. map would not do: over a Tuple of more than 32 entries Base's map collects into a Vector{Any} and boxes, and a product of two nested operands (innerₕ(D̽ₓ(D₋ₓ(u)), D̽ₓ(D₋ₓ(v)))) is that long at every grid point.

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.

_each_entry unrolls and stays type-stable at any stencil length, with zero allocations, so code generation would buy 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. Called from a ShiftNode, inner is that stencil only for a bare TrialFunction or TestFunction, and nothing for any other operand, whose paths below re-evaluate it at the shifted point and never read inner.

The one place the "shift by relabelling" assumption is made, so the one place a node that cannot be relabelled has to be handled. The default dispatches on stencil_shift_trait: TranslationInvariantStencil evaluates inner_op at the shifted point and relabels its offsets (a bare trial or test leaf relabels inner itself); PointDependentStencil evaluates inner_op at the shifted point, relabelling its offsets only when they are offsets (_has_relative_offsets): not for an interpolation's absolute columns or a source's value.

Five node types override this default outright rather than answering through the trait alone. Two are below; the others are RegionRestriction and OperatorScale (operators/regionrestriction.jl) and TappedNode (assembly/stencileval.jl). The trait's two stock branches cannot express what the two below need. Re-evaluating the whole subtree at the shifted point is exact for a source, which carries only a value, but wrong for anything that also carries a trial or test column, since that column has to move by relabelling while whatever multiplies it has to be read at the new point:

• GridFunctionScale shifts its operand by the operand's own rule (a recursive call to this same function), then reads its coefficient at the shifted point and scales by it, so D₋ₓ(cₕ * u) moves the trial column and re-reads cₕ at the tap.

• OperatorAdd recurses into each summand and concatenate_stencils the two results, since a sum can mix a point-dependent summand (one carrying a GridFunctionScale) with a translation-invariant one, and the generic point-dependent branch would re-evaluate the translation-invariant summand at the shifted point too, discarding its trial column.

Every branch – the two default ones and the overrides – produces 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 (operators/inner.jl) for a source-only subtree's own local_stencil. The offsets 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, since 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. That is exactly sum(f, itr; init = false)'s own behavior, which is what this calls. A @generated unrolled fold would buy nothing: sum already gives zero allocations and the same inferred type.

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, extended and second differences, the three 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.

Fourteen 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 fourteen 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.

The stencil nodes are a member as StencilNode{D}, not one by one: a method on StencilNode (_widened, bilinear.jl) must stay more specific than one on this union, and a list of the eleven nodes outranks the abstract type (gpena/Bramble.jl#930).

_face_margins(term, ax::NTuple{D, AbstractUnitRange{Int}}) -> NTuple{D, NTuple{2, Int}}

The walk's guarded margin on each face of the grid with axes ax, (low, high) per axis: _stencil_margin(term) on every face. A marked function of term sits only inside a mixed difference (_mark), so it reaches at most m - 1 along any one axis, m = _stencil_margin(term), and an interior point m in from a marked face reads the first node in from that face, never its ghost layer: the interior takes the term unmarked (_interior_term) without widening the margin on a marked facet.

_is_source_only(op::LazyOp) -> Bool

Whether a LazyOp subtree is source-only: built entirely from sources (SourceFunction, SourceVector, SourceConstant, DiracSource) 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, since a BilinearProduct's stencil carries a pair of offsets (trial and test), where _scatter_term! (assembly/linear.jl) expects one.

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

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.

The index is a value, not part of the type. A literal index such as v(1) + v(2) still gives form one concrete return type. A runtime component index (v(i) with i a function argument) makes the form's type a Union when a sum of indexed terms sits inside one inner product (innerₕ(u(i), v(i) + D₋ₓ(v(i)))), or under a scalar or grid-function coefficient shared by the whole sum (3 * (innerₕ(u(i), v(i)) + innerₕ(u(i), D₋ₓ(v(i))))), even when every term names the same i: whether the sum is split into one term per component is then decided at run time. Separate inner products, each carrying its own coefficient (3 * innerₕ(u(i), v(i)) + 3 * innerₕ(u(i), D₋ₓ(v(i)))), infer concretely. The form assembles the same either way; where a concrete type matters (automatic differentiation, a function barrier), use literal components or separate inner products.

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). row_offset always comes from test_leaf's offset and col_offset from trial_leaf's, so a caller reads them off the type instead of unpacking a positional (tc, sc) tuple.

_NamedBlock{TrialLeaves, TestLeaves}

The block a term names when its leaves differ in type, a Float64 leaf beside a Float32 one say. It holds the leaf tuples and the runtime (trial, test) components, not yet resolved to a Block.

Why not resolve it at once? Indexing a heterogeneous tuple with a runtime component gives a Union of leaf types, and a Block built from that is abstract, so every consumer would dispatch on it per block and box what it computes. _foldl_blocks walks the tuples statically instead. Homogeneous leaves never come here: a runtime index into Tuple{Vararg{T}} is concrete and compiles one block walk, not one per leaf pair.

_at_leaf(f, leaves, c) -> f(leaves[c])

f applied to leaf c of a leaf_spaces_offsets tuple, c already range-checked. A runtime index on homogeneous leaves; a static walk on leaves of different types, so f is called on a concrete leaf. The walk's recursion is @inline, so a literal component (v(2)) reaches the comparisons as a constant and folds the walk to one leaf: without it, constant propagation stopped at the recursion and form inferred a Union of the per-leaf results on three leaves or more.

_bind_marker_ids(op, Ωₕ::AbstractMeshType) -> LazyOp

op rebuilt with each RegionRestriction region replaced by the column of its label in Ωₕ's word matrix (_marker_id): a Symbol becomes an Int, a tuple of them an NTuple{N, Int}. Recurses exactly where _collect_region_labels does, so every label validated is bound. Labels are read as stored, no boundary alias resolved.

_bind_walk(term, Ωₕ::AbstractMeshType) -> Tuple{bound_term, Union{AbstractMatrix{UInt64}, Nothing}}
_bind_walk(term, sp::ScalarGridSpace) -> Tuple{bound_term, Union{AbstractMatrix{UInt64}, Nothing}}

term with every RegionRestriction it holds bound to Ωₕ's marker ids (_bind_marker_ids), and the marker table the walk over Ωₕ passes to local_stencil: Ωₕ's word matrix, any AbstractMatrix{UInt64} (_marker_words), or nothing when term restricts nothing, decided from its type.

Called at every walk entry, once per walk, so the ids are read from Ωₕ's label table as it is when the walk runs: assemble on a form built before a markers! call binds against the new table, and a label removed since throws the ArgumentError of _validate_term_markers. Editing markers(Ωₕ) in place is unsupported, and markers! changes the labels. A label written in place has no id and throws too (_throw_marker_not_bound). Given the walked leaf sp instead, it binds against mesh(sp) after checking sp's weights (weights), since its stencils read the weights unchecked.

_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.

_foldl_block_pairs(f, acc, b1, b2) -> acc

acc = f(acc, blk, blk2) for the k-th Blocks of b1 and b2, two equal-length blocks results: _foldl_blocks for a transposed pair, whose consumers need both blocks at once. A _NamedBlock on either side is resolved by nesting the two walks.

_foldl_blocks(f, acc, bs) -> acc

acc = f(acc, blk) for every Block of bs, a blocks result, in order. A homogeneous tuple of Blocks is a plain loop; any other tuple (one Block type per leaf pair, or a _NamedBlock) is walked by tail recursion, so f sees a concrete Block at every call. The accumulator is threaded by return value, never captured, so a counter such as the replay's next is not boxed.

_foldl_leaves(f, acc, leaves) -> acc

acc = f(acc, c, leaf) for each (c, leaf) of a leaf_spaces_offsets tuple. Homogeneous leaves loop with a runtime index; leaves of different types are walked by tail recursion, so f sees each leaf's own type rather than their Union.

_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/_coord_walk!), 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

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. A named block over leaves of different types is a _NamedBlock, so the result is concrete either way; walk it with _foldl_blocks to reach each Block concretely.

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 (assembly/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.

One exception: a marked function (_FacetExtrapolated) answers its function's margin, as its reads by extrapolation go inward from a ghost to the three nodes _marked_ghost requires; the walk's interior reads no marked ghost without widening (_face_margins).

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. A RegionRestriction is the exception: outside its region it contributes no entry at all, so the set is the reach wherever it contributes. 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.

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.

inner is the weight, an instance of InnerType; the walks read it, and InnerType stays the type parameter dispatch reads.

InnerGamma{M} <: AbstractInnerProduct

Quadrature weights for the \((D-1)\)-dimensional surface integral over the grid faces mask names (gpena/Bramble.jl#157).

Building the node holds the labels, mask::_SurfaceLabels{L}; form resolves them once, against the mesh of the test space (_resolve_surfaces), into the per-axis (min, max) face mask _surface_mask (mesh/queries.jl) answers, mask::NTuple{D, NTuple{2, Bool}}. The mask is a field, not a type parameter, so the resolved type is fixed by D and form infers whatever faces the labels name (gpena/Bramble.jl#895).

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.

InnerPlusSet{S} <: AbstractInnerProduct

Quadrature weights for the modified \(L^2_+\) inner product staggered in every direction the set S names at once, for |S| \geq 2 (gpena/Bramble.jl#115, #234).

The empty set is InnerH and a one-element set is InnerPlus: those two keep their own node types rather than becoming a special case of this one, because other code matches on their literal types (src/assembly/kronecker.jl's separability match, and the typeof(inner₊ₓ(id, id)).parameters[2] === InnerPlus{1}-style pin in test/form/inner_products.jl) and widening what they resolve to would change what those match. inner₊(u, v, Val(S)) is what builds this node; it is never constructed for S shorter than 2.

S is stored sorted (_canonical_set, this file) so two callers naming the same set in a different order – Val((1,2)) and Val((2,1)) – produce the identical singleton type, which is what lets which(inner₊, ...) resolve one way regardless of the order the caller wrote S in.

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

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

inner is the weight, an instance of InnerType, as for BilinearProduct.

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.

RegionRestriction{D, RegionType, OpType <: LazyOp{D}, Z} <: 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), or a tuple of them. Assembly binds it to the walked mesh's marker id (Int, or NTuple{N, Int} for a tuple; see _bind_marker_ids).

• inner_op::OpType: Underlying operator being restricted.

• zero_stencil::Z: what the node answers outside its region: nothing for the empty stencil (), or, for a restricted source a form has probed (_typed_restrictions, assembly/linear.jl), a zero stencil of the type the source has inside its region.

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

Restrict the operator op to a mesh region named by a marker label.

region is one of two forms:

• a Symbol, the label of a domain marker (:interior, :boundary, :left, a name the domain was given);

• a tuple of Symbols, the union of those labels.

Examples

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

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

# Restrict to the union of two regions
restrict_to((:left, :right), U)

Throws

• ArgumentError: if region is neither a Symbol nor a tuple of Symbols. A region given by a predicate is not supported: name it as a domain marker instead (domain(S, :name => x -> ...)) and restrict to :name.

_node_alias_expr(alias_name, dispatch_name, direction, direction_index, what; source = nothing)

Returns the expression defining one per-coordinate alias, alias_name(op) as dispatch_name(op, Val(direction_index)), with its docstring attached.

_node_bind_expr(node_name; source = nothing)

Returns the expression defining _bind_interp_spaces for one node family.

The interpolation binding pass (operators/interpolation.jl) walks a term the way resolve_ast does, and a directional node's share of that walk is the same line per family as its resolve_ast: rebuild with the bound operand inside.

_node_dimensional_expr(node_name, alias_name, alias_stem, what; source = nothing)

Returns the expressions defining a node family's dimensional entry point, alias_name(op::LazyOp{D}, ::Val{Dim}), as a Vector{Expr}.

Val only, and deliberately: Dim lands in the node's type, so an Int direction would make the return type a Union of the three nodes it could be. The space layer's Int and Symbol forms have no counterpart here.

_node_resolve_expr(node_name; source = nothing)

Returns the expression defining resolve_ast for one node family: resolve the inner operator once and rebuild the node around it.

_node_vectorial_expr(alias_name, dispatch_name, what; componentwise = false,
                     source = nothing)

Returns the expressions defining a node family's tuple-valued alias, as a Vector{Expr}: alias_name(op) giving a D-tuple of nodes, and the bare node rather than a one-tuple in one dimension, as ∇ₕ and ∇₊ₕ have always done.

componentwise = true adds the method that maps the alias over a tuple of scalar symbolic functions, which is how a composite space's components reach it. Only the two gradients carry it; nothing calls Mₕ((u1, u2)).

ntuple with a Val length unrolls, so dim is a literal in each unrolled body and the Val(dim) it builds is a compile-time constant – not the boxed Val of gpena/Bramble.jl#146, which came from a Val(i) under a runtime loop bound.

@node_family(kwargs...)

Defines one AST node family's whole alias surface and its resolve_ast, declaratively.

Keywords, all optional except node, stem and what:

_bind_interp_spaces is generated unconditionally: it is the same rebuild as resolve_ast, and a family that opts out of the latter still has to pass a bound operand through.

The space-layer counterpart is @operator_family in operators/stencil.jl, and the generators it shares – _relocate! for line attribution, and _BRAMBLE_var2symbol / _BRAMBLE_var2label for the coordinate suffixes – are the same ones used here.

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. A marked row with no diagonal column (A has fewer columns than rows) is zeroed but gets no identity entry, as in the sparse method.

Costs the boundary cardinality, not ndofs: the marked indices are walked with _each_marked rather than scanned for. A matrix in device memory takes _dirichlet_bc_device! instead.

_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.

index_in_marker has one entry per row, and A may have more columns than rows (a scalar test space against a composite trial space), so each lookup is bounds-checked through _row_marked, as in _dirichlet_bc_rows!. A marked row with no diagonal column is zeroed but gets no identity entry.

_has_device_csr_mirror(::Type{<:AbstractMatrix}) -> Bool

Whether a device matrix type stores CSR arrays rowPtr/colVal/nzVal plus a host mirror (_DeviceSparseMirror), so Dirichlet rows can be rewritten by a device kernel. false unless an extension (BrambleMetalExt for MetalSparseMatrixCSR) says so.

_launch_dirichlet_rows_csr!(rowPtr, colVal, nzVal, rows::Vector) -> Nothing

Zero each row r in rows of the device CSR matrix stored as rowPtr, colVal, nzVal and set its stored diagonal to one, one work item per row. rows is a host vector, uploaded once; the call blocks until the kernel has run, since the upload reads from rows.

Requires using KernelAbstractions; the real method is supplied by BrambleKernelAbstractionsExt.

Throws

• ErrorException: if KernelAbstractions is not loaded.

_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, α = true) -> 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!(policy, b, sp, term, ax, bidx, nbands, rest, lin_indices, mesh_markers, offset, α = true) -> 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.

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:

using Bramble: D₋ₓ
β = 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

(a::BilinearForm)(u, v) -> Number

The form evaluated at a trial function u and a test function v: vᴴ A u for A = assemble(a), without A. One walk of the form's stencil, the one matrix_free_operator's serial mul! makes, adds conj(v[row]) * weight * u[col] for every entry into a single scalar. No Dirichlet rows are applied, as assemble(a) applies none.

u and v are VectorElements or plain vectors, of lengths ndofs(trial_space(a)) and ndofs(test_space(a)). The result's type is promoted from the element types of A, u and v, so a ForwardDiff.Dual in either argument comes through.

On a CPU backend the call allocates only the one cell its sum runs in (a ContractionSink's), whatever the grid size. A product with three or more region-restricted summands is split by form into products of at most two (gpena/Bramble.jl#727); one it leaves whole, as it would split into more than 16 or its constant scales keep several types, allocates per point, as in assemble!. Each call has its own cell, so tasks may call one form at once. On a GPU backend the matrix is assembled and multiplied, as the walk reads its vectors on the host.

Throws

• DimensionMismatch: u or v has the wrong length.

• ArgumentError: u or v is not indexed from 1 (an offset-axes view, say).

Examples

using Bramble, LinearAlgebra
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
uₕ = Rₕ(Wₕ, x -> sin(x[1]))
vₕ = Rₕ(Wₕ, x -> x[1] + 1)
a(uₕ, vₕ) ≈ dot(parent(vₕ), assemble(a) * parent(uₕ))

# output
true

Segment{D}

One term's recorded nzval positions for one block, on a D-dimensional structured grid. A single concrete struct, with is_diagonal discriminating the two shapes (gpena/Bramble.jl#240): Enzyme's strict-aliasing type analysis rejects a Union return outright (IllegalTypeAnalysisException), regardless of whether the segment itself is ever the target of differentiation, so _try_diagonal_segment returns a type Enzyme can reason about, never a two-branch Union.

is_diagonal selects which shape is populated:

• false (flat). point_ptr/positions cover every entry; base/stride/interior are unused placeholders (empty, zero).

• true (diagonal, gpena/Bramble.jl#160). point_ptr/positions cover only the boundary shell; base/stride/P/interior carry the interior's per-tap stride arithmetic – interior entries are base[k] + stride[k] * n for the n-th point interior's own iteration order visits (n zero-based), rather than one stored Int per entry, so positions never carries the interior's O(N * P) share at all. 1D forms only (_diagonal_replay), since from 2D up no difference term's interior has a constant stride.

Built by _try_diagonal_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 the flat shape instead, the general case 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, Segment{D} is concrete. A CartesianIndices{D} alone, or a Segment 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.

positions_t is empty except on a segment recorded for a transposed pair ⟨Au, Bv⟩ + ⟨Bu, Av⟩ (_segments_from_positions): there it holds, entry for entry, the nzval position of the transposed entry the second term writes, so one walk of ⟨Au, Bv⟩ fills both. The half a pair's unit does not write (_PairReplaySink) leaves positions or positions_t empty. A pair segment is always flat.

See also: DiagonalReplaySink.

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::AbstractMatrix, interior::CartesianIndices, base::Vector{Int}, stride::Vector{Int}, P::Int, α)

Add a term's values to A's interior core using Segment'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 an incrementing counter on the sink itself) rather than by lin_idx, matching the order _try_diagonal_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!.

ReplaySink(A::AbstractMatrix, point_ptr::AbstractVector{Int}, positions::AbstractVector{Int}, α)

Add a term's values to A using slots the coordinate walk recorded (_CoordSink, searched in A by _coordinates_to_positions!).

Holds A's storage (_scatter_storage: the nzval vector of a sparse matrix) and its positions as any AbstractVector{Int}, so a sparse sink is plain vectors throughout (gpena/Bramble.jl#437).

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.

_CoordSink(point_ptr, I, J, base, tbase, dr, dc, half, fill, n)

The coordinate walk's sink (_coord_walk!, bilinearpattern.jl): opens each grid point's slice of the unit's entries in `pointptr(counting pass) or writes each entry's(row, col)intoI[base + n]/J[base + n](filling pass,fill = true),ncounting the unit's entries. For a transposed pair's unit (half >= 0) the filling pass also writes the transposed coordinate(col + dr, row + dc)attbase + n;half = 1writes the direct coordinate only andhalf = 2the transposed only, as_PairReplaySink` scatters them.

One concrete type for both passes and every unit, so a term's walk compiles once. Entries are written with setindex! into vectors sized by the counting pass, never grown with push!: push! in a recording pass is what made Enzyme unable to compile it when a term's own coefficient is differentiated (gpena/Bramble.jl#240).

See also: visit_bilinear_stencil, ReplaySink.

_DeviceSparseMirror{Tv, Ti}(rowptr::Vector{Ti}, colval::Vector{Ti}, nzval::Vector{Tv})

The host staging area behind a device sparse matrix's own mirror field.

rowptr/colval are host copies of the matrix's CSR structure, carried in the matrix's own index type Ti rather than widened to Int: on a device where the index and value arrays share physical memory with the host (Apple Silicon's unified DRAM), a mirror hardcoded to Int costs more host memory than the device arrays it mirrors the moment Ti is narrower than Int (Int32, say – half the footprint per index). nzval is what every scatter into that matrix actually accumulates into – the device array itself is left untouched until _flush_device_scatter! copies nzval across in one bulk copyto!. Any device sparse matrix type that wants to plug into assembly carries exactly this shape in a field named mirror: _scatter_position/_scatter_add! below duck-type it off A.mirror directly, so no file under src/ needs to name the concrete matrix type (BrambleMetalExt.MetalSparseMatrixCSR today). See docs/notes/internals/gpu.md.

_PairReplaySink(A, point_ptr, positions, positions_t, α1, α2, half::Int)

ReplaySink for a transposed pair ⟨Au, Bv⟩ + ⟨Bu, Av⟩, walking ⟨Au, Bv⟩ only. An entry adds α1 * weight at its recorded position and α2 * weight at its recorded transposed one, the entry (col + dr, row + dc) the second term writes (_CoordSink): at a grid point the second term's stencil is the first's with trial and test offsets exchanged.

half is 0 for both entries, 1 for the first term's only and 2 for the transposed only: a pair whose two blocks sit on different leaf objects walks ⟨Au, Bv⟩ once per leaf, one half each, since the second term's stencil is the first's exchanged only on its own leaf. The skipped half's positions are never read.

_WalkRows

The part of the padded grid a form on ghost spaces walks: the points points whose stencil writes a row of rows, the padded rows its test closure reads, and the columns cols those points read (_ghost_walk, ghostassembly.jl). whole when the closure reads every row, a ghost box: the walk is then the plain one. salt tells a recording of this walk from one of the whole grid on one matrix (`walkstamp).startsmoves the fused matrix-free product's bands so that each holds as many walked points (walkown), andbands[span]holds the starts of the replay's balanced bands for a unit of row reachspan(cut_slices`); both are empty on a whole walk.

Mutable, so that a threaded sweep's task captures it as one pointer: inline, its three _GridLines grew each Threads.@threads closure of a CpuThreaded refill by 144 bytes. A CpuPolyester task, which takes no mutable value, receives its _WalkCut instead.

_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. This is the reference form the tests pin; it is not on the hot path, where _step_entry! fuses it.

Returns

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

_flush_device_scatter!(A::AbstractMatrix) -> Nothing

Copy A's own host-staged mirror (A.mirror.nzval) back to its device storage in one bulk copyto!, and block until that copy has actually landed.

A no-op for a host matrix. bilinear_execution.jl calls this exactly once, after its band-coloured sweep returns.

One copyto! of the whole array per assembly, not one per entry (the first version scattered each entry straight to A.nzVal piecemeal; correct, but hundreds of thousands of 4-byte device transfers on a real matrix).

The copyto! alone was not enough either: a device write queues asynchronously exactly like a @kernel launch does, and every kernel launch in this milestone is followed by a synchronize – this one was not, so assemble/assemble! could return before the copy had actually landed, and a caller reading A right after got whatever was there so far. ka_synchronize is the fix, following that same kernel-launch idiom for a caller that has a device array but no kernel of its own.

_replay_add!(nzval::AbstractArray, pos::Integer, val) -> Nothing

Add val at nzval[pos] (linear indexing), where nzval is a matrix's _scatter_storage: the write every replay sink makes.

_scatter_storage(A::AbstractMatrix) -> AbstractArray

The array that _scatter_add! writes into, indexed by the same positions _scatter_position answers. Storage and position agree type by type: A.nzval for a SparseMatrixCSC (and for SparseMatrixCSR{1}, in its extension), whose position is an nzval index; A.mirror.nzval for a device matrix (gpena/Bramble.jl#313); and A itself for every other host matrix, whose position is a linear index (A[pos]). A backend that answers nzval indices from its own _scatter_position adds its own method here too.

The host fallback returns A rather than vec(A): on a Matrix, vec builds a new array header, one allocation per sink. Any other host sparse matrix (FixedSparseCSC, say) takes that fallback too, since its _scatter_position is the linear-index one.

What the replay sinks (ReplaySink, _PairReplaySink, _StrideReplaySink, _DiagonalReplayTarget) hold instead of A, so a sparse sink carries plain vectors only and a threaded sweep can hand it across a task boundary as is (gpena/Bramble.jl#437). Each sink adds through _replay_add!, the same += as _scatter_add!, so a replay stores the same bits as before.

_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. _CoordSink 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) 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.

_test_inbounds(lin_indices, I::CartesianIndex, off_v) -> Bool

Whether off_v's test index survives, without computing the row _test_row would return. An AbsoluteRow always survives, matching _test_row.

_test_row(lin_indices, I::CartesianIndex, off_v) -> Int

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

The mirror of _trial_column. An ordinary offset is bounds-checked against the grid and answers 0 outside it; a test-side interpolation carries an AbsoluteRow instead, naming a row of the other mesh outright (gpena/Bramble.jl#263).

_test_row_unguarded(lin_indices, I::CartesianIndex, off_v) -> Int

_test_row, without the bounds check: the interior-core counterpart, matching _trial_column_unguarded.

_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::AbstractMatrix, 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.

Dispatches through _scatter_position/_scatter_add! on A's own locality: a host matrix searches/writes A directly, a device one searches/writes A.mirror (gpena/Bramble.jl#313) – either way this function itself carries no device-specific branch or parameter.

See also: allocate_system_matrix, whose positions search raises the same way on the serial path.

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 (_step_entry!, which fuses _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 the margins of _face_margins 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. The one exception is a term under the facet rule: the core walks it unmarked (_interior_term), without the zero-weight entries of the marked function's extra stencil parts. The same values land, and the recording (_coord_walk!) and every replay split the same way, so their entries still agree. ReplaySink is unaffected by that: it addresses each point by its own linear index, not by visit order (see its docstring).

For the recording and replay sinks a product term is walked factor by factor, never built into its product stencil (_factored_walk): the same entries, in the same order and to the bit, handed to the sink as they are formed.

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. The sinks are _CoordSink (pattern and positions), ReplaySink, _PairReplaySink and DiagonalReplaySink (values); the contract is in _sink_entry! and _sink_point!.

• term: The AST node whose stencil is evaluated at each point.

• sp: The test leaf. Its grid is walked, and the stencil sees its markers.

• 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!.

_ReplayTarget

What a threaded replay sweep writes through, one per (term, block) unit: a ReplaySink (flat segment), a _PairReplaySink (transposed pair) or a _DiagonalReplayTarget (1D diagonal segment). Each is immutable and only reads its recording, so one value is shared by every thread of the sweep.

See also: _replay_point!.

_replay_point!(target, term, sp, I::CartesianIndex, lin_indices, mesh_markers, row_offset::Int, col_offset::Int) -> Nothing

Replay one grid point I of one unit: evaluate term's stencil there (weights are live) and add each entry at the position target's recording holds for it, with no search.

The per-point step of the threaded replay, which a threaded sweep calls once per point of its band or colour. Correct in any visit order and from any thread, since the recording is addressed by lin_indices[I], and race-free under the caller's colouring, which keeps two concurrently-replayed points off the same entry. The guarded entry walk (_visit_entries) keeps exactly the entries the recording pass kept, interior points included.

See also: _ReplayTarget, _batch_bilinear_band_replay!.

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

Scatter one band colour of term into A across threads. A is a _ReplayTarget, whose recorded positions are written, or an action target.

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.

_sweep_coloured!(policy, point, grid_inds::CartesianIndices, strides) -> Nothing

Call point(I) once for every I in grid_inds, across threads, so that no two points closer than strides run at once. policy is the effective policy (_effective_parallel_policy), and point a _BilinearPoint or a _LinearPoint, whose _sweep_bands! and _sweep_points! methods pick the mechanism.

Bands come first: a colour is then a slab of whole rows, walked contiguously, and there are two of them however wide the stencil is. Point colouring strides every axis and costs prod(strides) barriers, which on a nine-colour form was slower than not threading at all (Dcₓ + Dcᵧ at 1000x1000: 16.6 ms threaded against 16.4 ms serial, and 8.1 ms banded). A grid too short to band (a small grid, or a stencil reaching most of the last axis) is point coloured, which has no width requirement.

DiracSource{D, P, S} <: LazyOp{D}

An AST node representing a point (Dirac delta) distribution or collection of point sources.

A DiracSource represents the linear functional:

Discretized on a grid space \(W_h\), an off-node source \(x_0\) falling into cell \(\text{cell}(idx)\) with fractional coordinates \(t \in [0, 1]^D\) distributes its strength \(S\) across the \(2^D\) surrounding cell vertices via multilinear corner weights:

The resulting load vector \(b = \text{assemble}(l)\) satisfies \(b_{idx + c} = S w_c\) with \(\sum_c b_{idx + c} = S\), reproducing exact on-node delta placement (when \(x_0 = x_I\), \(b_I = S\)) and second-order accurate functional contraction \(\ell(v_h) = S v_h(x_0)\) against discrete grid functions.

See also: dirac, innerₕ, inner₊, LinearForm.

issymmetric(a::BilinearForm) -> Bool

Whether a is symmetric by construction: a sum or scaling of terms innerₕ(L(u), L(v)), with the same L written once and applied to both the trial and test argument, and of transposed pairs innerₕ(A(u), B(v)) + innerₕ(B(u), A(v)), which may sit anywhere in the sum. The two terms of a pair must carry the same coefficients: the same object (a Ref, a grid function), or numbers whose exact products are equal ((A / 2) / 3 partners B / 6), or none on either. A term with no partner makes the answer false.

Purely structural: this walks a's expression and never assembles a matrix. It is also conservative: a term that happens to produce a symmetric matrix through some other route answers false, the same as one that is not symmetric at all.

This describes the unconstrained operator. dirichlet_bc! zeros a row without touching its column, so a matrix assembled with dirichlet is not symmetric even when issymmetric(a) is true, until symmetrize! restores it. true here is a claim about a's expression, not about whatever matrix a particular call to assemble produced.

Examples

using Bramble: D₋ₓ, D₋ᵧ, inner₊ₓ
using LinearAlgebra: issymmetric
a = form(Wₕ, Wₕ, (u, v) -> inner₊ₓ(D₋ₓ(u), D₋ₓ(v)))
issymmetric(a)  # true: the same D₋ₓ on both sides

b = form(Wₕ, Wₕ, (u, v) -> inner₊(u, D₋ₓ(v)))
issymmetric(b)  # false: different operators either side, and no transposed partner

c = form(Wₕ, Wₕ, (u, v) -> innerₕ(D₋ₓ(u), D₋ᵧ(v)) + innerₕ(D₋ᵧ(u), D₋ₓ(v)))
issymmetric(c)  # true: a transposed pair

issymmetric(Matrix(assemble(a)))                                # true
issymmetric(Matrix(assemble(a; dirichlet = :boundary)))   # false: rows zeroed, columns not

isposdef(a::BilinearForm) -> Bool

Whether a is symmetric positive semi-definite by the LᵀWL construction issymmetric checks: true only when every term is innerₕ(L(u), L(v)) (a transposed pair need not be semi-definite, so one makes the answer false) and every scaling along the way is by a positive number, which is what keeps that positivity from being flipped or collapsed. A Ref scale is read when isposdef is called, so the answer follows its value at that moment and may change after the Ref is updated. A non-real scale or coefficient, complex included, answers false.

Purely structural, like issymmetric, and for the same reason conservative: this does not prove positive-definite (which also needs L to have trivial kernel), only that the assembled matrix is symmetric positive semi-definite, enough to make cholesky worth attempting first rather than a general factorization.

Describes the unconstrained operator, exactly as issymmetric does: a matrix assembled with dirichlet needs symmetrize! after dirichlet_bc! before either symmetry or positive-definiteness holds of it, isposdef(a) being true notwithstanding.