Grid spaces and discrete functions#

You will put a function on a mesh, measure it, and stack copies into a composite space with vector elements. The mesh tutorial builds the mesh and the geometry tutorial the domain it sits on. Difference operators act on the elements built here.

A grid space is the discrete function space over a mesh. It fixes how many degrees of freedom a field has and which quadrature weight each one carries. A VectorElement is a function in that space.

The running problem is small: sample $f(x) = \sin(\pi x)$ on the unit interval and measure it, first as one field, then as a pair of fields.

A grid space over a mesh#

gridspace builds a grid space over a mesh. ndofs counts its degrees of freedom, and with Tuple it gives the grid shape:

using Bramble
import Bramble: component_range, component_ranges, weights

Ω = domain(interval(0.0, 1.0))
Ωₕ = mesh(Ω, 11)
Wₕ = gridspace(Ωₕ)

ndofs(Wₕ), ndofs(Wₕ, Tuple)

(11, (11,))

There is one degree of freedom per mesh point, and points(Wₕ) is points(Ωₕ).

Discrete functions#

Rₕ restricts a function to the grid: it evaluates it at every point and returns a VectorElement, the discrete function.

uₕ = Rₕ(Wₕ, x -> sin(π * x))

length(uₕ), uₕ[6]

(11, 1.0)

VectorElement <: AbstractVector, so indexing, length and broadcasting work on it, and broadcasting keeps the parent space:

wₕ = 3.0 .* uₕ .+ element(Wₕ, 2.0)

space(wₕ) === Wₕ, wₕ[6]

(true, 5.0)

element allocates an uninitialised field; element(Wₕ, α) fills it with α and element(Wₕ, v) copies a vector. Rₕ! writes into an element that already exists.

Restriction and cell averaging#

Nodal restriction $$R_h$$ evaluates a function at the grid points; cell averaging $\mathrm{avg}_h$ integrates it over the cell around each point and divides by the cell measure:

\mathrm{avg}_h f(x_i) = \frac{1}{|\square_i|} \int_{\square_i} f(x) \, dx

Rₕ allocates and Rₕ! writes in place. avgₕ and avgₕ! mirror them, using a tensor-product Gauss-Legendre rule (AVG_QUAD_POINTS = 6, exact through degree eleven):

u_avg = avgₕ(Wₕ, x -> exp(-x))
avgₕ!(u_avg, x -> exp(-x))

u_avg[1]

0.9754115099857197

For a rough source term, avgₕ averages out the roughness that Rₕ samples.

A plain Function is not a grid function, so multiplying by one restricts it first (Rₕ(space(uₕ), f)) and scales pointwise. This lets a spatial condition multiply a field directly, including inside a form:

oneₕ = Rₕ(Wₕ, x -> 1.0)
below_half = (x -> x[1] < 0.5) * oneₕ

sum(below_half)

5.0

Inner products and norms#

Each degree of freedom carries a quadrature weight, the cell measure around its point. innerₕ weights each point by it (uₕ is conjugated for complex data, and markers = (...) restricts the sum), and normₕ is the norm it induces:

(u_h, v_h)_h = \sum_i |\square_i| \, u_h(x_i) v_h(x_i), \qquad \|u_h\|_h = \sqrt{(u_h, u_h)_h}.

vₕ = Rₕ(Wₕ, x -> cos(π * x))

innerₕ(uₕ, vₕ), normₕ(uₕ), normₕ(uₕ)^2 ≈ innerₕ(uₕ, uₕ)

(-5.148621912560444e-18, 0.7071067811865476, true)

The norm of $\sin(\pi x)$ on $$[0, 1]$$ is $1/\sqrt{2} \approx 0.707$ , which normₕ reproduces here. The discrete Sobolev norms are built on the backward gradient ∇ₕ:

|u_h|_{1,h}^2 = \sum_{d=1}^D \|D_{-x_d} u_h\|_+^2, \qquad \|u_h\|_{1,h}^2 = \|u_h\|_h^2 + |u_h|_{1,h}^2.

snorm₁ₕ(uₕ), norm₁ₕ(uₕ), norm₁ₕ(uₕ)^2 ≈ normₕ(uₕ)^2 + snorm₁ₕ(uₕ)^2

(2.212317420824743, 2.322573652327229, true)

inner₊ is the staggered counterpart of innerₕ: it weights by the backward spacing along the chosen direction (zero at the first node) and by the cell width across the others. A trailing :x/:y/:z (or 1/2/3) argument picks a single direction. That pairing makes summation by parts exact; the operators tutorial derives it.

D = ∇ₕ(uₕ)

inner₊(D, D, :x) ≈ snorm₁ₕ(uₕ)^2

true

Try this

Build the mesh as mesh(Ω, 11, false), a non-uniform one with the same number of points, and rerun the blocks above. ndofs is unchanged, but the interior weights are no longer equal (on a uniform mesh only the two end weights are halved), so normₕ and innerₕ change while the weights still sum to the interval length. The mesh tutorial builds such meshes on purpose.

Composite spaces and vector elements#

A composite space stacks copies of a space for a vector field or a coupled system. Its elements are VectorElements too: one flat vector holding every component, one after another. Wₕ^Val(N) stacks N copies:

Vₕ = Wₕ^Val(2)

ncomponents(Vₕ), ndofs(Vₕ), ndofs(Vₕ, Tuple)

(2, 22, (11, 11))

ndofs(Vₕ, Tuple) has one entry per component (on a scalar space it is the grid shape). A component is one of the stacked fields, and calling a composite element with i returns the i-th component as a view:

uvec = element(Vₕ)
uₓ = uvec(1)
uᵧ = uvec(2)

uₓ .= 1.5
uᵧ .= -2.0

parent(uvec)[[1, 12]]

2-element Vector{Float64}:
  1.5
 -2.0

A component is a view; writing to it writes to the parent. For a scalar space uₕ(1) is uₕ itself. components destructures all of them:

c₁, c₂ = components(uvec)

c₁ === uₓ, c₂ === uᵧ

(true, true)

On a composite space Rₕ! takes a tuple of scalar functions, or one function returning a tuple, and innerₕ sums over components:

Rₕ!(uvec, (x -> sin(π * x), x -> cos(π * x)))

normₕ(uvec)^2 ≈ normₕ(uvec(1))^2 + normₕ(uvec(2))^2

true

The storage is flat, and reshape gives each component a grid-shaped view of the same memory. For a 2D space on a 5 × 5 mesh the layout is:

Try this

Replace Val(2) with Val(3) and add a third function to the tuple. ncomponents and ndofs follow.

Ghost spaces#

A grid space can also carry the ghost points of a mesh beyond named facets, through the ghosts keyword of gridspace. The ghost point tutorial covers their degrees of freedom, weights and time loops.