Linear and bilinear forms

The operators in the previous tutorial act on grid functions. A form is the other half: an expression written in the test function, which Bramble assembles into the vector or the matrix a solver wants. This tutorial covers:

1. Writing a linear form and assembling its vector.
2. Refilling that vector inside a time loop, and why the pattern is built once.
3. Contracting a form against a grid function without building a vector at all.
4. Bilinear forms and the system matrix.
5. Imposing Dirichlet conditions, and solving a Poisson problem.
6. Coupled systems, where a term names which block it belongs to, including a term that reads from a leaf built over a different mesh.
7. Threaded assembly, and why it needs no locks.
8. Restricting a term to part of the mesh.
9. What form(...) simplifies automatically, and how to write a form so it can.

Every number below was produced by the code shown.

1. A form is an expression in the test function

A linear form is a function of one argument, and that argument stands for the test function rather than for any particular grid function:

using Bramble
using SparseArrays

Ωₕ = mesh(domain(interval(0.0, 1.0)), 33, true)
Wₕ = gridspace(Ωₕ)
fₕ = Rₕ(Wₕ, x -> sin(π * x))

l = form(Wₕ, v -> innerₕ(fₕ, v))

LinearForm {1D, Float64}:
  Test space: ScalarGridSpace{1D, Float64, 33 dofs}
  Vector: 33

Nothing has been computed. v is symbolic, so innerₕ(fₕ, v) builds a description of

\[\ell(v) = (f_h, v)_h\]

and the form stores the expression, not a vector. Building one is free, which is what makes it reasonable to write a form inside a function that is called repeatedly.

The asymmetry is worth naming early, because it decides what an expression costs. The source side is eager and the test side is symbolic: D₋ₓ(fₕ) computes a grid function, while D₋ₓ(v) adds a node to an expression. So a term is as cheap as its test side is symbolic, however elaborate the coefficient in front of it.

2. Assembling, and refilling

assemble allocates the vector and fills it:

b = assemble(l)
length(b), sum(b)

(33, 0.6361083632808496)

In a time loop the allocation is the part worth avoiding. assemble! refills a vector that already exists with zero allocations:

assemble!(b, l)
sum(b)

0.6361083632808496

Forms store their resolved abstract syntax tree (ast) directly upon creation, retaining direct references to the underlying coefficient arrays.

Live grid coefficients and dynamic scalars

• Grid functions: Overwrite a coefficient element in-place with Rₕ!(fₕ, ...) or parent(fₕ) .= ... between steps, and the next assemble!(b, l) evaluates the new values live with 0 bytes allocated, without needing to reconstruct the form.
• Scalar coefficients: Constant scalar factors can be written directly as plain numbers (e.g. 2.5 * innerₕ(fₕ, v)). A Ref(val) is only needed when you want a dynamic scalar coefficient that changes across loop iterations:

α = Ref(1.0)
l_dyn = form(Wₕ, v -> α * innerₕ(fₕ, v))
b_dyn = assemble(l_dyn)
α[] = 2.0
assemble!(b_dyn, l_dyn) # 0 bytes allocated, live 2x scaling
sum(b_dyn) ≈ 2 * sum(b)

true

3. Contracting without a vector

Often the vector is not wanted, only the number $\ell(v_h)$. A form is callable, and takes that shortcut:

oneₕ = Rₕ(Wₕ, x -> 1.0)
l(oneₕ), sum(b)

(0.6361083632808495, 0.6361083632808496)

Against the all-ones grid function a linear form is the sum of its assembled vector, which is the check above. The difference is that l(oneₕ) builds no vector: it contracts as it walks the grid, and allocates nothing at all. Where the result is a scalar, prefer it.

evaluate! is the middle case — it wants the assembled vector and the number, so it takes a scratch vector, fills it, and returns the contraction:

scratch = zeros(length(b))
evaluate!(scratch, l, oneₕ)

0.6361083632808496

A form contracts against an element of its test space, never against a bare vector. The length of a vector says nothing about whether its blocks line up with the components a form routes to, so accepting one would make a coupled mismatch silent rather than loud.

4. Bilinear forms and the system matrix

A bilinear form takes two symbolic arguments, trial first and test second:

a = form(Wₕ, Wₕ, (u, v) -> inner₊ₓ(D₋ₓ(u), D₋ₓ(v)))
A = assemble(a)
size(A), nnz(A)

((33, 33), 97)

which is the stiffness matrix of

\[a(u, v) = (D_{-x} u, D_{-x} v)_{+x}.\]

Ninety-seven nonzeros in a 33-by-33 matrix is the tridiagonal band, and the band is the point: the sparsity pattern follows from the stencil, so it is known before any value is computed.

That is what makes the two-step idiom worth using. allocate_system_matrix builds the pattern and nothing else; assemble! then fills a matrix whose structure already exists:

A2 = allocate_system_matrix(a)
assemble!(A2, a)
A2 ≈ A

true

Inside a time loop, build the pattern once outside it and call assemble! within. Refilling a matrix whose pattern is fixed allocates nothing, where assemble allocates a new matrix every step.

a above is innerₕ(L(u), L(v)) with the same D₋ₓ on both sides, which is symmetric — and, since the quadrature weight inner₊ₓ carries is positive, positive semi-definite — purely by that construction. issymmetric/isposdef answer this from the expression alone, without assembling anything:

using LinearAlgebra: issymmetric, isposdef
issymmetric(a), isposdef(a)

(true, true)

c = form(Wₕ, Wₕ, (u, v) -> inner₊(u, D₋ₓ(v)))
issymmetric(c)  # different operators either side — not this pattern

false

Knowing this before assembling is what makes a positive answer worth something: it says cholesky is worth trying on the result rather than a general factorization, at a cost — a few nanoseconds, against tens of microseconds to assemble even this small a matrix — close enough to free that there is no reason not to check.

5. Dirichlet conditions, and a Poisson problem

Boundary conditions come in two pieces, because a matrix and a right-hand side need different things done to them. Labels on the domain say where:

Ω = domain(interval(0.0, 1.0), :left => :left, :right => :right)
Ωd = mesh(Ω, 33, true)
Wd = gridspace(Ωd)

fd = Rₕ(Wd, x -> π^2 * sin(π * x))
ad = form(Wd, Wd, (u, v) -> inner₊ₓ(D₋ₓ(u), D₋ₓ(v)))
ld = form(Wd, v -> innerₕ(fd, v))

Ad = assemble(ad)
bd = assemble(ld)

dirichlet_constraints records the values, dirichlet_bc! applies them — to the matrix by replacing the constrained rows, and to the vector by writing the boundary values in. dirichlet_constraints takes the mesh (or a Domain/grid space) directly — no need to extract the underlying CartesianProduct first:

bcs = dirichlet_constraints(Ωd, :left => (x -> 0.0), :right => (x -> 0.0))
dirichlet_bc!(Ad, Ωd, :left, :right)
dirichlet_bc!(bd, Ωd, bcs, :left, :right)

Solving $-u'' = \pi^2 \sin(\pi x)$ with $u(0) = u(1) = 0$ gives $u = \sin(\pi x)$:

uh = Ad \ bd
exact = Rₕ(Wd, x -> sin(π * x))
maximum(abs, uh .- parent(exact))

0.0008035776793697824

Eight parts in ten thousand on 33 points, which is second order behaving itself.

using AppleAccelerate
Aa = AppleAccelerate.AASparseMatrix(Ad)
factor = AppleAccelerate.SparseFactor(AppleAccelerate.SparseFactorizationCholesky, Aa.matrix)
uh_fast = copy(bd)
AppleAccelerate.SparseSolve(factor, uh_fast)   # solves in place

Imposing conditions by replacing rows destroys symmetry, and a symmetric solver will want it back. symmetrize! moves the constrained columns onto the right-hand side, restoring symmetry and leaving the solution unchanged:

issymmetric(ad)          # true — the form is symmetric by construction, before any boundary condition

true

issymmetric(Matrix(Ad))  # false — dirichlet_bc! zeroed rows, not columns

false

symmetrize!(Ad, bd, Ωd, :left, :right)
issymmetric(Matrix(Ad))  # true again, and the solution above is unchanged

true

issymmetric(ad) is a claim about the expression ad, not about any one matrix that gets assembled from it — it says nothing about what dirichlet_bc! alone leaves behind, which is exactly why the middle line above answers false even though the first one answers true.

6. Coupled systems

A composite space stacks copies of a space, and a form over one addresses its blocks by component: u[1] (or u(1)) is the trial function of the first block, v[2] (or v(2)) the test function of the second. Both functor indexing u(i) and standard bracket indexing u[i] are supported on trial and test functions, as well as on compound operators ((D₋ₓ(u))[i]).

In addition, trial and test functions support tuple destructuring via components or direct iteration:

Vₕ = Wₕ^Val(2)
ac = form(Vₕ, Vₕ, (u, v) -> begin
    u₁, u₂ = components(u)
    v₁, v₂ = components(v)
    innerₕ(u₁, v₁) + inner₊ₓ(D₋ₓ(u₂), D₋ₓ(v₂))
end)
Ac = assemble(ac)
size(Ac)

(66, 66)

Sixty-six by sixty-six: two blocks of 33, assembled into one matrix. A term naming u[i] and v[j] lands in block $(j, i)$, so off-diagonal coupling is written the same way — innerₕ(u[1], v[2]) fills the block that couples the first unknown to the second equation.

Component indices are checked against the number of blocks at form construction time: accessing u[3] or u(3) on a 2-component space raises an immediate ArgumentError. Furthermore, a term must name both components or neither:

form(Vₕ, Vₕ, (u, v) -> innerₕ(u[1], v))   # ArgumentError

Naming one and leaving the other open has no reading as mathematics — the term would belong to every equation at once — so it is refused rather than guessed at. Naming neither is fine and means the diagonal, applied to every block.

Constraining one block, leaving another free

dirichlet on its own binds to every leaf sharing the named marker — fine when every block wants the same treatment, not when they don't. A Stokes-style system prescribing velocity while leaving pressure unconstrained needs dirichlet_components too: 1-based leaf positions, the same order u(1)/u(2) addressing already uses.

Ωc = domain(interval(0.0, 1.0), :left => :left, :right => :right)
Ωdc = mesh(Ωc, 21, true)
Vc = gridspace(Ωdc)^Val(2)           # 1: velocity-like, 2: pressure-like
ac2 = form(Vc, Vc, (u, v) -> innerₕ(u(1), v(1)) + innerₕ(u(2), v(2)))
Ac2 = assemble(ac2; dirichlet = (:left, :right), dirichlet_components = 1)

Block 1 (rows 1:21) has its boundary rows pinned; block 2 is untouched — still the plain assembled operator, no rows replaced at all. Leaving dirichlet_components at its default (nothing) applies the labels to every leaf, exactly as before this keyword existed; call assemble!/dirichlet_bc! again with a different dirichlet/dirichlet_components pair to constrain another block differently.

Interpolating between the leaves of a heterogeneous composite space

The composite spaces above stack copies of one space — every leaf shares a mesh. A composite space can also be built directly from a tuple of leaves over different meshes, and then a term coupling two leaves needs a way to move a value from one leaf's grid to the other's: πₕ, one argument fewer than the numeric πₕ/πₕ! pair (see the operators tutorial for the numeric side and a diagram of the interpolant itself) — the same name, told apart by dispatch rather than a different one.

πₕ(uₕ) reads exactly like any other source — it is one, an AST leaf wrapping x -> interpolate_at(uₕ, x) — so it composes with D₋ₓ, M₋ₓ, and the rest the same way sin, a VectorElement, or any other source does, and can sit on the left of innerₕ inside a coupled form:

Ωbig = mesh(domain(box((0.0, 0.0), (1.0, 1.0))), (8, 8), (true, true))
Ωsmall = mesh(domain(box((0.0, 0.0), (1.0, 1.0))), (4, 4), (true, true))
Wbig, Wsmall = gridspace(Ωbig), gridspace(Ωsmall)
Vh = CompositeGridSpace((Wbig, Wsmall))
uv = Rₕ(Vh, (x -> 0.0, x -> x[1] + x[2]))   # only the small leaf (2) carries data

lh = form(Vh, v -> innerₕ(πₕ(uv(2)), v(1)) + innerₕ(D₋ₓ(πₕ(uv(2))), D₋ₓ(v(1))))
b = assemble(lh)

# the differenced term is not a no-op: dropping it changes the answer
b_plain = assemble(form(Vh, v -> innerₕ(πₕ(uv(2)), v(1))))
maximum(abs, b .- b_plain)

0.14285714285714285

The two terms land in the same block (leaf 1, Wbig) even though the source they read from lives on leaf 2's own, coarser mesh — πₕ is what makes that a well-posed expression rather than a size mismatch. This is exactly what makes a heterogeneous composite space useful for more than indexing: leaf 2 can represent one field at a resolution the problem calls for, and a term over leaf 1 can still read it.

That last line is the check worth keeping, not length(b) == ndofs(Vh). An earlier draft of this page showed the shape instead, and the differenced term contributed exactly zero: an operator wrapped around a source had its offsets discarded, so D₋ₓ's +s/h and −s/h cancelled. The page was green either way, because a zero vector has the right length and is perfectly finite. A worked example should show what it computes.

An operated source is worth a word on what it means. innerₕ(D₋ₓ(f), v) is $\sum_i |\square_i| \, (D_{-x}f)_i \, v_i$ — the operator acts on the source, producing another grid function, which is then integrated against the test function. It agrees entry for entry with applying the numeric operator first: assemble(form(Wₕ, v -> innerₕ(D₋ₓ(fₕ), v))) equals parent(D₋ₓ(fₕ)) .* weights(Wₕ, Innerh()). That equivalence is what test/form/source_operators.jl pins, for every operator, against the numeric layer.

A bilinear term coupling two leaves over different meshes is a different matter, and it is refused:

assemble(form(Vh, Vh, (u, v) -> innerₕ(u(2), v(1))))   # ArgumentError: ... over different meshes ...

(The refusal is raised when the matrix is built, not when the form is written: form resolves the expression, and which leaves a term couples is a question about the spaces it is assembled against. allocate_system_matrix refuses it too, so neither entry point can be reached around.)

A coupled block is assembled by walking the test leaf's grid and reading the trial column out of that same index space, so it needs the two leaves to agree on what an index means. Two leaves over meshes of different sizes do not: index (3, 3) on an 8×8 grid and on a 4×4 grid name different points, and nothing in the term says how to get from one to the other. So there is no assembly to give, and the error says so rather than guessing — in one direction it used to overrun the trial block and throw from deep inside sparse!, and in the other it quietly filled in-range but wrong columns.

Coupling leaves that share a mesh is unaffected, which is every composite space built by repeating one space (Wₕ^Val(2)), including off-diagonal blocks.

7. Threading, chosen once on the backend

Rₕ!, avgₕ!, gridspace construction and form assembly all thread the same way: Serial() or Parallel(), chosen once when the backend is built, rather than decided per call. See the backend tutorial for backend construction in general — vector/matrix types included — and for how to choose between the two policies; this section only covers what the choice means for assembly specifically.

Wₕ_par = gridspace(mesh(domain(interval(0.0, 1.0)), 33, true;
    backend = backend(policy = Parallel())))
execution_policy(Wₕ_par)

Parallel()

There is no automatic size threshold. A Parallel() backend threads every eligible call, however small, however often a time loop repeats it — asking for Parallel() and getting it is the point, rather than a heuristic guessing on the caller's behalf whether a given call is big enough to be worth it. Pick Serial() (the default backend() already is) for small, frequently repeated calls instead.

assemble!/assemble follow test_space(form)'s (or, for a BilinearForm, trial_space's) policy directly, so the ordinary entry points already thread when the backend says to:

l_par = form(Wₕ_par, v -> innerₕ(Rₕ(Wₕ_par, x -> sin(π * x)), v))
b_par = assemble(l_par)      # threads, because Wₕ_par's backend says Parallel()

assemble_parallel! still exists underneath, as a lower-level entry point that always threads regardless of the backend's policy — useful for a one-off forced comparison or a benchmark, not the everyday call:

bp = similar(b)
assemble_parallel!(bp, l)    # always threads, whatever l's own backend says
bp ≈ b

false

Threading takes no locks, and needs none. Assembly partitions the grid by stride: the offsets a stencil reaches give the width of the footprint one point writes, and two points separated by at least that width cannot overlap. Points sharing a stride are therefore written concurrently with nothing to coordinate.

The common case is one colour. A form whose test argument carries no difference — innerₕ(fₕ, v) above — reaches only its own point, so the stride is 1 in every direction and the whole grid is swept in a single flat parallel pass. A gradient term in two dimensions reaches one point back along each axis, giving four colours swept in turn.

Whether threading pays depends on the size, and not always in the obvious direction: assembly is memory-bound, so the gain flattens well before the thread count does. The benchmarks page carries the measurements — that is what should decide which policy a backend is built with, not a guess.

8. Restricting a term to part of the mesh

innerₕ, inner₊ and the directional products all take a markers keyword, restricting the sum to the union of the regions the labels name — the same idea as restrict_to, spelled at the call site rather than wrapping an argument:

a_left = form(Wd, Wd, (u, v) -> innerₕ(u, v; markers = (:left,)))
size(assemble(a_left))

(33, 33)

Every mesh also carries :boundary and :interior automatically, computed from its own shape rather than needing any label set up in domain(...):

a_boundary = form(Wd, Wd, (u, v) -> innerₕ(u, v; markers = (:boundary,)))
size(assemble(a_boundary))

(33, 33)

This is a masked sum of the existing cell measures — not a surface integral, and the two are not interchangeable; a masked innerₕ scales like h and vanishes under refinement, where a true boundary integral does not. markers is for the former; a Neumann or Robin term needing the latter is a separate, not-yet-built piece (inner_Γ).

A marker that does not exist anywhere the term reaches is a loud error rather than a silent all-zero contribution — RegionRestriction's own per-point check cannot tell "nothing here is marked" from "no such marker", so this is caught once, before assembling anything:

try
    assemble(form(Wd, Wd, (u, v) -> innerₕ(u, v; markers = (:nope,))))
catch e
    println(e)
end

ArgumentError("the marker :nope is not defined on the form's space. A marker named in restrict_to or markers = (...) must exist on every space a term reaches; if it is only defined on some of a composite space's leaves, write the term per component instead, one innerₕ(u(i), v(i)) per leaf with that leaf's own markers, rather than one term naming a marker not every leaf it reaches has.")

On a composite space, a marker used without naming a component reaches every diagonal block, and has to exist on every leaf that reaches — write the term per component, each with its own markers, if it does not.

9. What form(...) simplifies automatically

form(Wₕ, Vₕ, f)/form(Wₕ, f) resolve the expression once and then run it through an algebraic simplification pass before storing it. This matters for how you write a form, because the assembler routes a form's summands one at a time: every + in the expression is a separate sweep over the mesh, so an expression with fewer top-level summands assembles faster, for exactly the same matrix or vector.

Most of the rewrites touch only +, * and / — never the operators inside them (D₋ₓ, inner₊, innerₕ, and the rest) — so they apply to whatever is built from those, coupled systems and restricted terms included. A few reach one layer deeper, into an inner product's own arguments and into shift_op, because leaving them out would mean either a correctness gap or a documented dead end; §"What stays as written, and why" below draws the exact line.

Identical terms combine into one term

Two summands that are the same expression merge into a single scaled term, rather than being routed and assembled separately:

a_dup = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + innerₕ(u, v))
Bramble.resolve_form_ast(a_dup)  # 2 * innerₕ(u, v) — one term, not two

Bramble.OperatorScale{1, Int64, Bramble.BilinearProduct{1, Bramble.InnerH, Bramble.TrialFunction{1, 1}, Bramble.TestFunction{1, 1}}}(2, Bramble.BilinearProduct{1, Bramble.InnerH, Bramble.TrialFunction{1, 1}, Bramble.TestFunction{1, 1}}(Bramble.TrialFunction{1, 1}(), Bramble.TestFunction{1, 1}()))

Matrix(assemble(a_dup)) ≈ 2 .* Matrix(assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v))))

true

A form built up piece by piece — accumulating one contribution per physical effect, some of which may coincide — pays nothing for the duplication once assembled: write the terms separately if that is the clearer expression of the model, rather than checking by hand whether two of them happen to repeat.

A shared scalar factors out of a sum

c * A + c * B, for two different A and B, becomes c * (A + B): still two operators to evaluate, but one routed term instead of two.

Ω2 = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (12, 10), (true, true))
W2 = gridspace(Ω2)

a_split = form(W2, W2, (u, v) -> 2 * inner₊ₓ(D₋ₓ(u), D₋ₓ(v)) + 2 * inner₊ᵧ(D₋ᵧ(u), D₋ᵧ(v)))
Bramble.resolve_form_ast(a_split)  # 2 * (inner₊ₓ(...) + inner₊ᵧ(...)) — one routed term

Bramble.OperatorScale{2, Int64, Bramble.OperatorAdd{2, Bramble.BilinearProduct{2, Bramble.InnerPlus{1}, Bramble.BackwardDifference{2, 1, Bramble.TrialFunction{2, 1}}, Bramble.BackwardDifference{2, 1, Bramble.TestFunction{2, 1}}}, Bramble.BilinearProduct{2, Bramble.InnerPlus{2}, Bramble.BackwardDifference{2, 2, Bramble.TrialFunction{2, 1}}, Bramble.BackwardDifference{2, 2, Bramble.TestFunction{2, 1}}}}}(2, Bramble.OperatorAdd{2, Bramble.BilinearProduct{2, Bramble.InnerPlus{1}, Bramble.BackwardDifference{2, 1, Bramble.TrialFunction{2, 1}}, Bramble.BackwardDifference{2, 1, Bramble.TestFunction{2, 1}}}, Bramble.BilinearProduct{2, Bramble.InnerPlus{2}, Bramble.BackwardDifference{2, 2, Bramble.TrialFunction{2, 1}}, Bramble.BackwardDifference{2, 2, Bramble.TestFunction{2, 1}}}}(Bramble.BilinearProduct{2, Bramble.InnerPlus{1}, Bramble.BackwardDifference{2, 1, Bramble.TrialFunction{2, 1}}, Bramble.BackwardDifference{2, 1, Bramble.TestFunction{2, 1}}}(Bramble.BackwardDifference{2, 1, Bramble.TrialFunction{2, 1}}(Bramble.TrialFunction{2, 1}()), Bramble.BackwardDifference{2, 1, Bramble.TestFunction{2, 1}}(Bramble.TestFunction{2, 1}())), Bramble.BilinearProduct{2, Bramble.InnerPlus{2}, Bramble.BackwardDifference{2, 2, Bramble.TrialFunction{2, 1}}, Bramble.BackwardDifference{2, 2, Bramble.TestFunction{2, 1}}}(Bramble.BackwardDifference{2, 2, Bramble.TrialFunction{2, 1}}(Bramble.TrialFunction{2, 1}()), Bramble.BackwardDifference{2, 2, Bramble.TestFunction{2, 1}}(Bramble.TestFunction{2, 1}()))))

a_x = form(W2, W2, (u, v) -> inner₊ₓ(D₋ₓ(u), D₋ₓ(v)))
a_y = form(W2, W2, (u, v) -> inner₊ᵧ(D₋ᵧ(u), D₋ᵧ(v)))
Matrix(assemble(a_split)) ≈ 2 .* (Matrix(assemble(a_x)) .+ Matrix(assemble(a_y)))

true

An isotropic operator written out direction by direction — the common way to build one before reaching for a name like inner₊/∇₋ₕ that already sums over every direction — assembles as cheaply as writing it the terser way by hand.

A zero-scaled term leaves no trace

A term scaled by the literal number 0 — a coefficient set to zero for a particular run, common in continuation methods and IMEX schemes toggling a physical effect on and off — contributes nothing to the sparsity pattern, rather than reserving space for the stencil it would otherwise have:

a_full = form(W2, W2, (u, v) -> innerₕ(u, v) + 0.0 * inner₊ₓ(D₋ₓ(u), D₋ₓ(v)))
a_mass = form(W2, W2, (u, v) -> innerₕ(u, v))
nnz(assemble(a_full)) == nnz(assemble(a_mass))  # the stiffness term left no entries at all

true

Without this, the zero-scaled stiffness term would still reserve its full band in the pattern — nonzero positions holding the value 0.0 — which costs both memory and a wasted sweep computing them. Toggling a term off is free to leave in the expression; there is no need to branch in Julia code around it.

A dynamic coefficient (a Ref, §2's "Live grid coefficients and dynamic scalars") combines and factors the same way a static number does, and keeps tracking its own updates afterwards — the rewrite only ever moves the Ref around, never reads the value inside it:

β = Ref(1.0)
a_ref = form(Wₕ, Wₕ, (u, v) -> β * innerₕ(u, v) + β * innerₕ(u, v))
Matrix(assemble(a_ref)) ≈ 2 .* Matrix(assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v))))

true

β[] = 3.0
Matrix(assemble(a_ref)) ≈ 6 .* Matrix(assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v))))

true

Two different Refs, or a Ref alongside a plain number, never combine — a rewrite that assumed two independent dynamic coefficients were the same value would be a correctness bug the first time they diverged, so it is not attempted; write the shared coefficient as one Ref, used on every term it scales, if two terms are meant to move together.

A scalar inside an inner product's argument is lifted back out

The three rules above stop at innerₕ/inner₊/... itself — but a scalar written inside one of their arguments is lifted back out to wrap the whole product, exposing it to exactly those rules:

a_hidden = form(Wₕ, Wₕ, (u, v) -> innerₕ(2 * D₋ₓ(u), D₋ₓ(v)) + innerₕ(3 * D₋ₓ(u), D₋ₓ(v)))
Bramble.resolve_form_ast(a_hidden)  # 5 * innerₕ(D₋ₓ(u), D₋ₓ(v)) — one term, not two

Bramble.OperatorScale{1, Int64, Bramble.BilinearProduct{1, Bramble.InnerH, Bramble.BackwardDifference{1, 1, Bramble.TrialFunction{1, 1}}, Bramble.BackwardDifference{1, 1, Bramble.TestFunction{1, 1}}}}(5, Bramble.BilinearProduct{1, Bramble.InnerH, Bramble.BackwardDifference{1, 1, Bramble.TrialFunction{1, 1}}, Bramble.BackwardDifference{1, 1, Bramble.TestFunction{1, 1}}}(Bramble.BackwardDifference{1, 1, Bramble.TrialFunction{1, 1}}(Bramble.TrialFunction{1, 1}()), Bramble.BackwardDifference{1, 1, Bramble.TestFunction{1, 1}}(Bramble.TestFunction{1, 1}())))

Matrix(assemble(a_hidden)) ≈ 5 .* Matrix(assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(D₋ₓ(u), D₋ₓ(v)))))

true

This is not only a routing question. issymmetric/isposdef (§4) recognise innerₕ(L(u), L(v)) — the same operator on both sides — structurally, and a scalar sitting inside one argument used to hide that shape from the check, because 2 * D₋ₓ(u) and D₋ₓ(v) are different node types even though the pattern is exactly the symmetric one:

issymmetric(form(Wₕ, Wₕ, (u, v) -> innerₕ(2 * D₋ₓ(u), D₋ₓ(v))))

true

Write the scalar wherever reads best — 2 * innerₕ(D₋ₓ(u), D₋ₓ(v)) and innerₕ(2 * D₋ₓ(u), D₋ₓ(v)) now assemble, and are checked for symmetry, identically.

A component-mixing sum inside one inner product

innerₕ(fₕ, v(1) + v(2)) names two different components of a coupled test space inside one product — asking, in effect, for fₕ's contribution to land in two different equations at once. There is no single routed term that means that, so this distributes into two, the same shape as writing them separately:

Vₕ = Wₕ^Val(2)
fₕ = Rₕ(Wₕ, x -> sin(π * x[1]))

l_mixed = form(Vₕ, v -> innerₕ(fₕ, v(1) + v(2)))
Bramble.resolve_form_ast(l_mixed)  # innerₕ(fₕ, v(1)) + innerₕ(fₕ, v(2))

Bramble.OperatorAdd{1, Bramble.LinearProduct{1, Bramble.InnerH, Bramble.SourceVector{1, Vector{Float64}}, Bramble.IndexedTestFunction{1}}, Bramble.LinearProduct{1, Bramble.InnerH, Bramble.SourceVector{1, Vector{Float64}}, Bramble.IndexedTestFunction{1}}}(Bramble.LinearProduct{1, Bramble.InnerH, Bramble.SourceVector{1, Vector{Float64}}, Bramble.IndexedTestFunction{1}}(Bramble.SourceVector{1, Vector{Float64}}([0.0, 0.0980171403295606, 0.19509032201612825, 0.29028467725446233, 0.3826834323650898, 0.47139673682599764, 0.5555702330196022, 0.6343932841636455, 0.7071067811865475, 0.7730104533627369  …  0.7730104533627371, 0.7071067811865476, 0.6343932841636455, 0.5555702330196022, 0.47139673682599786, 0.3826834323650899, 0.2902846772544624, 0.1950903220161286, 0.09801714032956083, 1.2246467991473532e-16]), Bramble.IndexedTestFunction{1}(1)), Bramble.LinearProduct{1, Bramble.InnerH, Bramble.SourceVector{1, Vector{Float64}}, Bramble.IndexedTestFunction{1}}(Bramble.SourceVector{1, Vector{Float64}}([0.0, 0.0980171403295606, 0.19509032201612825, 0.29028467725446233, 0.3826834323650898, 0.47139673682599764, 0.5555702330196022, 0.6343932841636455, 0.7071067811865475, 0.7730104533627369  …  0.7730104533627371, 0.7071067811865476, 0.6343932841636455, 0.5555702330196022, 0.47139673682599786, 0.3826834323650899, 0.2902846772544624, 0.1950903220161286, 0.09801714032956083, 1.2246467991473532e-16]), Bramble.IndexedTestFunction{1}(2)))

l_split = form(Vₕ, v -> innerₕ(fₕ, v(1)) + innerₕ(fₕ, v(2)))
assemble(l_mixed) ≈ assemble(l_split)

true

Before this rule, l_mixed assembled to an ArgumentError naming the two mismatched components rather than a vector — writing the sum by hand, as l_split does, was the only way to couple one source to two equations. Both spellings work now; write whichever reads better at the call site.

This is the one rule that can turn a single sweep back into two rather than the reverse: a sum naming the same component on both sides (v(1) + D₋ₓ(v(1)), say) is left as the one term it already was, since nothing forces it apart — only a genuine mismatch, which had no valid single-term routing to begin with, triggers the split. A coefficient wrapping a mixed sum distributes along with it, for the same reason: 2 * innerₕ(fₕ, v(1) + v(2)) assembles 2 * innerₕ(fₕ, v(1)) + 2 * innerₕ(fₕ, v(2)), not an OperatorScale hiding the same unroutable shape from view.

Nested grid-function scalings fuse into one array

u_h * (v_h * A) — two grid functions scaling the same operator, one wrapping the other — precomputes their elementwise product once, at construction, rather than reading both arrays at every point of every assembly:

vₕ = Rₕ(Wₕ, x -> x[1] + 1.0)
wₕ = Rₕ(Wₕ, x -> 2.0)
a_fused = form(Wₕ, Wₕ, (u, v) -> vₕ * (wₕ * innerₕ(u, v)))
Bramble.resolve_form_ast(a_fused).grid_function ≈ parent(vₕ) .* parent(wₕ)

true

Unlike the rules above, this one is not free: fusing two coefficients into one is an elementwise multiply over the whole array, paid once when the form is built rather than once per assembly. Worth it whenever the form outlives a single assembly (a time loop, a residual evaluated repeatedly), which is the common case; if a form is truly built and assembled once, the two scalings would have cost the same either way.

Two nested shifts combine, and a zero shift disappears

shift_op composes the way integer addition does: two shifts along the same dimension combine their amounts, and a net shift of zero is the identity — including a shift undone by its own inverse:

using Bramble: shift_op
a_shift = form(Wₕ, Wₕ, (u, v) -> innerₕ(shift_op(shift_op(u, 1, 2), 1, -2), v))
Bramble.resolve_form_ast(a_shift)  # innerₕ(u, v) — the two shifts cancelled

Bramble.BilinearProduct{1, Bramble.InnerH, Bramble.TrialFunction{1, 1}, Bramble.TestFunction{1, 1}}(Bramble.TrialFunction{1, 1}(), Bramble.TestFunction{1, 1}())

A shift along a different dimension never combines with one it wraps: Shift_x and Shift_y are different operations, not two amounts of the same one, so nesting them stays exactly as written.

What stays as written, and why

Every rule above stops at BilinearProduct/LinearProduct/ShiftNode: none of it descends into a difference, an average, a jump, a restriction or an interpolation. A scalar or a shift buried one layer further in —

# NOT lifted: the `2` sits inside D₋ₓ's own argument, one layer past where this pass looks
form(Wₕ, Wₕ, (u, v) -> innerₕ(D₋ₓ(2 * u), v))

BilinearForm {1D, Float64}:
  Trial space: ScalarGridSpace{1D, Float64, 33 dofs}
  Test space: ScalarGridSpace{1D, Float64, 33 dofs}  (same as trial)
  Matrix: 33 × 33
  Symmetric: no

— is invisible to it, the same way innerₕ(2 * u, v) used to be before the rule above: write the scalar where the pass can see it, 2 * innerₕ(D₋ₓ(u), v) or innerₕ(2 * D₋ₓ(u), v), rather than nested inside the difference's own argument.

πₕ(Wsrc, u) is never folded away, even when Wsrc happens to be exactly the space u is assembled against — a rewrite that could fire would need to know the trial space a term is about to be assembled into, which an expression built before form sees any space does not have. This is rarely a real cost: coupling two leaves that already share a mesh needs no πₕ at all (§6, "Interpolating between the leaves of a heterogeneous composite space").

Where to go next

The internals page on forms documents the colouring and the stencil algebra underneath all of this, including how the matrix path colours on the test-side span alone, and the exact rewrite rules behind §9's automatic simplification.