Linear and bilinear forms#

What you will learn. How to write a linear or bilinear form, assemble it into a vector or a matrix, impose Dirichlet conditions and solve a Poisson problem, restrict a term to a marked part of the mesh, integrate over a boundary face, add point sources, and what form simplifies.

What you need first. The mesh tutorial, the space tutorial for grid spaces and elements, and the operators tutorial for D₋ₓ, innerₕ and inner₊.

Where next. Ghost points for Neumann and Robin conditions, and coupled systems, where one form addresses several unknowns.

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. The running problem is the one-dimensional Poisson equation on the unit interval.

A form is an expression in the test function#

A linear form is a function of one argument, the test function:

using Bramble
using SparseArrays
import Bramble: allocate_system_matrix, evaluate!, reaction

Ωₕ = 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
  Expression: innerₕ(vec, v)

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

and the form stores the expression only. Building one is free, so it is reasonable to write a form inside a function that is called repeatedly.

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.

Assembling and refilling#

assemble allocates the vector and fills it:

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

(33, 0.6361083632808495)

In a time loop, assemble! refills an existing vector with zero allocations:

assemble!(b, l)
sum(b)

0.6361083632808495

Live grid coefficients and dynamic scalars#

A form keeps direct references to its coefficient arrays. Overwrite a grid function in place with Rₕ!(fₕ, ...) or parent(fₕ) .= ... between steps, and the next assemble!(b, l) reads the new values with no allocation and no rebuilt form.

A constant factor is written as a plain number, as in 2.5 * innerₕ(fₕ, v). A Ref is for a scalar that changes across iterations:

α = Ref(1.0)
l_dyn = form(Wₕ, v -> α * innerₕ(fₕ, v))
b_dyn = assemble(l_dyn)
α[] = 2.0
assemble!(b_dyn, l_dyn)
sum(b_dyn) ≈ 2 * sum(b)

true

Contracting without a vector#

To get only the number \(\ell(v_h)\), call the form:

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

(0.6361083632808495, 0.6361083632808495)

Against the all-ones grid function a linear form is the sum of its assembled vector. l(oneₕ) builds no vector and allocates nothing; prefer it when the result is a scalar.

evaluate! is the middle case. It takes a scratch vector, fills it with the assembled vector and returns the contraction:

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

0.6361083632808496

A bare vector is refused: its length does not show whether its blocks match the components the form routes to.

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₊(∇ₕ(u), ∇ₕ(v), :x))
A = assemble(a)
size(A), nnz(A)

((33, 33), 97)

This is the stiffness matrix of

Ninety-seven nonzeros in a 33-by-33 matrix is the tridiagonal band. The sparsity pattern follows from the stencil, so it is known before any value is computed. allocate_system_matrix builds the pattern and nothing else, and assemble! then fills it:

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

true

Build the pattern once outside a time loop and call assemble! inside it.

a pairs the same D₋ₓ on both sides, so it is symmetric. The weight inner₊ₓ carries is positive, so it is also positive semi-definite. issymmetric and isposdef answer 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, ∇ₕ(v)))
issymmetric(c)  # different operators either side

false

A positive answer means cholesky is the factorization to try first. isposdef answers for semi-definiteness: this matrix has the constants in its kernel, so it is singular until the Dirichlet conditions below are applied.

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₊(∇ₕ(u), ∇ₕ(v), :x))
ld = form(Wd, v -> innerₕ(fd, v))

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

dirichlet_constraints records the values and 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, a Domain or a grid space directly.

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

The error is 8e-4 on 33 points.

The figure draws uh over the exact solution, and below it the pointwise error with its own scale:

HTML(svg)

Imposing conditions by replacing rows destroys symmetry, and a symmetric solver wants it back. symmetrize! moves the constrained columns onto the right-hand side, which restores symmetry and leaves the solution unchanged:

issymmetric(ad)          # true: the form is symmetric before any boundary condition

true

issymmetric(Matrix(Ad))  # false: dirichlet_bc! replaced the constrained rows

false

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

true

issymmetric(ad) describes the expression, so it stays true after dirichlet_bc! changes the matrix.

To constrain one block of a coupled system and leave another free, see Coupled systems.

Boundary fluxes#

By the time uh exists, dirichlet_bc! has overwritten the constrained rows of Ad, and the flux they carried is gone. reaction recovers it. It reassembles the unconstrained ad and ld and reads the flux off the residual A*uh - F, which is approximately zero on every free row and, on a constrained row, exactly the flux the condition had to supply:

uh_elt = element(Wd, uh)
reaction(ad, ld, uh_elt; marker = :left), reaction(ad, ld, uh_elt; marker = :right)

(3.139068950903201, 3.139068950903206)

Both come out near \(\pi\). For \(u = \sin(\pi x)\) the flux \(-u'\) leaving the domain is \(\pi\) at each end, and the two together recover the discrete net source \((f_h, 1)_h\) to round-off at any resolution, which approaches \(2\pi\) at second order (\(6.2781\) here). The sign convention is in the docstring of reaction; reaction_density gives the pointwise quantity, suitable for export_vtk.

Restricting a term to part of the mesh#

innerₕ, inner₊ and the directional products take a markers keyword, which restricts the sum to the union of the regions the labels name:

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, computed from its shape and needing no label 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 cell measures. A masked innerₕ scales like \(h\) and vanishes under refinement, where a boundary integral does not. A Neumann or Robin term needs the integral, which is inner_Γ:

β, g = 1.7, x -> 2.0 + x[1]
a_robin = form(Wd, Wd, (u, v) -> inner_Γ(β * u, v; markers = (:xmax,)))
l_neumann = form(Wd, v -> inner_Γ(g, v; markers = (:xmax,)))
sum(assemble(l_neumann))

3.0

That is \(g(x_N)\) exactly, whatever the refinement. In 1D a face is a point of measure 1, so the term is the endpoint pairing \(g(x_N)\,v(x_N)\). In 2D the same expression is an edge integral weighted by the transverse half-spacing, and in 3D a face integral. inner_Γ integrates over whole coordinate faces. A label is read from its definition: a marker in domain(...) selects the faces its points cover, so under :left => :ymin the label :left is the bottom edge. A label that is not a marker is a face name (:xmin to :zmax, or a viewpoint alias) and needs no marker; :boundary is a marker every mesh holds, the union of its faces.

A marker that exists nowhere the term reaches is an error:

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

ArgumentError: the marker :nope is not defined on the form's space. Available marker labels on this space are: :boundary, :interior, :left, :right. 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 must exist on every leaf. Write the term per component otherwise.

The outward normal#

A flux term \(\int_\Gamma \mathbf{F} \cdot \boldsymbol{\eta}\, v\) needs the outward normal at each boundary point. η destructures and indexes into its per-coordinate components like ∇ₕ:

Ω2 = mesh(domain(box((0.0, 0.0), (1.0, 1.0)),
    :xmin => :xmin, :xmax => :xmax, :ymin => :ymin, :ymax => :ymax), (6, 6), (true, true))
W2 = gridspace(Ω2)
ηₓ, ηᵧ = η
F1 = Rₕ(W2, x -> 1.0 + x[1])
F2 = Rₕ(W2, x -> 2.0 + x[2])
l_flux = form(W2, v -> inner_Γ(F1 * ηₓ + F2 * ηᵧ, v; markers = (:xmax,)))
sum(assemble(l_flux))

2.0

Point sources#

A continuous source \(f(x)\) is a density, integrated over cell volumes as \(\ell(v) = (f, v)_h\). A point source is a singular functional, the Dirac distribution \(S\,\delta(x - x_0)\):

Writing dirac(x0, strength) inside innerₕ or inner₊ evaluates it without cell-measure weighting:

l_pt = form(Wₕ, v -> innerₕ(dirac(0.5, 2.0), v))
sum(assemble(l_pt))

2.0

When \(x_0\) is a grid node, the strength \(S\) lands in that degree of freedom. When it falls inside a cell, \(S\) is spread over the \(2^D\) surrounding vertices with multilinear weights, and the total is conserved:

l_off = form(Wₕ, v -> innerₕ(dirac(0.35, 1.0), v))
sum(assemble(l_off)) ≈ 1.0

true

Several sources take a vector of coordinates and a vector of strengths. Each coordinate is a tuple (a 1-tuple in 1D), because a bare vector of reals is read as the coordinates of one point:

l_multi = form(Wₕ, v -> innerₕ(dirac([(0.2,), (0.7,)], [1.0, -1.0]), v))
sum(assemble(l_multi)) ≈ 0.0

true

A strength wrapped in a Ref updates live, without rebuilding the form:

S_live = Ref(1.0)
l_dyn_pt = form(Wₕ, v -> innerₕ(dirac(0.5, S_live), v))
b_dyn_pt = assemble(l_dyn_pt)
S_live[] = 3.5
assemble!(b_dyn_pt, l_dyn_pt)
sum(b_dyn_pt) ≈ 3.5

true

Automatic simplification#

form resolves the expression once and simplifies it before storing it. The assembler routes summands one at a time, so every top-level + is a separate sweep over the mesh. The pass is small and keeps terms as written: it does not combine like terms or factor a shared argument. Two identical terms stay two terms and are summed at assembly:

a_dup = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + innerₕ(u, v))
expression(a_dup)     # two terms, as written

"innerₕ(u, v) + innerₕ(u, v)"

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

true

Write 2 * innerₕ(u, v) when one sweep is wanted. Nested scales fold into one:

a_fold = form(Wₕ, Wₕ, (u, v) -> 2 * (3 * innerₕ(u, v)))
expression(a_fold)    # 6 * innerₕ(u, v)

"6 * innerₕ(u, v)"

A zero-scaled term leaves nothing behind, so a coefficient can switch a term off without a branch around the form:

a_off = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + 0 * inner₊(∇ₕ(u), ∇ₕ(v)))
expression(a_off)     # innerₕ(u, v): the zero term is gone

"innerₕ(u, v)"

Matrix(assemble(a_off)) ≈ Matrix(assemble(form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v))))

true

A rule whose outcomes have different node types fires only when the compiler settles it from types, so the scalar of a dropped zero term must be an Integer. A Float64 known at run time, or a Ref, leaves the term in: same numbers, one extra sweep. And the pass stops at an inner product's own arguments, so a scalar buried in a difference such as innerₕ(D₋ₓ(2 * u), v) is invisible to it. Write 2 * innerₕ(D₋ₓ(u), v) instead.

The internals page on forms has the docstrings of the stencil algebra and of simplify_ast.

Ghost-point problems#

Neumann and Robin conditions on ghost points, the relations that carry them and the solutions they give are in the ghost point tutorial. Two how-to pages build on it: time-dependent boundary data and pure Neumann problems.

Threading#

Assembly follows the execution policy of the space's backend, Serial() by default. Build the mesh with mesh(Ω, n; backend = backend(policy = Parallel())) to assemble on several threads.