Mesh tutorial

Bramble.jl provides structured, zero-allocation Cartesian and tensor-product mesh representations optimized for finite difference, finite volume, and mimetic discretization schemes.

In this tutorial, you will learn how to:

1. Construct 1D meshes (Mesh1D) and multi-dimensional tensor-product meshes (MeshnD).
2. Configure uniform and non-uniform coordinate distributions.
3. Query mesh geometric properties: coordinates, half-points (cell centers), spacings, cell measures, and $h_{\max}$.
4. Query mesh boundaries and interiors using CartesianIndices.
5. Access and evaluate boundary and region markers on meshes.
6. Perform in-place mesh refinement (iterative_refinement!) and coordinate relocation (change_points!).

1. Constructing meshes

Meshes in Bramble.jl are built on top of computational Domains. The primary entry point is the mesh function.

Every mesh also carries a linear-algebra Backend — chosen with the backend keyword below, or left to mesh's own default of backend(eltype(Ω)) — that fixes its vector/matrix types and whether threading-capable operations run serially or in parallel. See the backend tutorial for the full picture; nothing below depends on it.

1.1 One-dimensional meshes

To construct a 1D mesh with $N$ points over an interval $[a, b]$:

using Bramble

# 1. Define a domain
Ω = domain(interval(0.0, 1.0))

# 2. Build a uniform mesh with 11 grid points (step h = 0.1)
Ωₕ = mesh(Ω, 11)

By default, mesh generates a uniform grid. You can also specify non-uniform point distributions:

# Explicit non-uniform 1D mesh
Ωₕ_nonunif = mesh(Ω, 11, false)

For 1D meshes, is_uniform checks whether all cell widths are identical:

is_uniform(Ωₕ)          # true
is_uniform(Ωₕ_nonunif)   # false

1.2 Multi-dimensional tensor-product meshes

For 2D and 3D domains, Bramble.jl constructs a MeshnD as a Cartesian product of 1D submeshes. This allows $O(N_x + N_y + N_z)$ coordinate storage while providing full $O(N_x \times N_y \times N_z)$ grid traversal:

# 2D Unit square: [0, 1] × [0, 2]
Ω_2d = domain(interval(0.0, 1.0) × interval(0.0, 2.0))

# Create a 2D mesh with 10 × 20 grid points (uniform in both directions)
Ωₕ_2d = mesh(Ω_2d, (10, 20))

# Create a 2D mesh with mixed uniformity (uniform in x, non-uniform in y)
Ωₕ_mixed = mesh(Ω_2d, (10, 20), (true, false))

You can retrieve the underlying 1D submesh along any coordinate axis using functor call syntax:

x_mesh = Ωₕ_2d(1)  # 1D submesh in x-direction
y_mesh = Ωₕ_2d(2)  # 1D submesh in y-direction

1.3 Direct Cartesian set input and isotropic resolution

When custom boundary markers are not required, you can pass a CartesianProduct geometric set directly to mesh without wrapping it in domain. The default geometric markers :boundary and :interior are provisioned automatically:

# Direct discretization of intervals, products, and boxes
X = interval(0.0, 1.0) × interval(0.0, 2.0)
Ωₕ_direct = mesh(X, (10, 20))

# Default boundary and interior markers are available immediately
:boundary in keys(markers(Ωₕ_direct)) # true
:interior in keys(markers(Ωₕ_direct)) # true

In any dimension $D \ge 1$, passing a single integer npts::Int creates an isotropic grid with the same resolution across all coordinate axes:

# Isotropic 20 × 20 grid directly from the geometric set
Ωₕ_iso = mesh(X, 20)
size(Ωₕ_iso) # (20, 20)

# Isotropic resolution on a Domain
Ω = domain(X)
Ωₕ_iso_domain = mesh(Ω, 20)
size(Ωₕ_iso_domain) # (20, 20)

Both positional unif and keyword uniform accept a single boolean for isotropic uniformity (e.g. mesh(X, 20, false) or mesh(Ω, (10, 20); uniform = false)) or an NTuple{D, Bool} for per-axis control.

2. Accessing grid coordinates and metric properties

The mesh above is [0.0, 0.2, 0.6, 1.0], deliberately non-uniform. Four conventions are worth reading off it, because they are the ones that most often surprise:

• half_points has N + 1 entries, not N. They are the cell interfaces, and the first and last coincide with x₁ and x_N rather than being extrapolated outside the domain. Here they are [0.0, 0.1, 0.4, 0.8, 1.0].
• The cell around xᵢ spans half_points[i] .. half_points[i+1], and its width is exactly half_spacing(Ωₕ, i), which is what cell_measure(Ωₕ, i) returns. The four cells here measure [0.1, 0.3, 0.4, 0.2] and sum to the domain length.
• Boundary cells are half-width. x₁ and x_N sit on the edge of their own cell, not at its centre, which is why the first and last measures are the smallest.
• spacing looks backward and forward_spacing looks forward, so spacing(Ωₕ, i) = xᵢ − xᵢ₋₁ and forward_spacing(Ωₕ, i) = xᵢ₊₁ − xᵢ. Each has one special case at the boundary where the neighbour is missing: spacing(Ωₕ, 1) returns x₂ − x₁ and forward_spacing(Ωₕ, N) returns x_N − x_{N−1}.

2.1 Points and coordinates

• points(Ωₕ): Returns the coordinate vector (1D) or tuple of coordinate vectors (nD).
• point(Ωₕ, idx) or direct indexing Ωₕ[idx]: Evaluates the coordinate at linear index i, coordinate tuple (i, j), or CartesianIndex(i, j).

# 1D mesh point access
p3 = Ωₕ[3]          # Coordinate x₃
p3_alt = point(Ωₕ, 3)

# 2D mesh point access
p_ij = Ωₕ[2, 5]     # Coordinate tuple (x₂, y₅)

2.2 Half-points and cell centers

Finite volume and staggered-grid methods frequently require cell midpoints $x_{i+1/2}$:

# Pre-computed cell centers
hp = half_points(Ωₕ)
hp_i = half_point(Ωₕ, 3)  # x_{3+1/2}

2.3 Spacings and cell measures

A 1D mesh stores its backward spacings rather than recomputing them, so spacings(Ωₕ) hands back the whole vector and spacing(Ωₕ, i) is a single array read. forward_spacing(Ωₕ, i) reads the same vector one entry along, since $x_{i+1} - x_i$ is the backward spacing at $i+1$. The cache is rebuilt by set_points!, and so by iterative_refinement! and change_points! as well, meaning it always matches the current points.

Ωₕ = mesh(domain(interval(0.0, 1.0)), 5, false)

spacings(Ωₕ)                       # every hᵢ at once
spacings(Ωₕ)[3] == spacing(Ωₕ, 3)  # true, the accessor just indexes it

This matters for the difference operators, which need one spacing per grid point: reading the cached vector directly, rather than boxing the spacing accessor as a generic callable passed into the inner stencil loop, keeps that loop allocation-free.

# Maximum grid stepsize
h = hₘₐₓ(Ωₕ_2d)

# Control volume measure at cell (3, 4)
vol = cell_measure(Ωₕ_2d, (3, 4))

An n-dimensional mesh is a tensor product of 1D meshes, and every quantity above is built the same way. The cell around (xᵢ, yⱼ) is the rectangle spanned by the two per-axis intervals, so its measure is the product of the per-axis widths:

cell_measure(Ωₕ, CartesianIndex(3, 2))          # 0.2
half_spacing(Ωₕ(1), 3) * half_spacing(Ωₕ(2), 2)  # 0.2 — the same number

Ωₕ(k) is the 1D submesh along axis k, so anything documented for a 1D mesh applies to it directly. As in one dimension the cells tile the domain exactly — here the twelve cell measures sum to the area 1.0 — and the cells touching a boundary are correspondingly thinner along that axis.

3. Boundary and interior indexing

Bramble.jl uses Julia's native CartesianIndices for zero-overhead, multi-dimensional grid navigation:

# Complete Cartesian grid indices
idxs = indices(Ωₕ_2d)  # CartesianIndices((1:10, 1:20))

# Interior indices (excluding all boundaries)
interior = interior_indices(Ωₕ_2d)  # CartesianIndices((2:9, 2:19))

# Boundary facets as a tuple of CartesianIndices
facets = boundary_indices(Ωₕ_2d)

# Test whether an index lies on the domain boundary
is_boundary = is_boundary_index(Ωₕ_2d, CartesianIndex(1, 5))  # true

4. Markers on meshes

When creating a mesh from a labeled Domain, markers are projected onto the grid points as highly efficient BitVectors:

# Domain with boundary and obstacle markers
I = interval(0.0, 1.0)
Ω = domain(I × I,
           :left_inlet => :xmin,      # or legacy alias :left
           :right_outlet => :xmax,    # or legacy alias :right
           :walls => (:ymin, :ymax),  # or legacy alias (:top, :bottom)
           :obstacle => x -> (x[1]-0.5)^2 + (x[2]-0.5)^2 < 0.15^2)

# Generate mesh
Ωₕ = mesh(Ω, (20, 20))

# Query markers
m_dict = markers(Ωₕ)

# Retrieve bit-vector for a specific label
is_wall = index_in_marker(Ωₕ, :walls)
is_obs  = index_in_marker(Ωₕ, :obstacle)

5. Mesh adaptation and modification

Meshes in Bramble.jl are mutable structures designed for adaptive algorithms:

5.1 In-place mesh refinement

Halves every cell by inserting new points at each cell midpoint, simultaneously updating indices and reapplying domain markers:

# Refine mesh in-place
iterative_refinement!(Ωₕ)

# Point count increases: (2N_x - 1) × (2N_y - 1)
npoints(Ωₕ, Tuple)  # (39, 39)

5.2 Relocating mesh coordinates

For moving-boundary problems or non-uniform smoothing:

# Supply new coordinates for a 1D mesh (matching the 11 points of Ωₕ)
new_pts = range(0.0, 1.0, length=npoints(Ωₕ)) |> collect
change_points!(Ωₕ, new_pts)

# Or update points and re-evaluate markers for a multi-dimensional mesh
nx, ny = npoints(Ωₕ, Tuple)
new_x_pts = range(0.0, 1.0, length=nx) |> collect
new_y_pts = range(0.0, 1.0, length=ny) |> collect
change_points!(Ωₕ, markers(Ω), (new_x_pts, new_y_pts))