GitHub

Grid spaces and discrete functions

Discrete function spaces in Bramble connect continuous mathematical functions to finite-dimensional discrete degrees of freedom defined over computational meshes.

This tutorial introduces:

  • Scalar grid spaces defined on discrete meshes
  • Multi-component composite spaces for vector and tensor fields
  • Vector elements representing discrete grid functions
  • Component indexing and zero-copy field views
  • Logical grid layouts via reshaped matrix views
  • Nodal projection and restriction ($R_h$)
  • Cell averaging operators ($\mathrm{avg}_h$)
  • Discrete inner products and norms (innerₕ, normₕ, norm₁ₕ, snorm₁ₕ)

1. Constructing scalar grid spaces

A scalar grid space represents discrete scalar fields over a mesh. To create a scalar space, call gridspace on a mesh $\Omega_h$:

using Bramble

# 1. Define a 2D domain: [0, 1] × [0, 1]
Ω = domain(box((0.0, 0.0), (1.0, 1.0)))

# 2. Discretize into a uniform 5 × 5 mesh with periodic boundary conditions
Ωₕ = mesh(Ω, (5, 5), (true, true))

# 3. Construct a scalar grid space on the mesh
Wₕ = gridspace(Ωₕ)

The resulting space Wₕ is an instance of ScalarGridSpace.

Degrees of freedom and quadrature weights

The number of degrees of freedom in Wₕ corresponds to the total number of grid points in the mesh:

ndofs(Wₕ)  # returns 25 (5 * 5)

To perform numerical integration and compute inner products, each degree of freedom has an associated quadrature weight given by the cell measure around that point:

w = weights(Wₕ)
length(w)  # 25

2. Multi-component composite spaces

Many physical problems involve vector-valued quantities such as velocities $\mathbf{u} = (u_x, u_y)$, displacement fields, or coupled state variables. In Bramble, multi-component spaces are represented by CompositeGridSpace.

Constructing vector spaces with power notation

The simplest way to create a vector grid space of dimension $D$ is using exponentiation:

# A 2-component vector space (e.g., 2D velocity space)
Vₕ = Wₕ^2

Alternatively, vector_gridspace can be constructed directly from a mesh:

Vₕ = vector_gridspace(Ωₕ, 2)

Inspecting composite spaces

ncomponents(Vₕ)  # 2
ndofs(Vₕ)        # 50 (2 * 25)
spaces(Vₕ)       # (Wₕ, Wₕ)

ndofs(Vₕ, Tuple) also works, but means something different here than it did for Wₕ above: on a ScalarGridSpace it is the grid's shape, one entry per spatial dimension (Nₓ, Nᵧ); on a CompositeGridSpace it is one entry per component instead — (25, 25) for Vₕ, not a shape. weights, by contrast, is not defined at all for a CompositeGridSpace: its components can sit on different meshes, so there is no single weight vector to hand back for the whole space — call weights on a components(Vₕ) leaf instead.

Composite spaces can also be constructed from distinct constituent spaces:

V_custom = CompositeGridSpace((Wₕ, Wₕ))

3. Vector elements and grid functions

A VectorElement represents a discrete field in a given grid space. It wraps a coefficient vector together with a reference to its parent function space.

Instantiating elements

You can instantiate uninitialized elements, elements filled with a constant, or wrap existing coefficients:

# Uninitialized vector element
uₕ = element(Wₕ)

# Element initialized to a constant value
u_zero = element(Wₕ, 0.0)
u_ones = element(Wₕ, 1.0)

Array operations and broadcasting

Because VectorElement <: AbstractVector, it supports standard vector indexing, length queries, and arithmetic:

uₕ[1] = 42.0
length(uₕ)  # 25

# Broadcasting preserves the parent space without unnecessary allocations
vₕ = element(Wₕ, 2.0)
wₕ = 3.0 .* uₕ .+ vₕ

Scaling by a continuous function

A plain Function has no meaning as a grid function on its own — f * uₕ restricts f to uₕ's own space first (Rₕ(space(uₕ), f)) and scales elementwise, giving back an ordinary VectorElement:

uₕ = Rₕ(Wₕ, x -> 1.0)
below_half = (x -> x[1] < 0.5) * uₕ   # same as (uₕ * (x -> x[1] < 0.5))

Useful for a spatial condition multiplying a grid function directly, including as a form's source: innerₕ((x -> x[1] < 0.5) * uₕ, v).

4. Component indexing and field extraction

When working with vector fields in a CompositeGridSpace, you often need to inspect or manipulate individual physical components (such as velocity in the $x$ or $y$ direction).

Functor call syntax and component views

Calling a vector element as a function uₕ(i) returns a VectorElement representing the $i$-th component:

uₕ = element(Vₕ)

# Extract component views using coordinate subscripts ("ₓ", "ᵧ", "₂")
uₓ = uₕ(1)
uᵧ = uₕ(2)

Degree-of-freedom ranges

To retrieve the degree-of-freedom index ranges occupied by components in the underlying flat vector, use component_range or component_ranges:

# Range of component 1: 1:25
rng1 = component_range(Vₕ, 1)

# All component ranges as a tuple: (1:25, 26:50)
rngs = component_ranges(Vₕ)

Zero-copy view semantics

Component extraction uses zero-copy array views into the parent degree-of-freedom vector. Modifying a component modifies the parent element in-place:

# Assign values directly to components
uₓ .= 1.5
uᵧ .= -2.0

# The parent vector reflects the updates immediately
parent(uₕ)

For scalar spaces, uₕ(1) cleanly returns uₕ itself.

Tuple destructuring

All components can be extracted simultaneously as a tuple using components:

uₓ, uᵧ = components(uₕ)
# or using numeric index subscripts:
u₁, u₂ = components(uₕ)

5. Logical grid layouts and reshaped matrix views

While degrees of freedom are stored internally as flat 1D vectors for linear algebra operations, finite difference stencils and field evaluations often require querying points in physical grid coordinates.

Direct multidimensional and Cartesian indexing

For any ScalarGridSpace element, degrees of freedom can be read and mutated directly by grid coordinates (uₕ[i, j] in 2D, uₕ[i, j, k] in 3D) or by CartesianIndex without calling reshape:

u_scal = element(Wₕ, 0.0)

# Direct 2D grid coordinate access
u_scal[2, 3] = 10.0
@assert u_scal[2, 3] == 10.0

# CartesianIndex indexing across dimensions
I = CartesianIndex(2, 3)
u_scal[I] = 20.0
@assert u_scal[2, 3] == 20.0

These indexing operations translate spatial grid coordinates directly into linear coefficient offsets via the mesh's LinearIndices with zero heap allocations and full @inbounds transparency.

Reshaped array views

When full multidimensional matrix operations (such as size(M) == (nx, ny), matrix factorizations, or external plotting) are needed, the zero-argument reshape(uₕ) returns a Base.ReshapedArray view of the flat coefficient vector matching the mesh geometry:

u_grid = reshape(u_scal)
size(u_grid)  # (5, 5)

# Access or mutate value through the reshaped view
u_grid[2, 3] = 10.0

Because reshape(u_scal) returns a Base.ReshapedArray view of the underlying vector, mutating u_grid modifies u_scal in-place with zero memory allocation.

Multi-component elements

For multi-component vector elements, reshape returns a tuple of reshaped arrays, one for each component:

mats = reshape(uₕ)
# mats is a Tuple containing (reshape(uₓ), reshape(uᵧ))

size(mats[1])  # (5, 5)
size(mats[2])  # (5, 5)
Flat 1D degree of freedom storage: uₕ.data (length 50) Component 1 (uₓ): DOFs 1:25 Component 2 (uᵧ): DOFs 26:50 uₕ(1) uₕ(2) VectorElement view (SubArray of length 25) VectorElement view (SubArray of length 25) reshape(uₓ) reshape(uᵧ) 5 × 5 ReshapedArray (uₓ[i, j]) 5 × 5 ReshapedArray (uᵧ[i, j])

6. Nodal restriction and projection

The nodal restriction operator $R_h$ evaluates a continuous function $f(x)$ at the discrete grid points of a mesh and stores the resulting values in a VectorElement.

Projecting scalar functions

# Define a continuous function of spatial coordinates x = (x₁, x₂)
f(x) = sin(2π * x[1]) * cos(2π * x[2])

# Allocate and project
u_proj = Rₕ(Wₕ, f)

# In-place projection into an existing element
Rₕ!(u_proj, f)

Projecting vector-valued functions

For multi-component spaces, $R_h$ accepts either a tuple of scalar functions or a vector function:

# Tuple of coordinate functions
fx(x) = x[1]
fy(x) = 2 * x[2]

Rₕ!(uₕ, (fx, fy))

# Or a function returning a tuple/vector
f_vel(x) = (sin(x[1]), cos(x[2]))
Rₕ!(uₕ, f_vel)

7. Numerical cell averaging

When discretizing conservation laws or finite volume formulations, quantities often represent cell averages rather than pointwise values.

The cell averaging operator $\mathrm{avg}_h$ integrates a function $f$ over each computational cell $\square_i$ around grid point $x_i$, normalized by the cell volume $|\square_i|$:

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

Nodal restriction: Rₕ(Wₕ, f) Pointwise evaluation: Rₕ f(xᵢ) = f(xᵢ) f(x) xᵢ₋½ xᵢ₊½ xᵢ f(xᵢ) Single evaluation at node xᵢ Cell averaging: avgₕ(Wₕ, f) Integral mean: avgₕ f(xᵢ) = |□ᵢ|⁻¹ ∫_{□ᵢ} f(x) dx f(x) xᵢ₋½ xᵢ₊½ xᵢ 6-point Gauss-Legendre quadrature (N_q = 6)

In Bramble, $\mathrm{avg}_h$ uses a tensor-product Gauss-Legendre rule (by default AVG_QUAD_POINTS = 6, exact for polynomials up to degree eleven):

# Compute cell-averaged element
u_avg = avgₕ(Wₕ, x -> exp(-x[1] - x[2]))

# In-place version
avgₕ!(u_avg, x -> exp(-x[1] - x[2]))

For multi-component spaces, averages can likewise be computed across components:

vₕ = avgₕ(Vₕ, (x -> 1.0, x -> 2.0 * x[1]))

8. Discrete inner products and norms

In continuous analysis, function spaces like $L^2(\Omega)$ and $H^1(\Omega)$ are equipped with inner products and norms:

\[(u, v)_{L^2} = \int_\Omega u(x) v(x) \, dx, \quad \|u\|_{L^2} = \sqrt{(u, u)_{L^2}}, \quad |u|_{H^1}^2 = \int_\Omega |\nabla u|^2 \, dx.\]

In Bramble, discrete functions in a ScalarGridSpace or CompositeGridSpace have direct discrete counterparts that weight grid values by cell measures and quadrature weights.

The discrete $L^2$ inner product and norm

The primary discrete inner product is innerₕ(uₕ, vₕ). It weights each point by its cell measure $w_i = |\square_i|$:

\[(u_h, v_h)_h = \sum_i w_i \, u_h(x_i) v_h(x_i).\]

The discrete $L^2$ norm normₕ(uₕ) is induced by innerₕ:

\[\|u_h\|_h = \sqrt{(u_h, u_h)_h}.\]

using Bramble

Ωₕ = mesh(domain(interval(0.0, 1.0)), 100, true)
Wₕ = gridspace(Ωₕ)
uₕ = Rₕ(Wₕ, sin)
vₕ = Rₕ(Wₕ, cos)

# Discrete L2 inner product and norm
innerₕ(uₕ, vₕ)
normₕ(uₕ)

# Exact norm identity: ‖uₕ‖ₕ² == (uₕ, uₕ)ₕ
normₕ(uₕ)^2 ≈ innerₕ(uₕ, uₕ)

Discrete Sobolev norms: $H^1$ seminorm and full $H^1$ norm

Bramble provides discrete $H^1$ Sobolev norms based on the backward discrete gradient $\nabla_{-h}$:

  • snorm₁ₕ(uₕ): the discrete $H^1$ seminorm $|u_h|_{1,h}$, defined as:

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

  • norm₁ₕ(uₕ): the full discrete $H^1$ norm, satisfying the Pythagorean identity:

    \[\|u_h\|_{1,h}^2 = \|u_h\|_h^2 + |u_h|_{1,h}^2.\]

# Discrete H¹ seminorm and full H¹ norm
snorm₁ₕ(uₕ)
norm₁ₕ(uₕ)

# Verification of identity
norm₁ₕ(uₕ)^2 ≈ normₕ(uₕ)^2 + snorm₁ₕ(uₕ)^2

Inner products on composite (vector) spaces

For vector-valued grid functions in a CompositeGridSpace (such as velocities or gradients), innerₕ sums the discrete inner products across all components:

\[(\mathbf{u}_h, \mathbf{v}_h)_h = \sum_{c=1}^{\mathrm{NC}} (u_{h,c}, v_{h,c})_h.\]

Vₕ = Wₕ^2
u_vec = Rₕ(Vₕ, (x -> sin(x), x -> cos(x)))

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

Staggered weights and directional inner products

Energy estimates in finite difference schemes often balance flux differences against intermediate values at cell faces. Bramble provides staggered inner products:

  • inner₊(uₕ, vₕ): inner product using staggered forward weights.
  • Coordinate-specific forms: inner₊ₓ, inner₊ᵧ, inner₊₂.

These staggered inner products form the exact algebraic pairing needed for summation by parts (discussed in detail in the Difference, jump and average operators tutorial).