Mesh tutorial#

What you will learn. How a domain becomes a mesh, how to make the points non-uniform, what the spacings, half points and cell measures that difference operators read mean.

What you need first. The geometry tutorial, for sets, markers and domains.

Where next. The space tutorial puts discrete functions on a mesh.

A first mesh#

Take the unit interval and ask for eleven points. mesh takes a domain and a point count:

using Bramble
import Bramble: cell_measure, change_points!, half_point, half_points, half_spacing,
    index_in_marker, indices, interior_indices, is_boundary_index, is_uniform, point,
    spacings

Ω = domain(interval(0.0, 1.0))
Ωₕ = mesh(Ω, 11)          # 11 uniformly spaced points

points(Ωₕ)

11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.30000000000000004
 0.4
 0.5
 0.6000000000000001
 0.7000000000000001
 0.8
 0.9
 1.0

The points are equally spaced, \(h = 0.1\), and the largest spacing confirms it:

hₘₐₓ(Ωₕ), is_uniform(Ωₕ)

(0.10000000000000009, true)

A geometric set can be passed directly when no custom labels are needed. :boundary and :interior are provisioned either way:

Ωₕ_direct = mesh(interval(0.0, 1.0) × interval(0.0, 2.0), 20)      # isotropic: 20 × 20

:boundary in keys(markers(Ωₕ_direct)), :interior in keys(markers(Ωₕ_direct))

(true, true)

The mesh carries more than coordinates. Spacings, half points and cell measures, which the difference operators read, are computed from the points and cached.

Non-uniform meshes#

Boundary layers and corners need points where the solution changes fast. The third positional argument of mesh, or the uniform keyword, chooses the distribution per axis. false draws interior points at random and sorts them, which the convergence studies in this manual use to expose order reduction on non-uniform grids:

Ωₕ_rand = mesh(Ω, 11, false)

is_uniform(Ωₕ), is_uniform(Ωₕ_rand)

(true, false)

To place the points yourself, build any mesh and replace its coordinates. change_points! keeps the point count, so it can grade a mesh toward the left end:

Ωₕ_graded = mesh(Ω, 11)
change_points!(Ωₕ_graded, collect(range(0.0, 1.0, length = npoints(Ωₕ_graded)) .^ 2))

points(Ωₕ_graded), hₘₐₓ(Ωₕ_graded)

([0.0, 0.010000000000000002, 0.04000000000000001, 0.09, 0.16000000000000003, 0.25, 0.36, 0.48999999999999994, 0.6400000000000001, 0.81, 1.0], 0.18999999999999995)

Squaring the uniform points packs them near \(x = 0\) and stretches the cells near \(x = 1\), so hₘₐₓ is no longer \(0.1\). Bramble.set_points! is the variant that also accepts a different point count: it rebuilds the :boundary and :interior markers for the new grid, and refuses to change the point count of a mesh carrying custom labels. It also refuses any move of a mesh carrying a label defined by a predicate of position, and of a submesh of one, since those sets would stay as the old points gave them; call change_points!(Ωₕ, markers(Ω), pts) on the whole mesh instead.

A grid space built before set_points!, change_points! or iterative_refinement! is stale: it keeps the weights of the old points, and its inner products and norms throw until gridspace(Ωₕ) is called again.

Half points, spacings and cell measures#

One small non-uniform mesh exposes every convention at once. Four points at \(0, 0.2, 0.6, 1\):

Ωₕ_fig = mesh(domain(interval(0.0, 1.0)), 4, true)
Bramble.set_points!(Ωₕ_fig, [0.0, 0.2, 0.6, 1.0])

points(Ωₕ_fig)

4-element Vector{Float64}:
 0.0
 0.2
 0.6
 1.0

Conventions:

• half_points has N + 1 entries. They are the cell interfaces, and the first and last coincide with x₁ and x_N.

• The cell around xᵢ spans half_points[i] to half_points[i+1], with width half_spacing(Ωₕ, i), which is what cell_measure(Ωₕ, i) returns.

• Boundary cells are half-width, since x₁ and x_N sit on the edge of their cell.

• spacing looks backward, forward_spacing forward: spacing(Ωₕ, i) = xᵢ - xᵢ₋₁ and forward_spacing(Ωₕ, i) = xᵢ₊₁ - xᵢ, each falling back to its neighbour at the end where the stencil runs out.

half_points(Ωₕ_fig), [cell_measure(Ωₕ_fig, i) for i in 1:4]

([0.0, 0.1, 0.4, 0.8, 1.0], [0.1, 0.3, 0.4, 0.2])

The four cells sum to the domain length.

Meshes in two dimensions#

Above one dimension a mesh is a tensor product of 1D submeshes, stored per axis, so coordinate storage is \(O(N_x + N_y + N_z)\) while the grid it addresses is the full product. The point count is one integer per axis, or a single integer for the same resolution along every axis. The distribution flag is per axis too: (true, false) would be uniform in \(x\) and random in \(y\).

Ωₕ_2d = mesh(domain(interval(0.0, 1.0) × interval(0.0, 2.0)), (10, 20), (true, true))

npoints(Ωₕ_2d, Tuple), size(Ωₕ_2d)

((10, 20), (10, 20))

Ωₕ_2d(k) is the submesh along axis k, and everything said above about a 1D mesh applies to it:

Ωₕ_2d(1)

Mesh1D {1D, Float64}
10 points
Domain: [0.0, 1.0]
Spacing: uniform • hₘₐₓ = 0.111111
Markers: 2 labels • :boundary (2), :interior (8)

The cell around \((x_i, y_j)\) is the rectangle spanned by the two per-axis intervals. The picture shows a \(4 \times 3\) mesh with non-uniform points on both axes:

The cell around (x₃, y₂) is highlighted. Its measure is the product of the two widths:

Ω_prod = domain(interval(0.0, 1.0) × interval(0.0, 1.0))
Ωₕ_prod = mesh(Ω_prod, (4, 3), (true, true))
change_points!(Ωₕ_prod, markers(Ω_prod), ([0.0, 0.2, 0.6, 1.0], [0.0, 0.5, 1.0]))

cell_measure(Ωₕ_prod, CartesianIndex(3, 2)),
half_spacing(Ωₕ_prod(1), 3) * half_spacing(Ωₕ_prod(2), 2)

(0.2, 0.2)

The cells tile the domain exactly: the twelve measures sum to the area.

sum(cell_measure(Ωₕ_prod, I) for I in indices(Ωₕ_prod))

0.9999999999999999

Above 1D, pts is a tuple of per-axis coordinate vectors. In every dimension, change_points!(Ωₕ, markers(Ω), pts) re-evaluates the labels at the new positions, while the two-argument form leaves them as they are.

Boundary, interior and markers#

Indices are Julia's own CartesianIndices:

indices(Ωₕ_prod), interior_indices(Ωₕ_prod)

(CartesianIndices((4, 3)), CartesianIndices((2:3, 2:2)))

is_boundary_index(Ωₕ_prod, CartesianIndex(1, 2)),
is_boundary_index(Ωₕ_prod, CartesianIndex(2, 2))

(true, false)

Building a mesh over a labelled domain projects each marker onto the grid points as a BitVector, so membership is one lookup. Treat them as read-only: assembly does not see an edit made to one in place. Here the domain has an inlet, an outlet, walls, and a disc-shaped obstacle given by a predicate:

I = interval(0.0, 1.0)
Ω_marked = domain(I × I,
    :left_inlet => :xmin,
    :right_outlet => :xmax,
    :walls => (:ymin, :ymax),
    :obstacle => x -> (x[1] - 0.5)^2 + (x[2] - 0.5)^2 < 0.15^2)

Ωₕ_marked = mesh(Ω_marked, (20, 20))

sum(index_in_marker(Ωₕ_marked, :walls)), sum(index_in_marker(Ωₕ_marked, :obstacle))

(40, 24)

Ghost points#

A mesh padded with one point beyond each boundary facet, for Neumann and Robin conditions, is the subject of the ghost point tutorial.

Refinement#

iterative_refinement! inserts a point at every cell midpoint and updates the indices. Refinement is uniform and dyadic: N points become 2N - 1 along every axis. Above 1D, refine the whole mesh: iterative_refinement!(Ωₕ(1)) on one submesh leaves the mesh's index set and markers sized for the old grid, and gridspace refuses it.

iterative_refinement!(Ωₕ) rebuilds :boundary and :interior and throws an ArgumentError on a mesh carrying any other label. The two-argument form re-evaluates every label on the new points. A grid space built before the call is stale, as after change_points!.

Ωₕ_ref = mesh(Ω_marked, (20, 20))
iterative_refinement!(Ωₕ_ref, markers(Ω_marked))

npoints(Ωₕ_ref, Tuple), sum(index_in_marker(Ωₕ_ref, :obstacle))

((39, 39), 101)