Spaces#

Coordinate axis labels used in documentation and error messages.

Corresponds to the x, y, and z spatial dimensions (1st, 2nd, and 3rd dimensions).

See also: _BRAMBLE_var2symbol

Subscript Unicode symbols for x, y, z coordinates used in operator notation.

These symbols are used to generate directional operator aliases via metaprogramming.

Examples

• D₊ₓ - forward difference in x-direction

• Mᵧ - backward average in y-direction

• jump₂ - jump in the z-direction

See also: _BRAMBLE_var2label

AbstractSpaceType{N}

Abstract supertype for all function spaces defined on a mesh.

This is the top-level abstraction for a grid-based function space. The parameter N represents the number of components of the field (e.g., N=1 for a scalar field, N=3 for a 3D vector field).

Innerh

Selector for the standard discrete \(L^2\) inner product weighted by cell measures.

The weights are the volumes (lengths in 1D, areas in 2D, volumes in 3D) of grid cells, denoted \(|\square_k|\). This is the most common inner product for finite difference methods and corresponds to the trapezoid rule for integration on non-uniform grids.

Mathematical form

For a 2D grid:

where \(|\square_{i,j}|\) is the area of the cell centered at \((x_i, y_j)\).

Examples

# Compute L² inner product
result = innerₕ(uₕ, vₕ)  # Uses Innerh() internally

# Compute L² norm
norm_value = normₕ(uₕ)  # Equivalent to sqrt(innerₕ(uₕ, uₕ))

See also: Innerplus, innerₕ, normₕ

Innerplus

Selector for modified discrete \(L^2\) inner products using staggered grid spacings.

These inner products use a combination of backward spacings \(h_i = x_i - x_{i-1}\) (zero at the first node) and cell widths \(h_{i+1/2}\), appearing naturally in energy estimates for finite difference operators. Different spatial directions may have different weight formulas.

The modified inner products are used for:

• Proving discrete energy stability

• Analyzing discrete conservation properties

• Constructing stable finite difference schemes

Mathematical form

For a 2D grid in the x-direction:

Examples

using Bramble: inner₊ₓ, inner₊ᵧ
# These functions use Innerplus internally
result_x = inner₊ₓ(uₕ, vₕ)  # Modified inner product, x-direction
result_y = inner₊ᵧ(uₕ, vₕ)  # Modified inner product, y-direction

See also: Innerh, inner₊ₓ, weights

mesh_type(Wₕ::AbstractSpaceType) -> Type{<:AbstractMeshType}
mesh_type(::Type{<:AbstractSpaceType}) -> Type{<:AbstractMeshType}

Returns the type of the mesh associated with the function space Wₕ. Also works if the argument is the type of the space.

SeparableWeights{D, T, VT}(factors::NTuple{D, VT}, dims::NTuple{D, Int})

A lazy, separable weight vector over a D-dimensional grid: entry I (linear or CartesianIndex{D}) is \(\prod_{d=1}^D\) factors[d][I[d]], computed on every access rather than stored once for the whole grid.

Every family SpaceWeights offers is one of these (gpena/Bramble.jl#115, #234; O(nD) storage rather than O(n^D)): innerh and each entry of innerplus are built once, when the space is constructed, from the per-axis aligned and cellfactor factors, and stored on SpaceWeights – so reading them costs no more than the field access plus whatever the caller's own indexing does. weights(Wₕ, Val(S)) for |S| ≥ 2 builds a fresh one on every call instead, from those same two tuples, since nothing yet asks for the same larger set twice in a hot loop.

A linear getindex converts to a CartesianIndex first (one division per axis); a caller that already holds the CartesianIndex – an assembly loop over local_stencil, for instance – should use it directly and skip that cost. src/space/inner_product.jl's _dot/_dot_masked do this for every weight family alike, and so does src/operators/inner.jl's compute_weight for InnerH, InnerPlus{Dim} and InnerPlusSet (|S| ≥ 2): a symbolic innerₕ/inner₊ₓ/etc. term inside a form reads its weight by the CartesianIndex the assembly loop already holds, with no division per point (gpena/Bramble.jl#428).

On device-backed factors (gpena/Bramble.jl#310), getindex raises an error naming host_weights rather than letting GPUArraysCore throw its own scalar-indexing message: Array(w) and host_weights(w) are the two bulk-transfer escapes, each pulling the D per-axis factor vectors to the host once and evaluating the tensor product (or keeping it lazy, for host_weights) there instead.

Fields

• factors::NTuple{D, VT}: the per-axis vectors multiplied together at each index.

• dims::NTuple{D, Int}: the grid shape, npoints(Ωₕ, Tuple).

See also: weights, SpaceWeights.

SpaceWeights(innerh::SeparableWeights, innerplus::NTuple{D}, aligned::NTuple{D, VT}, cellfactor::NTuple{D, VT}, built_version::Int)
SpaceWeights{D, T, VT}(innerh, innerplus, aligned::NTuple{D, VT}, cellfactor::NTuple{D, VT}, built_version::Int)

Holds the diagonal weight vectors for a grid space's discrete inner products, both the standard \(L^2\) weights and the staggered ones, precomputed once rather than recomputed on every call.

A mutable struct whose fields are all const: built once by space_weights and never changed, it is held by reference, so a ScalarGridSpace stores a pointer to it rather than the weights inline. Every copy of a space (a form's, an operator's result element, an assembly walk's rebuild) carries that pointer, not the weights, which keeps the 3D form tree under Enzyme's type-analysis limit (gpena/Bramble.jl#1154). Two spaces are === only when they share one mesh and one SpaceWeights object.

Fields

• innerh::SeparableWeights{D, T, VT}: weight for the standard discrete \(L^2\) inner product (:innerₕ), based on cell measures (\(|\square_k|\)). Built once, straight from cellfactor (no axis is in S = ()), and stored here rather than recomputed on every access – it is a lazy per-axis product rather than a full-grid vector (gpena/Bramble.jl#115, #234), so weights(Wₕ, Innerh()) costs the field read plus whatever the caller's own indexing does; see SeparableWeights.

• innerplus::NTuple{D, SeparableWeights{D, T, VT}}: one such lazy weight per spatial direction, for the modified, staggered inner products (:inner₊ₓ, :inner₊ᵧ, etc.).

• aligned::NTuple{D, VT}: per-axis factor \(h_d(i)\) (length npoints(Ωₕ, d), not the full grid) – the same values innerplus uses on the axis aligned with its own difference direction. Kept so weights(Wₕ, Val(S)) can build any staggered set from O(D) numbers instead of a fresh O(n^D) vector (gpena/Bramble.jl#115, #234).

• cellfactor::NTuple{D, VT}: per-axis factor \(h_d(i+1/2)\) (also length npoints(Ωₕ, d)). This is innerh's own per-axis cell measure – innerh's factors tuple is cellfactor, not a copy of it – and doubles as every staggered direction's transverse factor: the two coincide for any axis with more than one point (both read the submesh's own cached half-spacings), so one vector serves both roles rather than two.

• built_version::Int: the mesh's _mesh_version at the moment these weights were computed (gpena/Bramble.jl#221) – weights re-checks it against the mesh's current version on every access, so a space built before an in-place mutation (set_points!, change_points!, iterative_refinement!) throws naming the mismatch rather than silently returning weights for a mesh that no longer exists. A form's assembly walk checks it once, where the walk starts, not at every point.

For a detailed explanation of the mathematical formulas corresponding to these weights, please refer to the documentation for ScalarGridSpace.

eltype(Wₕ::ScalarGridSpace) -> Type
eltype(::Type{<:ScalarGridSpace}) -> Type

Returns the element type of vectors in this space (e.g., Float64).

See also: backend

_innerh_weights!(u, Ωₕ::AbstractMeshType)

Builds the weights for the standard discrete \(L^2\) inner product, \(inner_h(\cdot, \cdot)\), on the space of grid functions, following the order of the points provided by indices(Ωₕ). The values are stored in vector u.

_innerplus_weights!(u::VT, Ωₕ, component = 1) where VT

Builds a set of weights based on the spacings, associated with the component-th direction, for the modified discrete \(L^2\) inner product on the space of grid functions, following the order of the points provided by indices(Ωₕ). The values are stored in vector u.

host_weights(Wₕ::ScalarGridSpace) -> ScalarGridSpace
host_weights(w::SeparableWeights) -> SeparableWeights

Return a host-resident mirror of Wₕ/w, one bulk transfer per underlying array rather than one scalar read per grid point (gpena/Bramble.jl#94, #310).

On a ScalarGridSpace, mesh and weights alike come along: allocate_system_matrix's sparsity-pattern walk evaluates a form's AST through local_stencil and so reads both this space's weights and its mesh's spacings one grid point at a time; see the comment above _host_mirror_mesh for why both, not only the weight vectors the name alone might suggest. On a bare SeparableWeights (gpena/Bramble.jl#310), only its D per-axis factor vectors move – map(Array, w.factors) – since that is everything a SeparableWeights owns; Array(w) calls this first and then forms the tensor product on the host.

A no-op on a space or weight that already lives on the host: host_weights(Wc) === Wc / host_weights(w) === w, not a copy, the same guarantee host_spacings gives on a CPU mesh (a bug host_spacings once had – Array(v) === v is false for a Vector, so a copy-on-host would silently repeat it here).

See also: host_spacings, weights, SeparableWeights.

^(Wₕ::ScalarGridSpace, ::Val{N}) where N -> CompositeGridSpace{N}
^(Wₕ::ScalarGridSpace, N::Int) -> CompositeGridSpace{N}

Constructs an N-component vector grid space from a scalar grid space Wₕ using mathematical exponentiation syntax: Vₕ = Wₕ^2 or Vₕ = Wₕ^dim(mesh).

Wₕ^1 is Wₕ, for both the Int and Val spellings.

The Int spelling only accepts 1 <= N <= 3, each branch-unswitched to a literal Wₕ^Val(N) call: Julia's return-type inference only union-splits up to 3 concrete types, so this is the largest range that stays inferrable even when N is a runtime value the compiler cannot constant-fold (e.g. threaded through a generic function argument). It also covers every mesh dimension this package supports (Wₕ^dim(mesh)). For N > 3, call Wₕ^Val(N) directly, which has no upper bound.

A ghost space Cₕ (a GhostGridSpace) takes Cₕ^N the same way, giving a composite of ghost spaces; its methods live with the ghost composites (assembly/ghost_composite.jl).

getindex(Wₕ::CompositeGridSpace, i::Int) -> AbstractSpaceType

Returns the i-th top-level sub-space of Wₕ. It is itself a composite when Wₕ is nested. leaf_count counts the scalar leaves across all levels.

length(Wₕ::CompositeGridSpace) -> Int

Returns the number of top-level sub-spaces of Wₕ. A nested composite counts its direct children; leaf_count counts the scalar leaves.

_shares_one_mesh(comps::Tuple) -> Bool

Whether every leaf element in comps (a composite's components(uₕ)) sits on the exact same mesh object.

A composite's leaves need not share a mesh: heterogeneous composites, whose leaves are built over differently-sized meshes, are a supported pattern (see the interpolation tutorial). When they do share one, a single evaluation of a vector-valued function at one leaf's grid points is valid for every leaf, which is what the scatter paths of Rₕ!/avgₕ! use this to decide; when they do not, each leaf needs its own evaluation, since there is no shared "grid point i" across differently-sized meshes.

leaf_count(Wₕ::AbstractSpaceType) -> Int

Returns the number of leaf scalar spaces in Wₕ. For a scalar space, returns 1. For a composite space, returns the total flattened number of scalar leaves.

leaf_spaces_offsets(Wₕ) -> Tuple

The scalar spaces underneath Wₕ paired with their offsets into the global degree of freedom vector, depth first and left to right, as a tuple of (space, offset).

A scalar space is its own only leaf, at offset zero.

copyto!(dest::VectorElement, src::VectorElement) -> VectorElement

Copies the coefficients of src into dest in place, and returns dest.

dest and src must share one grid space (===); the size check alone is not enough.

fill!(uₕ::VectorElement, α) -> VectorElement

Fills VectorElement uₕ with the value α and returns uₕ.

getindex(uₕ::VectorElement{<:ScalarGridSpace}, I::CartesianIndex)

Access field degrees of freedom by grid coordinates through a CartesianIndex.

Translates spatial grid coordinates directly into flat linear coefficient offsets using the mesh's LinearIndices with zero heap allocations and full @inbounds transparency.

VectorElement is an AbstractVector, so an Integer multi-index does not address grid coordinates: it follows Base's own trailing-index rule for an AbstractArray with fewer dimensions than indices given (uₕ[k, 1] == parent(uₕ)[k], and any trailing index other than 1 throws a BoundsError). Grid-coordinate access always goes through CartesianIndex instead, which is what generic AbstractVector code (e.g. SparseArrays, which reads u[k, 1] while building a sparse column) needs uₕ[i, j] to not mean.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (10, 10))
Wₕ = gridspace(Ωₕ)
uₕ = element(Wₕ, 0.0)

# Set and get via CartesianIndex grid coordinates
I = CartesianIndex(2, 3)
uₕ[I] = 42.0
uₕ[I] == 42.0

# A bare Integer multi-index is linear (trailing-index) indexing, not grid coordinates
uₕ[5, 1] == parent(uₕ)[5]

See also: VectorElement, ScalarGridSpace, reshape

reshape(uₕ::VectorElement)

Reshapes the flat coefficient vector of uₕ into a multidimensional array that matches the logical layout of the grid points.

• For a scalar space, this returns a D-dimensional array.

• For an N-component vector space, it returns an N-tuple of arrays, one for each component.

This zero-argument form is specific to VectorElement and does not conflict with reshape(A, dims); it does shadow Base.reshape(A) = reshape(A, ())'s 0-dimensional result for this type specifically, in favor of the shape a grid function actually has.

setindex!(uₕ::VectorElement{<:ScalarGridSpace}, val, I::CartesianIndex) -> VectorElement

Mutate field degrees of freedom by grid coordinates through a CartesianIndex, in-place.

Translates spatial grid coordinates directly into flat linear coefficient offsets using the mesh's LinearIndices with zero heap allocations and full @inbounds transparency.

VectorElement is an AbstractVector, so an Integer multi-index does not address grid coordinates: it follows Base's own trailing-index rule for an AbstractArray with fewer dimensions than indices given (uₕ[k, 1] = v writes parent(uₕ)[k], and any trailing index other than 1 throws a BoundsError). Grid-coordinate access always goes through CartesianIndex instead, which is what generic AbstractVector code (e.g. SparseArrays, which reads and writes u[k, 1] while building a sparse column) needs uₕ[i, j] = val to not mean.

Returns uₕ matching Base collection conventions.

_batch_broadcast!(v, bc, ax) -> Nothing

CpuPolyester's _polyester_broadcast!, filled by BramblePolyesterExt (one _broadcast_band! per band under Polyester.@batch). Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_find_vec_in_broadcast(bc)

Internal helper to extract a VectorElement from a broadcast expression.

Recursively searches through the arguments of a broadcast expression tree to find a VectorElement instance. This is used by the broadcasting machinery to determine which function space should be used for the result.

Arguments

• bc: A broadcast expression, tuple of arguments, or individual value

Returns

• The first VectorElement found in the expression tree

• nothing if no VectorElement is found

Implementation notes

Uses multiple dispatch to handle:

• Broadcasted objects: Extract and search arguments

• Tuples: Recursively search each element

• VectorElement: Return immediately (found!)

• Other types: Return nothing and continue searching

This enables broadcasts like uₕ .+ vₕ .* 2 to automatically preserve the space information.

space_type(::Type{<:VectorElement}) -> Type{<:AbstractSpaceType}
space_type(uₕ::VectorElement) -> Type{<:AbstractSpaceType}

Returns the concrete AbstractSpaceType associated with a VectorElement type or instance.

Examples

Wₕ = gridspace(Ωₕ)
uₕ = element(Wₕ)
space_type(typeof(uₕ)) === typeof(Wₕ)  # true
space_type(uₕ) === typeof(Wₕ)          # true

See also: space, VectorElement

CellAverage(f, nq::Val)

Project f by averaging it over each cell with nq quadrature points per direction. The rule behind avgₕ!.

PointValue(f)

Project f by evaluating it at each mesh point. The rule behind Rₕ!.

ProjectionRule

Supertype of the rules project! can apply: what a grid function's value at a point should be, given a continuous function. PointValue evaluates at the mesh point, CellAverage averages over the cell.

_device_backend(backend::Backend) -> Backend

The backend whose device a projection kernel launches on: the wrapped device backend for a GpuOffload policy, otherwise backend itself.

_device_masked_project!(loc, rule, raw, sp, markers) -> Bool
_device_masked_project!(loc, rule, raws::Tuple, sp, markers, ::Val{NC}) -> Bool

The masked counterparts of _device_project!/_device_scatter_project!: on DeviceLocality, zero the destination(s) on the device and fill only the points in the union of the markers regions, through the rule's own device kernel launched over _marker_index_list. Always false on HostLocality, where the generic sweep folds the mask into its kernel instead.

_device_project!(loc, rule, raw, sp, sel = nothing) -> Bool

Attempt a dedicated device kernel filling raw (a leaf's coefficient vector) according to rule on space sp. Returns true when it did; false when no such kernel exists for this rule/sp combination on this locality, so the caller runs the generic _sweep_for! sweep instead. Always false on HostLocality: the CPU sweep already indexes raw directly and never needs this path. Legality here is a locality question, not a policy question (gpena/Bramble.jl#298): a device kernel either exists for this destination's storage or it doesn't, regardless of which strategy the policy names within that locality. A device index vector sel restricts the fill to those linear indices and leaves every other entry of raw untouched.

_device_scatter_project!(loc, rule, raws::Tuple, sp, ::Val{NC}, sel = nothing) -> Bool

The scatter counterpart of _device_project!: a dedicated device kernel filling every leaf array in raws (NC of them, all sharing sp's mesh) in one launch. Same true/false contract, keyed on locality for the same reason.

_marker_index_list(like, Ωₕ, markers) -> AbstractVector{Int32}

The linear indices of Ωₕ's points lying in any of the markers regions, gathered on the host from index_in_marker's masks and copied into a new Int32 vector allocated like like (so on like's device).

Throws

• ErrorException: the mesh has more points than typemax(Int32).

_offload_backend(backend::Backend) -> Union{Backend, Nothing}

The device backend a GpuOffload policy wraps, or nothing for any other policy, so project! knows whether to fill through the device.

_offload_project!(db, rule, raw, sp, markers) -> Bool
_offload_project!(db, rule, raws::Tuple, sp, markers, ::Val{NC}) -> Bool

Fill the host destination(s) through device backend db: allocate a device buffer per destination, run the rule's device kernel on it (masked when markers is non-empty) and copy it back. false, with the destination untouched, when the rule has no device kernel or the destination's element type is not the device's (_offload_representable).

_on_device(like, x) -> AbstractVector or Tuple

x (a mesh coordinate vector, or a tuple of them) as device arrays allocated like like: returned as is when already device-resident, otherwise copied into a new device array.

_rule_component(rule, k) -> ProjectionRule

rule restricted to leaf k of a composite space whose leaves do not share one mesh, so there is no shared grid point to evaluate once and scatter.

_rule_kernel(rule, space) -> callable

A concretely typed callable mapping a linear grid index to the scalar value rule gives at that point. Kept a named struct per rule rather than a closure: measured, an anonymous closure over the captures here takes a miscompiled path that allocates per grid point (gpena/Bramble.jl#64).

_rule_scatter_kernel(rule, space, ::Val{NC}) -> callable

As _rule_kernel, for a rule whose function returns all NC leaf values at once, so it is evaluated once per point and scattered across the leaves.

project!(uₕ::VectorElement, rule, markers = ()) -> uₕ

Project onto uₕ in place according to rule, optionally restricted to the union of the labelled marker regions, leaving every other entry zero.

Handles the space's shape and the execution policy; rule supplies only the per-point value. See PointValue and CellAverage, and Rₕ! / avgₕ! for the public spellings.

_device_project!(::DeviceLocality, rule::PointValue, raw::AbstractVector, sp::ScalarGridSpace{1}, sel = nothing) -> Bool
_device_project!(::DeviceLocality, rule::PointValue, raw::AbstractVector, sp::ScalarGridSpace{D}, sel = nothing) where {D} -> Bool

Fills raw with rule.f evaluated at every point of the mesh mesh(sp), via a device kernel (gpena/Bramble.jl#94, #174). The 1D method reads the mesh's own coordinate vector directly; the D-dimensional method (D here is never 1, since the method above is strictly more specific and wins dispatch for it) builds each point from the D per-axis coordinate vectors instead. rule.f runs on the device either way: see _gpu_for! for what that requires of it. A device index vector sel (a masked call, see project!) restricts the fill to those linear indices, leaving the rest untouched.

_device_scatter_project!(::DeviceLocality, rule::PointValue, raws::Tuple, sp, ::Val{NC}, sel = nothing) -> Bool

The scatter counterpart of _device_project! above, for an NC-component composite space sp whose leaves share one mesh: rule.f is evaluated once per point and its NC components scattered into raws in the same device kernel. sp is typed generically (not ScalarGridSpace) because the caller is always the composite branch of project!; the dimension the mesh has determines which launcher runs (1D or the D-dimensional counterpart), checked on mesh(sp) at runtime since it cannot be expressed as a type constraint on the composite space itself. sel restricts the fill to a marker index list as in _device_project!.

_launch_restriction!(v::AbstractVector, pts::AbstractVector, f, dev, sel) -> Nothing

Fills v with f evaluated at every point of a 1D mesh's coordinate vector pts (v[i] = f(pts[i])), or only at the linear indices in sel when sel is a device index vector rather than nothing (a masked call, see project!), via a KernelAbstractions.@kernel launch on dev, filled by ext/BrambleKernelAbstractionsExt.jl. f runs on the device, so it must be GPU-compilable in the sense Rₕ!'s docstring describes.

Throws

• ErrorException: no KernelAbstractions extension is loaded, so there is no device kernel to reach (_throw_no_ka_projection_kernel).

_launch_restriction_nd!(v::AbstractVector, pts::Tuple, idxs, f, dev, sel) -> Nothing

The D >= 2 counterpart of _launch_restriction!: builds each grid point from the per-axis coordinate vectors pts at Cartesian index idxs[i] and fills v[i] = f(point), via a KernelAbstractions.@kernel launch on dev. sel restricts it to a marker index list as in _launch_restriction!. Same GPU-compilability requirement on f.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_projection_kernel).

_launch_restriction_scatter!(mats::Tuple, pts::AbstractVector, f, dev, sel) -> Nothing

The scatter counterpart of _launch_restriction!: evaluates f(pts[i]) once per point of a 1D mesh and scatters its components into the destination tuple mats, via a KernelAbstractions.@kernel launch on dev. sel restricts it to a marker index list as in _launch_restriction!. Same GPU-compilability requirement on f.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_projection_kernel).

_launch_restriction_scatter_nd!(mats::Tuple, pts::Tuple, idxs, f, dev, sel) -> Nothing

The D >= 2, scatter counterpart of _launch_restriction!: builds each grid point from the per-axis coordinate vectors pts at Cartesian index idxs[i], evaluates f at it once, and scatters its components into the destination tuple mats, via a KernelAbstractions.@kernel launch on dev. sel restricts it to a marker index list as in _launch_restriction!. Same GPU-compilability requirement on f.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_projection_kernel).

AVG_QUAD_POINTS

Default number of Gauss-Legendre points per direction, per cell, used by avgₕ. Six points are exact for polynomials up to degree eleven.

Unlike an adaptive rule, a fixed one does not tighten itself on coarse cells, so the default is chosen to be accurate on cells far coarser than any practical grid. Measured on 4 points spanning [-1, 4] with a function varying by a factor of e^5 across the domain (deliberately harsher than a real mesh), the worst error over 30 random grids was

points   1D        2D        3D        evaluations per cell (3D)
3        6.1e-5    1.6e-4    3.1e-4     27
4        2.8e-7    8.1e-7    3.4e-6     64
5        1.9e-9    4.1e-9    7.8e-9    125
6        5.5e-12   1.4e-11   1.6e-11   216

Cost is quad_points^D evaluations per cell. On a fine grid three points are usually ample; lower it with the quad_points keyword when the integrand is cheap to resolve and the cells are small.

_cell_average(f, x::AbstractVector, idx::CartesianIndex{1}, nodes::NTuple{NQ, T}, wts::NTuple{NQ, T}) where {NQ, T}
_cell_average(f, x::AbstractVector, i::Int, nodes::NTuple{NQ, T}, wts::NTuple{NQ, T}) where {NQ, T}
_cell_average(f, x::NTuple{D}, idx::CartesianIndex{D}, nodes::NTuple{NQ, T}, wts::NTuple{NQ, T}) where {D, NQ, T}

NQ-point Gauss-Legendre average of f over the cell spanned by the half points around idx, along every axis at once. Scalar or composite (NTuple) by f's own return type: broadcasting .+/.* over a Number accumulator computes exactly the scalar arithmetic +/* would, and over an NTuple accumulator the per-component sum, so no ::Val{NC} or caller-built seed is needed to tell the two apart (gpena/Bramble.jl#102, #148). The accumulator is seeded with the first weighted term itself rather than zero(T), and the loop runs over the remaining nodes, so f is called exactly NQ^D times per cell (gpena/Bramble.jl#693) and the result is the promotion of T (the mesh's element type) and f's return type: a ForwardDiff.Dual under AD stays a Dual (gpena/Bramble.jl#148), and an f narrower than the mesh type (Int, or Float32 on a Float64 mesh) gives a sum whose type does not change mid-loop (gpena/Bramble.jl#523).

A one-dimensional mesh answers half_points with a plain vector but indexes with CartesianIndex{1}; the CartesianIndex{1} method unwraps that to the Int the x::AbstractVector method takes, which is the only thing the 1D/D-dimensional split in this file ever did (gpena/Bramble.jl#69). The 1D, 2D and 3D methods are specialised for the meshes this package builds (D <= 3); the generic NTuple{D} method exists for dispatch correctness at any D and is tested directly rather than through a mesh.

Arguments

• f: Function evaluated at quadrature points, returning a scalar or an NTuple.

• x: Half points along each axis (AbstractVector in 1D, NTuple{D} otherwise).

• idx/i: The cell's CartesianIndex (or, in 1D, its Int).

• nodes, wts: Gauss-Legendre nodes and weights on [0, 1] from _gauss_rule.

Returns

• The cell average of f, matching f's own return shape.

_device_project!(::DeviceLocality, rule::CellAverage, raw::AbstractVector, sp::ScalarGridSpace{1}, sel = nothing) -> Bool
_device_project!(::DeviceLocality, rule::CellAverage, raw::AbstractVector, sp::ScalarGridSpace{D}, sel = nothing) where {D} -> Bool

Fills raw with the cell average of rule.f over every cell of the mesh mesh(sp), via a device kernel that calls the same _cell_average quadrature the CPU sweep uses (gpena/Bramble.jl#94, #174) – the 1D method reads half_points(Ωₕ) directly, the D-dimensional one (never reached for D == 1, since the method above is strictly more specific) hands it the per-axis tuple _cell_average's own 2D/3D methods already expect. rule.f runs on the device either way: see _gpu_for! for what that requires of it. A device index vector sel (a masked call, see project!) restricts the fill to those linear indices, leaving the rest untouched.

_device_scatter_project!(::DeviceLocality, rule::CellAverage, raws::Tuple, sp, ::Val{NC}, sel = nothing) -> Bool

The scatter counterpart of _device_project! above, for an NC-component composite space sp whose leaves share one mesh. sp is typed generically for the same reason restriction.jl's counterpart is: the mesh's own dimension, checked on mesh(sp) at runtime, picks the 1D or D-dimensional launcher. sel restricts the fill to a marker index list as in _device_project!.

_gauss_rule(::Val{N}, ::Type{T})

Returns (nodes, weights) for the N-point Gauss-Legendre rule on [0, 1] as NTuple{N,T}, so the per-cell loop that consumes them does not allocate.

The rule is built by QuadGK.gauss in the requested element type, so Float32 and BigFloat grids get a rule at their own precision rather than a rounded Float64 one. Weights sum to one, which makes the weighted sum over a cell the cell average directly.

_launch_cell_average!(v::AbstractVector, x::AbstractVector, nodes, wts, f, dev, sel) -> Nothing

Fills v[i] with the _cell_average of f over the 1D cell around grid point i, using the Gauss-Legendre nodes/wts and the mesh's half points x – at every i, or only at the linear indices in sel when it is a device index vector rather than nothing (a masked call, see project!) – via a KernelAbstractions.@kernel launch on dev, filled by ext/BrambleKernelAbstractionsExt.jl. f runs on the device: same GPU-compilability requirement as _gpu_for!.

Throws

• ErrorException: no KernelAbstractions extension is loaded, so there is no device kernel to reach (_throw_no_ka_projection_kernel).

_launch_cell_average_nd!(v::AbstractVector, x::Tuple, idxs, nodes, wts, f, dev, sel) -> Nothing

The D >= 2 counterpart of _launch_cell_average!: fills v[i] with the _cell_average of f over the cell at Cartesian index idxs[i], using the per-axis half-point vectors x and the Gauss-Legendre nodes/wts, via a KernelAbstractions.@kernel launch on dev. Same GPU-compilability requirement on f.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_projection_kernel).

_launch_cell_average_scatter!(mats::Tuple, x::AbstractVector, nodes, wts, f, dev, sel) -> Nothing

The scatter counterpart of _launch_cell_average!: computes the _cell_average of f over each 1D cell once and scatters its components into the destination tuple mats, via a KernelAbstractions.@kernel launch on dev. Same GPU-compilability requirement on f.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_projection_kernel).

_launch_cell_average_scatter_nd!(mats::Tuple, x::Tuple, idxs, nodes, wts, f, dev, sel) -> Nothing

The D >= 2, scatter counterpart of _launch_cell_average!: computes the _cell_average of f over the cell at Cartesian index idxs[i] once and scatters its components into the destination tuple mats, via a KernelAbstractions.@kernel launch on dev. Same GPU-compilability requirement on f.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_projection_kernel).

ShiftNode{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a stencil shift operation by shift_amount grid points in dimension Dim.

Built by shift_op and, with shift_amount \(\pm 1\), by the public shifts S₊ₓ(op), S₋ᵧ(op), ... on a symbolic operand. A neighbour off the grid reads 0, as for a grid function.

⊗(A, B)

Kronecker product operator (alias for kron).

Computes the Kronecker product (tensor product) of matrices A and B. This operator is used extensively in constructing multidimensional shift operators.

Examples

I₂ = I(2)
I₃ = I(3)
result = I₂ ⊗ I₃  # 6×6 identity matrix

See also: shift

_Eye(be::Backend, npts, ::Val{i})

Internal helper to create identity or shifted diagonal matrices, in the matrix type the backend be chose, routed through backend_eye/matrix_type(be) rather than a package-specific lazy type, so the result always matches whatever matrix_type the caller's backend picked.

Arguments

• be::Backend: the backend whose matrix_type the result is built in.

• npts::Int: Size of the square matrix.

• ::Val{i}: Diagonal offset (0 = main diagonal, 1 = superdiagonal, -1 = subdiagonal).

Returns

• For i=0: Identity matrix of size npts × npts.

• For i≠0: Matrix with ones on the i-th diagonal, zeros elsewhere.

See also: shift

_shift_band!(out, in_ref, dims::NTuple{D, Int}, dir::GridDirection, dim_val::Val, nbands::Int, b::Int) -> Nothing

Writes the index shift of in_ref along dim_val in direction dir (Forward reads u_{i+1}, Backward reads u_{i-1}) into out, on the b-th of nbands slabs of the grid cut along its last axis, as _average_band! does for the average. The slice with no neighbour gets 0. nbands == 1 is the whole grid, which is the serial engine.

kronecker_operator_matrix(Ωₕ, op)

The Kronecker-product construction stencil_matrix replaces in every operator family (gpena/Bramble.jl#185): op is one of the public per-axis aliases (D₋ₓ, jumpᵧ, M₂, ...) and this returns the same matrix built the old way, out of shift and its combinations (difference_shift, add_half_shift).

Kept, rather than removed with the rest of the Kronecker construction (gpena/Bramble.jl#185), as the oracle stencil_matrix is checked against: the equality test needs an independent construction, not a deprecated one. Declared here; a dispatch method is added next to each family's own implementation (difference.jl, average.jl, jump.jl), mapping its public aliases to the _kron_* construction that family kept.

shift(Ωₕ::AbstractMeshType, ::Val{SHIFT_DIM}, ::Val{i})

Returns the matrix that shifts a grid function by i points along direction SHIFT_DIM, in the backend's matrix_type, over the flattened degrees of freedom of Ωₕ.

Arguments

• SHIFT_DIM: the direction to shift along, 1 for \(x\), 2 for \(y\), 3 for \(z\).

• i: how far to shift. 1 is the superdiagonal, -1 the subdiagonal, 0 the identity. The stencil is truncated at the boundary rather than wrapped, so the matrix has n - |i| nonzeros per direction rather than n.

Tensor-product structure

A mesh is a tensor product of its one-dimensional meshes, and its degrees of freedom are flattened in column-major order, so a shift along one direction is the identity in every other direction. Writing \(E_k\) for the identity of size \(n_k\) and \(S_k(i)\) for the one-dimensional shift by i, the operator is a Kronecker product with \(S\) in one slot:

_recursive_shift builds exactly this, recursing from the outermost dimension inwards and placing \(S\) when it reaches SHIFT_DIM. The per-direction forms it generalises, each of which the test suite checks against shift:

# 1D
shift(Ωₕ, Val(1), Val(i))  ==  _Eye(be, nₓ, Val(i))

# 2D, on an nₓ × n_y grid
shift(Ωₕ, Val(1), Val(i))  ==  backend_eye(be, n_y) ⊗ _Eye(be, nₓ, Val(i))
shift(Ωₕ, Val(2), Val(i))  ==  _Eye(be, n_y, Val(i)) ⊗ backend_eye(be, nₓ)

# 3D, on an nₓ × n_y × n_z grid
shift(Ωₕ, Val(3), Val(i))  ==  _Eye(be, n_z, Val(i)) ⊗ backend_eye(be, nₓ * n_y)

i == 0 short-circuits to the identity of the whole grid without building any Kronecker product.

It builds the S₊/S₋ matrix forms and the Kronecker oracle kronecker_operator_matrix; the operator families' own matrices come from stencil_matrix.

See also: ⊗, diff₋ₓ.

shift_op(op::LazyOp{D}, ::Val{Dim}, amount::Int) where {D, Dim}
shift_op(op::LazyOp{D}, dim::Int, amount::Int) where D

Shifts the stencil of op by amount grid points in dimension Dim (or dim). A direction that is a Bool, not an integer, or outside 1:D throws an ArgumentError. The Int direction forwards to the Val form, so a literal direction still gives a concrete node type.

OperatorArgument

The three things a discrete operator accepts: a mesh, a grid space or a VectorElement. Exactly the three _op_mesh has a method for.

Named so that the dimensional entry points can say it. The per-coordinate aliases take an untyped arg, which is harmless at one argument, but D₋(arg, d) shares its name with the form layer's D₋(op::LazyOp, ::Val): an untyped two-argument method would swallow a LazyOp and hand it to the grid engines, and an Int direction on a symbolic node cannot be type-stable anyway (BackwardDifference{D, 1, …} and {D, 2, …} are different types). Spelling the union keeps that a MethodError (gpena/Bramble.jl#74).

StencilOp{Dim}

One operator family along axis Dim, as stencil_matrix takes it. Each concrete subtype answers _stencil_taps (its fixed neighbour offsets along Dim) and _stencil_weights (the coefficients matching those offsets at a grid point).

_alias_bang_expr(base_op_name, alias_name, dir_string, suffix, direction_index, what,
                 formula; opening_sentence = "", source = nothing)

Returns the expression defining alias_name(vₕ, uₕ) as base_op_name(vₕ, uₕ, Val(direction_index)), with its docstring attached.

The in-place sibling of _alias_expr. Two generators rather than one because the shapes differ: the allocating alias takes a single argument that may be a mesh, a space or a grid function, while this one takes a destination and a source and is only ever about grid functions.

opening_sentence, given non-empty, replaces the generic "The $dir_string $what of uₕ along the $suffix direction, \($formula\), written into vₕ." the same way it does for _alias_expr – Dc!/D̽! have no backward/forward adjective to put in dir_string either.

source is the LineNumberNode the generated method is attributed to; the macro passes its own __source__, so a generated method reports the family's call site rather than the quoted line in this file.

_alias_expr(base_op_name, alias_name, dir_string, suffix, direction_index, what,
            formula; opening_sentence = "", formula_note = "", alias_note = "",
            trailing_note = "", source = nothing)

Returns the expression defining alias_name(arg) as base_op_name(arg, Val(direction_index)), with its docstring attached.

what names the quantity, such as "finite difference", and formula is the LaTeX for it. Both are needed because the four operator families share this generator: describing every alias as a "difference" would be wrong for the averages, and would not separate the unscaled difference from the finite difference.

opening_sentence, given non-empty, replaces the generic "The $dir_string $what along the $suffix direction, \($formula\)." with the caller's own wording: Dc and D̽ have no backward/forward adjective to put in dir_string at all.

The three remaining keyword notes are each a sentence the docstring includes only when given (non-empty), one per insertion point a family may need: formula_note follows the opening sentence (the diff/finite-difference families use this to contrast the two, which does not apply to an average); alias_note follows the Alias for ... sentence, before arg is described (D̽ uses this to compare itself with Dc); trailing_note follows the description of arg, before the closing "Accepts a grid function..." paragraph (D̃, Dc and D̽ use this for their boundary-behaviour and precondition caveats, which differ both in what happens at the ends – D̃/Dc truncate, D̽ falls back to a one-sided difference (gpena/Bramble.jl#183) – and in whether a mesh needs at least three points along the direction).

source attributes the generated method to the family's own call site, as in _alias_bang_expr.

_dimensional_expr(base_op_name, alias_name, alias_stem, dir_string, what, formula;
                  opening_sentence = "", source = nothing)

Returns the expressions defining one family's dimensional entry point, as a Vector{Expr}: alias_name(arg, d) for d a Val, an Int or a Symbol.

Every family bottoms out in a Val-parameterised base function; this adds the entry point with the direction as a value, so sum(D₋(uₕ, d) for d in 1:D) is dimension-agnostic (gpena/Bramble.jl#74).

The Val method is omitted when alias_name is base_op_name itself, as it is for jump: the family's own file already wrote that method, and generating it again would define jump(arg, ::Val) = jump(arg, ::Val).

alias_stem is the family's subscript stem, which the docstring cross-references: it is alias_name for most families, but the averages put their entry point on the tuple-valued Mₕ/M₊ₕ while the per-coordinate aliases are still Mₓ/M₊ₓ, and Mₕₓ is not a name.

opening_sentence is the family's own prose template, substituted with d as the direction; the families that do not override it get the generic "The $dir_string $what along the d direction, \($formula\)." instead.

_grid_function_forms_expr(base_name, apply_fn, extra_args, dir_instance;
                          docstring = "", source = nothing)

Returns the expressions defining the three methods every directional operator family applies a grid function through, as a Vector{Expr}: base_name! on a scalar grid function, base_name! on a composite one, and the allocating base_name built on top of them.

The three are byte-identical across the families apart from which applicator they call and what it takes before the direction (gpena/Bramble.jl#101):

extra_args are spliced as bare identifiers, so each generated method names an ordinary top-level function rather than closing over one: that is what keeps every call site specialising to its own zero-allocation method, as the hand-written versions did.

The composite method recurses into the scalar one through apply_fn's own composite method rather than repeating the walk, so a leaf's own submesh is what each leaf is measured against (gpena/Bramble.jl#79).

docstring, given non-empty, is attached to the scalar base_name! method; the families whose prose lives here rather than on a separately hand-written matrix form use it.

_operator_aliases_expr(base_name, alias_stem, dir_string, what, formula;
                       dispatch_alias = alias_stem, vectorial_alias = nothing,
                       vectorial_dir_string = dir_string, vectorial_what = what,
                       vectorial_note = "", opening_sentence = "", formula_note = "",
                       alias_note = "", trailing_note = "", bang_opening_sentence = "",
                       source = nothing)

Returns the expressions defining one family's whole alias surface, as a Vector{Expr}: the per-coordinate alias_stem pair for every direction (Dcₓ/Dcₓ!, Mᵧ/Mᵧ!, …), the dimensional entry point dispatch_alias(arg, d), and, when vectorial_alias is given, the ₕ alias over every coordinate at once.

dispatch_alias defaults to alias_stem, which is what every difference family wants: D₋ₓ, D₋ᵧ, D₋₂ and D₋(uₕ, d). The averages pass their tuple-valued alias instead, so that the family gains a direction argument without minting a bare M – the single most common local name in finite-element code, and one using Bramble would then reserve. Passing nothing skips the entry point altogether.

The five note keywords are prose templates rather than finished sentences: each may carry {direction} and {suffix}, substituted per alias by _subst. They were closures taking (direction, suffix) until gpena/Bramble.jl#258 – a macro sees the expression that would build a closure, never the closure itself, so the notes a family needs ("over the averaged spacing", "second order on a non-uniform grid, where Dc{suffix} is first") travel as string literals instead.

vectorial_dir_string/vectorial_what default to dir_string/what and exist for the families that describe the tuple-valued alias differently from the per-coordinate ones: D̃/Dc/D̽ override every directional opening_sentence and so pass dir_string and what empty, while D̃ₕ/Dcₕ/D̽ₕ still want "the centered difference of arg along every coordinate".

_relocate!(ex, source)

Rewrites every LineNumberNode in ex to source, and returns ex.

A generated method otherwise reports the :(...) quote it was built from – one shared line in this file for every family – because that is the line info the quote carries. Rewriting it to the __source__ the macro was expanded at makes methods(D₋ₓ), a stacktrace and an editor's "go to definition" all land on the family's own @operator_family call.

source given as nothing leaves ex untouched, which is what a direct call outside the macro gets.

_subst(template, direction, suffix)

Substitutes the {direction} and {suffix} placeholders in a family's prose template.

Both placeholders are replaced in a single pass, so neither can be rewritten by the other's replacement text. An empty template stays empty, which is how a family says it needs no note at that insertion point.

_vectorial_expr(base_op_name, alias_name, dir_string, what; note = "", source = nothing)

Returns the expression defining the ₕ alias that applies base_op_name along every coordinate and returns a tuple, one entry per spatial dimension, as a Vector{Expr}. On a one-dimensional mesh the alias returns that single entry rather than a one-tuple.

The counterpart of _alias_expr for the tuple-valued aliases (∇ₕ, diff₋ₕ, Mₕ). The operator families generated the same three methods independently before this existed.

note, given non-empty, is an extra sentence appended after the worked 2D example – D̽ₕ uses this to place itself relative to ∇ₕ/∇₊ₕ, a comparison none of the other vectorial aliases need.

_vectorial_index_expr(alias_name, alias_stem, source)

Returns the expressions giving one family's vectorial alias (∇ₕ, Mₕ, jumpₕ, ...) Julia's iteration and indexing protocol, as a Vector{Expr} (gpena/Bramble.jl#340): Base.iterate (one- and two-argument), Base.length, Base.eltype, Base.firstindex, Base.lastindex and Base.getindex on an Integer or a Symbol (:x/:y/:z).

Every method reads off the same three per-coordinate aliases (D₋ₓ, D₋ᵧ, D₋₂, from alias_stem) that alias_name(arg) itself calls through _vectorial_apply, so dx, dy = ∇ₕ and D₋ₓ === ∇ₕ[1] === ∇ₕ[:x] are the same function object, not a copy. The three names are spliced in as literals, so ∇ₕ[2] is a call to a @inline function of a literal Int against literal comparisons – the compiler constant-folds it to the concrete coordinate operator, same as any other three-armed literal dispatch in this file.

Always three coordinates regardless of the mesh a caller eventually applies the result to: the alias itself is dimension-agnostic (typeof(∇ₕ), not tied to any mesh), so indexing it out of range along a direction the mesh does not have is caught downstream, by the operator the index yields, the same way ∇ₕ(u) itself is.

@operator_family(kwargs...)

Defines one operator family's grid-function forms and its whole alias surface, declaratively.

The methods arrive by macro expansion, so the parser sees each definition where it is written; nothing is Core.evaled at load time (gpena/Bramble.jl#258).

Keywords, all optional except base and stem:

Every value is a literal: a symbol for a name, a string for prose, a call such as Backward() for the direction tag. Nothing is evaluated at expansion time beyond assembling the docstrings, which is ordinary string building over those literals.

_stencil_taps(op::StencilOp) -> NTuple{K,Int}

The fixed neighbour offsets along op's axis that its stencil reads, in the same order _stencil_weights returns their coefficients. A tap whose neighbour falls outside the grid contributes no entry at all (checked in stencil_matrix). A tap that stays on the grid is always stored, even where its weight is exactly zero: the stencil, not the value, decides the pattern.

_stencil_weights(op::StencilOp{Dim}, Ωₕ::AbstractMeshType, I::CartesianIndex) -> NTuple{K}

The coefficients matching _stencil_taps(op) at I. The unscaled families' self-weight is always ±1 regardless of boundary – their truncation is entirely the "tap falls outside the grid" exclusion stencil_matrix applies – while the scaled and averaged families additionally zero their own weight at a boundary that a Kronecker construction would reach through a since-zeroed weight vector rather than exclude, which this reproduces explicitly.

Ωₕ is typed as Union{AbstractMeshType, _HostAxisSpacings} rather than AbstractMeshType alone: every method below reads its mesh through _axis_spacing/_axis_forward_spacing/ _axis_npoints/_axis_eltype only, and _HostAxisSpacings answers exactly those four without claiming to be a mesh (see the comment above it).

stencil_matrix(Ωₕ::AbstractMeshType, op::StencilOp)
stencil_matrix(Wₕ::AbstractSpaceType, op::StencilOp)

Builds op's operator matrix in one pass over Ωₕ's grid points, in the matrix type matrix_type(backend(Ωₕ)) picked: a SparseMatrixCSC backend writes colptr/rowval/ nzval directly (two passes, count then fill; no intermediate matrix, no broadcast), any other AbstractMatrix backend falls back to a plain dense fill, built on the host and handed to the backend's matrix type in one copyto! – never one scalar write per stored entry into device memory (gpena/Bramble.jl#94). Ωₕ may be a mesh or a grid space, taken as mesh(Wₕ).

Checked entrywise, nnz included, against kronecker_operator_matrix, the Kronecker-product construction every operator family used before (gpena/Bramble.jl#185).

BackwardDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a backward finite difference operator acting in dimension Dim.

CenteredDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the centered difference along Dim,

Truncated at both ends of Dim, having no neighbour on one side.

CrossWeightedDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the cross-weighted centered difference along Dim,

The same two one-sided differences the centered difference combines, weighted by the opposite spacings. That swap is what makes it second order on a non-uniform grid where Dc is first, and the two coincide when the spacing is constant. Neither end is truncated: each collapses to the one-sided difference its near side still defines, D₊ at the first point and D₋ at the last.

ForwardDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a forward finite difference operator acting in dimension Dim.

SecondDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the second difference along Dim, the composition \(\tilde{D}_{+}(D_{-}(u))\) as one three-point stencil,

Zero at the last point of Dim, where D̃ is truncated, and without the backward tap at the first, where D₋ is: entry for entry the composition D̃ₓ(D₋ₓ(u)), built by Δₕ once per axis. One node, so its operand is evaluated once per tap rather than once per tap of a nested inner stencil.

StarDifference{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node for the forward difference over the averaged spacing along Dim,

The forward difference over the averaged spacing rather than the forward one, which is what makes the discrete integration by parts close. Truncated at the far end of Dim.

StarSpacings(h)

Lazy view of the averaged spacings \((h_i + h_{i+1})/2\) over a mesh's cached backward spacings h, which is what D̃ₓ divides by.

Entry i reads h[i] and h[i+1], so it is defined for i < length(h). That is exactly the range the forward stencil's interior covers; the last point has no forward neighbour and the engine truncates it to zero without consulting this.

Bramble._CenteredStrainTensor{D}

Builder-only container for εcₕ's and ε̽ₕ's D × D entries, each held as its additive pieces. Consumed only by innerₕ on two of these, never part of an assembled AST.

Bramble._DivergenceTerms{D}

Builder-only container for the D per-component terms of divₕ and divₘₕ: the symbolic twin of Bramble._StrainTensor for the divergence, and with the same lifetime – consumed immediately by inner₊, never part of an assembled AST.

Bramble._GradientRows{D, N, T}

The N-tuple of D-tuples ∇ₕ and ∇₊ₕ return for a composite trial or test function, held so that inner₊ can name the component of each row. Indexing g[c][d], length, iteration, map and broadcasting give what the tuple it holds gives, and Tuple(g) is that tuple; anything else tuple-like needs Tuple(g). map takes only these containers: with a plain tuple among its arguments, pass Tuple(g). It is not a LazyOp and never reaches assembly: inner₊ on two of them builds the sum of the rows' products, held in a Bramble._ComponentSum that names each row's component pair.

Bramble._StrainTensor{D}

Builder-only container for εₕ's D × D entries, each held as its additive pieces (Bramble._strain_pieces) rather than their sum. Never reaches assembly or an AST walker: inner₊ on two of these expands immediately into the sum of single-block products a user would write out by hand, held in a Bramble._ComponentSum that names each term's component pairs.

_batch_difference_engine!(out, in_ref, h, dims::Tuple, dir::GridDirection, dim_val::Val) -> Nothing

CpuPolyester's _difference_engine!, filled by BramblePolyesterExt (one _difference_band! per band). Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_difference_band!(out, in_ref, h, dims::NTuple{D, Int}, dir::GridDirection, dim_val::Val, nbands::Int, b::Int) -> Nothing

Runs _difference_engine!'s loops for dir on the b-th of nbands slabs of the grid cut along its last axis (_band_range). The slabs partition the grid, so running every band, in any order or concurrently, writes every point of out exactly once with the value the serial engine gives it.

_difference_band_adjoint!(d_in, d_out, h, dims::NTuple{D, Int}, dir::GridDirection, dim_val::Val, nbands::Int, b::Int) -> Nothing

The transpose of _difference_band!, in gather form. For every point j of the b-th of nbands slabs, adds to d_in[j] what d_out carries back through each output that reads in_ref[j]: the one at j and its neighbours along DIM. Only the band owning j writes d_in[j], so every band may run concurrently. d_out is read, never cleared.

_difference_band_spacing_tangent!(d_out, in_ref, h, d_h, dims::NTuple{D, Int}, dir::GridDirection, dim_val::Val, nbands::Int, b::Int) -> Nothing

Adds to d_out, on the b-th of nbands slabs of the grid, the tangent of _difference_band!'s output along the spacing tangent d_h, read as h is (_get_h_val). Each point is written by the band owning it.

_difference_spacing_adjoint!(d_h, in_ref, d_out, h, dims::NTuple{D, Int}, dir::GridDirection, dim_val::Val, nbands::Int, b::Int) -> Nothing

The cotangent of the spacing values that _difference_band! reads, in gather form. For every q of the b-th of nbands bands of eachindex(d_h), adds to d_h[q] the sum, over the outputs reading h(q) (the slice at q along DIM, and for D̽ₕ the interior slice at q - 1), of d_out times the output's partial derivative in h(q). d_h is indexed as h is. Only the band owning q writes d_h[q].

_launch_difference_centered!(out::AbstractVector, in_ref, h, dims::Tuple, dir::CenteredStencil, dim_val::Val, dev) -> Nothing

Fills out with the centered (Centered/CrossWeighted) difference of in_ref along the axis dim_val, scaled by h. At the two boundary slices it defers to dir's own CenteredStencil boundary rule in _compute_difference – zero for Centered, the one-sided difference the near side still defines for CrossWeighted. Runs via a KernelAbstractions.@kernel launch on dev. The device counterpart of the centered branch of _difference_engine!'s CPU sweep above.

Throws

• ErrorException: no KernelAbstractions extension is loaded (_throw_no_ka_stencil_kernel).

_launch_difference_onesided!(out::AbstractVector, in_ref, h, dims::Tuple, dir::GridDirection, dim_val::Val, dev) -> Nothing

Fills out with the one-sided (Forward/Backward) difference of in_ref along the axis dim_val, scaled by h when given (nothing for the unscaled difference), truncating the one boundary slice with no neighbour to zero, matching _compute_difference's GridDirection methods. Runs via a KernelAbstractions.@kernel launch on dev, filled by ext/BrambleKernelAbstractionsExt.jl. The device counterpart of the one-sided branch of _difference_engine!'s CPU sweep above.

Throws

• ErrorException: no KernelAbstractions extension is loaded, so there is no device kernel to reach (_throw_no_ka_stencil_kernel).

_star_spacings_adjoint!(d_h, d_s, nbands::Int, b::Int) -> Nothing

Adds the cotangent d_s of a StarSpacings view's entries to the cotangent d_h of the spacings it averages, for every q of the b-th of nbands bands of eachindex(d_h). Entry i of the view is \((h_i + h_{i+1})/2\), so d_h[q] gathers half of d_s[q - 1] and of d_s[q], each where it exists.

_threaded_difference_engine!(out, in_ref, h, dims::Tuple, dir::GridDirection, dim_val::Val) -> Nothing

CpuThreaded's _difference_engine!: one _difference_band! per thread under Threads.@threads :static, or every band in turn where a :static loop cannot start (_static_or_serial). Kept in an isolated function, as _threaded_for! is, so the Threads.@threads closure is never built on a serial path.

backward_derivative_weights!(v::AbstractVector, Ωₕ::AbstractMeshType, diff_dim::Val)

Computes the geometric weights for the backward finite difference operator and stores them in-place in vector v.

backward_difference(arg, dim_val::Val)

Constructs the unscaled backward difference operator, representing the operation \(u_{i} - u_{i-1}\).

backward_difference_dim!(out, in, [h], dims, diff_dim)

In-place backward difference of in along diff_dim, written into out: \(u_{i} - u_{i-1}\), divided by h when it is given (the finite difference) and unscaled without it.

backward_finite_difference(arg, dim_val::Val)

Constructs the backward finite difference operator, which approximates the first derivative using the formula \(\frac{u_{i} - u_{i-1}}{h_i}\).

centered_difference(uₕ::VectorElement, dim_val::Val)

The centered difference of uₕ along dim_val:

Reached through Dcₓ and its siblings, and takes a mesh, a grid space or a grid function as the other difference families do.

The denominator is \(x_{i+1} - x_{i-1}\), so the operator reproduces the derivative of an affine function exactly on any grid, uniform or not. Both the first and the last point lack a neighbour on one side, so both are truncated to zero.

See also: star_spacings, D₋ₓ, D₊ₓ.

centered_difference(Ωₕ::AbstractMeshType, dim_val::Val)

The centered difference along dim_val, as a sparse matrix.

Reaches one point either side, so both end rows are zero and the mesh needs at least three points along the direction.

cross_weighted_difference(uₕ::VectorElement, dim_val::Val)

The cross-weighted centered difference of uₕ along dim_val:

Reached through D̽ₓ and its siblings, and takes a mesh, a grid space or a grid function as the other difference families do.

It is the same two one-sided differences Dcₓ combines, weighted by the opposite spacings. That is the combination which cancels the leading truncation term on a non-uniform grid, so this is second order where Dcₓ is first, and the two coincide when the spacing is constant.

The first and the last point each lack a neighbour on one side, but unlike Dcₓ neither is truncated: each collapses to the one-sided difference its near side still defines, D₊ₓ(uₕ) at the first point and D₋ₓ(uₕ) at the last.

See also: Dcₓ, D₋ₓ, D₊ₓ.

cross_weighted_difference(Ωₕ::AbstractMeshType, dim_val::Val)

The cross-weighted centered difference along dim_val, as a sparse matrix.

A three-point stencil in the interior, but neither end row is empty: with no neighbour on the far side, row 1 agrees with D₊ₓ(Ωₕ, dim_val) and row n with D₋ₓ(Ωₕ, dim_val), each under its own diagonal weight alongside the interior cross-weighting. The mesh still needs at least three points along the direction.

diff₊(arg, d)

The forward unscaled difference along the d direction, \(u_{i+1} - u_i\).

d names the direction three ways: a Val (Val(1)), an Int (1, 2, 3) or a Symbol (:x, :y, :z). The Val form is the one diff₊ₓ and its siblings forward through; the other two select between literal Vals, one arm per direction the mesh has, so the direction still reaches the stencil engine as a compile-time constant and a loop over d allocates nothing beyond its results. An Int outside 1:dim(mesh), or a Symbol that is not :x/:y/:z, throws an ArgumentError.

Alias for forward_difference(arg, Val(d)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₊ᵧ!(vₕ, uₕ)

The forward unscaled difference of uₕ along the y direction, \(u_{i+1} - u_i\), written into vₕ.

The in-place form of diff₊ᵧ: it allocates nothing, where the allocating form allocates its result. Returns vₕ, so it composes: normₕ(diff₊ᵧ!(vₕ, uₕ)).

vₕ and uₕ must be grid functions of the same space, and must not be the same object, as every stencil reads neighbours of the target coordinate; aliasing them would read values that have already been overwritten. Either violation throws an ArgumentError.

Alias for forward_difference!(vₕ, uₕ, Val(2)). Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter.

diff₊ᵧ(arg)

The forward unscaled difference along the y direction, \(u_{i+1} - u_i\). The unscaled difference is not divided by the grid spacing; the finite difference is.

Alias for forward_difference(arg, Val(2)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₊₂!(vₕ, uₕ)

The forward unscaled difference of uₕ along the z direction, \(u_{i+1} - u_i\), written into vₕ.

The in-place form of diff₊₂: it allocates nothing, where the allocating form allocates its result. Returns vₕ, so it composes: normₕ(diff₊₂!(vₕ, uₕ)).

vₕ and uₕ must be grid functions of the same space, and must not be the same object, as every stencil reads neighbours of the target coordinate; aliasing them would read values that have already been overwritten. Either violation throws an ArgumentError.

Alias for forward_difference!(vₕ, uₕ, Val(3)). Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter.

diff₊₂(arg)

The forward unscaled difference along the z direction, \(u_{i+1} - u_i\). The unscaled difference is not divided by the grid spacing; the finite difference is.

Alias for forward_difference(arg, Val(3)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₊ₓ!(vₕ, uₕ)

The forward unscaled difference of uₕ along the x direction, \(u_{i+1} - u_i\), written into vₕ.

The in-place form of diff₊ₓ: it allocates nothing, where the allocating form allocates its result. Returns vₕ, so it composes: normₕ(diff₊ₓ!(vₕ, uₕ)).

vₕ and uₕ must be grid functions of the same space, and must not be the same object, as every stencil reads neighbours of the target coordinate; aliasing them would read values that have already been overwritten. Either violation throws an ArgumentError.

Alias for forward_difference!(vₕ, uₕ, Val(1)). Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter.

diff₊ₓ(arg)

The forward unscaled difference along the x direction, \(u_{i+1} - u_i\). The unscaled difference is not divided by the grid spacing; the finite difference is.

Alias for forward_difference(arg, Val(1)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₊ₕ(arg)

The forward unscaled difference of arg along every coordinate, as a tuple with one entry per spatial dimension. On a one-dimensional mesh it returns that single entry rather than a one-tuple.

For a 2D space, diff₊ₕ(uₕ) is (forward_difference(uₕ, Val(1)), forward_difference(uₕ, Val(2))). arg is a mesh, a grid space or a VectorElement, as for forward_difference.

Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter: each entry of the tuple is then itself a composite grid function.

diff₋(arg, d)

The backward unscaled difference along the d direction, \(u_{i} - u_{i-1}\).

d names the direction three ways: a Val (Val(1)), an Int (1, 2, 3) or a Symbol (:x, :y, :z). The Val form is the one diff₋ₓ and its siblings forward through; the other two select between literal Vals, one arm per direction the mesh has, so the direction still reaches the stencil engine as a compile-time constant and a loop over d allocates nothing beyond its results. An Int outside 1:dim(mesh), or a Symbol that is not :x/:y/:z, throws an ArgumentError.

Alias for backward_difference(arg, Val(d)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₋ᵧ!(vₕ, uₕ)

The backward unscaled difference of uₕ along the y direction, \(u_{i} - u_{i-1}\), written into vₕ.

The in-place form of diff₋ᵧ: it allocates nothing, where the allocating form allocates its result. Returns vₕ, so it composes: normₕ(diff₋ᵧ!(vₕ, uₕ)).

vₕ and uₕ must be grid functions of the same space, and must not be the same object, as every stencil reads neighbours of the target coordinate; aliasing them would read values that have already been overwritten. Either violation throws an ArgumentError.

Alias for backward_difference!(vₕ, uₕ, Val(2)). Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter.

diff₋ᵧ(arg)

The backward unscaled difference along the y direction, \(u_{i} - u_{i-1}\). The unscaled difference is not divided by the grid spacing; the finite difference is.

Alias for backward_difference(arg, Val(2)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₋₂!(vₕ, uₕ)

The backward unscaled difference of uₕ along the z direction, \(u_{i} - u_{i-1}\), written into vₕ.

The in-place form of diff₋₂: it allocates nothing, where the allocating form allocates its result. Returns vₕ, so it composes: normₕ(diff₋₂!(vₕ, uₕ)).

vₕ and uₕ must be grid functions of the same space, and must not be the same object, as every stencil reads neighbours of the target coordinate; aliasing them would read values that have already been overwritten. Either violation throws an ArgumentError.

Alias for backward_difference!(vₕ, uₕ, Val(3)). Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter.

diff₋₂(arg)

The backward unscaled difference along the z direction, \(u_{i} - u_{i-1}\). The unscaled difference is not divided by the grid spacing; the finite difference is.

Alias for backward_difference(arg, Val(3)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₋ₓ!(vₕ, uₕ)

The backward unscaled difference of uₕ along the x direction, \(u_{i} - u_{i-1}\), written into vₕ.

The in-place form of diff₋ₓ: it allocates nothing, where the allocating form allocates its result. Returns vₕ, so it composes: normₕ(diff₋ₓ!(vₕ, uₕ)).

vₕ and uₕ must be grid functions of the same space, and must not be the same object, as every stencil reads neighbours of the target coordinate; aliasing them would read values that have already been overwritten. Either violation throws an ArgumentError.

Alias for backward_difference!(vₕ, uₕ, Val(1)). Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter.

diff₋ₓ(arg)

The backward unscaled difference along the x direction, \(u_{i} - u_{i-1}\). The unscaled difference is not divided by the grid spacing; the finite difference is.

Alias for backward_difference(arg, Val(1)). arg is a mesh, a grid space or a VectorElement: the first two give the operator as a sparse matrix, the third applies it and returns a VectorElement.

Accepts a grid function of a scalar or of a composite grid space. On a composite one the operator is applied to each component in turn, and the result is the composite grid function whose components are those results.

diff₋ₕ(arg)

The backward unscaled difference of arg along every coordinate, as a tuple with one entry per spatial dimension. On a one-dimensional mesh it returns that single entry rather than a one-tuple.

For a 2D space, diff₋ₕ(uₕ) is (backward_difference(uₕ, Val(1)), backward_difference(uₕ, Val(2))). arg is a mesh, a grid space or a VectorElement, as for backward_difference.

Accepts a grid function of a scalar or of a composite grid space, componentwise on the latter: each entry of the tuple is then itself a composite grid function.

forward_derivative_weights!(v::AbstractVector, Ωₕ::AbstractMeshType, diff_dim::Val)

Computes the geometric weights for the forward finite difference operator and stores them in-place in vector v.

forward_difference(arg, dim_val::Val)

Constructs the unscaled forward difference operator, representing the operation \(u_{i+1} - u_i\).

forward_difference_dim!(out, in, [h], dims, diff_dim)

In-place forward difference of in along diff_dim, written into out: \(u_{i+1} - u_i\), divided by h when it is given (the finite difference) and unscaled without it.

forward_finite_difference(arg, dim_val::Val)

Constructs the forward finite difference operator, which approximates the first derivative using the formula \(\frac{u_{i+1} - u_i}{h_i}\).

forward_star_difference(uₕ::VectorElement, dim_val::Val)

The forward difference of uₕ along dim_val, divided by the averaged spacing:

Reached through D̃ₓ and its siblings, and takes a mesh, a grid space or a grid function as the other difference families do.

The last point has no forward neighbour, so it is truncated to zero, as in D₊ₓ.

See also: star_spacings, D₊ₓ.

forward_star_difference(Ωₕ::AbstractMeshType, dim_val::Val)

The forward difference over the averaged spacing along dim_val, as a sparse matrix.

The forward difference scaled by the averaged spacing instead of the forward one. The last point along the direction has no forward neighbour, so its row is zero.

star_spacings(Ωₕ::Mesh1D)

Returns the averaged spacings \((h_i + h_{i+1})/2\) of Ωₕ as a StarSpacings view over its cached backward spacings. Allocates nothing.

Away from the first point this equals half_spacing(Ωₕ, i). At i = 1 it does not: the cached h₁ repeats the first interval, so this gives \(x_2 - x_1\) where the cell width gives half of it, the boundary cell being a half cell.

JumpNode{D, Dim, OpType <: LazyOp{D}} <: LazyOp{D}

AST node representing the jump across interfaces along dimension Dim, \(u_{i+1} - u_i\).

Not truncated at the far end: the absent u_{i+1} is taken as zero, yielding -uᵢ there. This matches the space-layer matrix convention, whose boundary row preserves -1.

jump!(vₕ, uₕ, dim_val::Val)

The jump across the interfaces along dim_val, written into vₕ, which is returned.

Forwards to forward_difference!, as the allocating form forwards to forward_difference: a jump across an interface and an unscaled forward difference are the same quantity.

jump_dim!(out, in, dims, jump_dim::Val)

In-place jump of in along jump_dim, written into out, computing \(u_{i+1} - u_i\).

Forwards to forward_difference_dim!; the jump and the unscaled forward difference are the same arithmetic.

BackwardAverage{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a backward spatial averaging operator acting in dimension Dim.

CenteredAverage{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing the centered spatial averaging operator acting in dimension Dim, \((u_{i-1} + 2 u_{i} + u_{i+1})/4\), zero on both end slices of Dim.

ForwardAverage{D,Dim,OpType<:LazyOp{D}} <: LazyOp{D}

An AST node representing a forward spatial averaging operator acting in dimension Dim.

_average_band!(out, in_ref, dims::NTuple{D, Int}, dir::GridDirection, dim_val::Val, nbands::Int, b::Int) -> Nothing

Fills the b-th of nbands slabs of the grid cut along its last axis with the two-point average of in_ref along dim_val in direction dir, zero on the one boundary slice with no neighbour, as _difference_band! does for the difference; 1, 1 is the whole grid, the serial _average_engine!. The traversal is shared with the difference engine; see _stencil_ranges in operators/stencil.jl.

_average_band_adjoint!(d_in, d_out, dims::NTuple{D, Int}, dir::Union{Forward, Backward}, dim_val::Val, nbands::Int, b::Int) -> Nothing

The transpose of _average_band!, in gather form. For every point j of the b-th of nbands slabs, adds to d_in[j] half of d_out at each output that reads in_ref[j]: j itself and its neighbour along DIM, where that output is not the zeroed boundary slice. Only the band owning j writes d_in[j], so every band may run concurrently. d_out is read, never cleared.

_batch_average_engine!(out, in_ref, dims::Tuple, dir::GridDirection, dim_val::Val) -> Nothing

CpuPolyester's _average_engine!, filled by BramblePolyesterExt (one _average_band! per band). Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_batch_centered_average_engine!(out, in_ref, dims::Tuple, dim_val::Val) -> Nothing

CpuPolyester's _centered_average_engine!, filled by BramblePolyesterExt (one _centered_average_band! per band). Its src/ method only throws, naming Polyester; BramblePolyesterExt defines the one that runs.

_centered_average_band!(out, in_ref, dims::NTuple{D, Int}, dim_val::Val, nbands::Int, b::Int) -> Nothing

Runs _centered_average_engine!'s loops on the b-th of nbands slabs of the grid cut along its last axis, as _difference_band! does for the difference.

_centered_average_band_adjoint!(d_in, d_out, dims::NTuple{D, Int}, dim_val::Val, nbands::Int, b::Int) -> Nothing

The transpose of _centered_average_band!, in gather form: for every point j of the b-th of nbands slabs, adds to d_in[j] the weights 1/4, 1/2, 1/4 times d_out at the interior outputs j - 1, j, j + 1 along DIM (both end slices are zero and read nothing). Only the band owning j writes d_in[j]; d_out is read, never cleared.

_launch_average_engine!(out::AbstractVector, in_ref, dims::Tuple, dir::GridDirection, dim_val::Val, dev) -> Nothing

Fills out with the two-point average (_compute_average) of in_ref along the axis dim_val, in direction dir (Forward/Backward), truncating the one boundary slice with no neighbour to zero, via a KernelAbstractions.@kernel launch on dev, filled by ext/BrambleKernelAbstractionsExt.jl. The device counterpart of _average_band!'s CPU sweep below.

Throws

• ErrorException: no KernelAbstractions extension is loaded, so there is no device kernel to reach (_throw_no_ka_average_kernel).

_threaded_average_engine!(out, in_ref, dims::Tuple, dir::GridDirection, dim_val::Val) -> Nothing

CpuThreaded's _average_engine!: one _average_band! per thread under Threads.@threads :static (every band in turn where a :static loop cannot start, see _static_or_serial), isolated so the closure is never built on a serial path.

_threaded_centered_average_engine!(out, in_ref, dims::Tuple, dim_val::Val) -> Nothing

CpuThreaded's _centered_average_engine!: one _centered_average_band! per thread under Threads.@threads :static, serially where a :static loop cannot start.

backward_average(arg, dim_val::Val)

Constructs or applies the backward averaging operator, representing the operation \(\frac{u_{i-1} + u_{i}}{2}\).

backward_average_dim!(out, in, dims, average_dim)

Low-level, in-place function to compute the backward average of vector in along dimension average_dim, storing the result in out. This function computes \(\frac{u_{i-1} + u_{i}}{2}\).

centered_average(arg, dim_val::Val)

Constructs or applies the centered averaging operator, representing the operation \(\frac{u_{i-1} + 2 u_{i} + u_{i+1}}{4}\).

centered_average_dim!(out, in, dims, average_dim)

Low-level, in-place function to compute the centered average of vector in along dimension average_dim, storing the result in out. This function computes \(\frac{u_{i-1} + 2 u_{i} + u_{i+1}}{4}\).

forward_average(arg, dim_val::Val)

Constructs or applies the forward averaging operator, representing the operation \(\frac{u_{i} + u_{i+1}}{2}\).

forward_average_dim!(out, in, dims, average_dim)

Low-level, in-place function to compute the forward average of vector in along dimension average_dim, storing the result in out. This function computes \(\frac{u_{i} + u_{i+1}}{2}\).

InterpolationNode{D, S, OpType} <: LazyOp{D}

The symbolic interpolation operator: inner_op lives on src_space, and this node evaluates it at points of whatever mesh the assembly is walking.

Side is TrialSide or TestSide, taken from the leaf πₕ wrapped. It decides which leaf _bind_interp_spaces stamps into src_space, whether the stencil names AbsoluteColumns or AbsoluteRows, and which leaf's grid the sweep walks: always the other one, the side that stays native.

src_space is nothing until assembly binds it. πₕ(op) cannot name the space itself: on a composite space the leaf a term interpolates from is only resolved block by block, so _bind_interp_spaces stamps that leaf in at that point, once per block, and the S === Nothing node exists only between construction and that binding.

Distinct from the source wrapper πₕ(uₕ), which carries a grid function's values. This node carries no values; it carries the map, and its stencil names degrees of freedom of the space it interpolates from.

outside is one of :error, :clamp or :extrapolate – never a fill value – and is stored as an isbits OutsidePolicy enum rather than the Symbol given to πₕ, so the node holds no Symbol. A fill value is refused since this node's stencil is a linear map (weighted trial columns), and a constant independent of the trial unknowns cannot be written that way; see interpolation_matrix's docstring for the same restriction.

TrialSide
TestSide

Which side of a bilinear term an InterpolationNode sits on.

A trial-side interpolation names absolute columns and leaves the rows to the walked mesh; a test-side one names absolute rows and leaves the columns to it. The side is a type parameter rather than a field so the node stays a singleton, and so every walk that reads it folds away at compile time.

The walked mesh is the native side's in both cases: it is the one carrying the quadrature weight the product integrates against.

TrialSide
TestSide

Which side of a bilinear term an InterpolationNode sits on.

A trial-side interpolation names absolute columns and leaves the rows to the walked mesh; a test-side one names absolute rows and leaves the columns to it. The side is a type parameter rather than a field so the node stays a singleton, and so every walk that reads it folds away at compile time.

The walked mesh is the native side's in both cases: it is the one carrying the quadrature weight the product integrates against.

_walked_leaf(term, trial_leaf, test_leaf)

The leaf whose grid the assembly sweep walks for term, and whose weights and markers its stencil sees.

The test leaf, unless the term interpolates on the test side: then the rows are named absolutely and the trial leaf is the one that stays native, so it supplies the grid, the quadrature weight and the columns.

backward_finite_difference!(out::AbstractVector, v::GhostVector, dim_val::Val) -> out
backward_finite_difference(v::GhostVector, dim_val::Val) -> Vector

The one-axis operators on a GhostVector: D₋ₓ, D₊ₓ, Mₓ (the backward average), M₊ₓ, Dcₓ, D̃ₓ, D̽ₓ and their ᵧ, ₂ and ! forms, and ∇ₕ(v), the tuple of D₋ along every axis. Each gives the values at the physical nodes, in the mesh's linear order, as a plain vector with no ghosts: those of the same operator on lift(Cₕ, u) at the physical nodes, when v holds the degrees of freedom u and update_ghosts! has filled its ghosts. A stencil reaching beyond a facet with ghosts reads the ghost there, at x₁ - h₂ or xₙ + hₙ; beyond a facet without ghosts it is truncated, as on a plain space.

out is a vector of the physical length, a GhostVector on the mesh of v (its physical part is written) or a grid function on the physical mesh; any other mesh, or an out sharing memory with v, refuses with an ArgumentError. The in-place forms allocate nothing on CpuSerial and CpuPolyester. A result composed with another operator is truncated at the facets, as on a plain space.

See also: Δₕ, update_ghosts!

backward_finite_difference!(out::AbstractVector, v::GhostVector, dim_val::Val) -> out
backward_finite_difference(v::GhostVector, dim_val::Val) -> Vector

The one-axis operators on a GhostVector: D₋ₓ, D₊ₓ, Mₓ (the backward average), M₊ₓ, Dcₓ, D̃ₓ, D̽ₓ and their ᵧ, ₂ and ! forms, and ∇ₕ(v), the tuple of D₋ along every axis. Each gives the values at the physical nodes, in the mesh's linear order, as a plain vector with no ghosts: those of the same operator on lift(Cₕ, u) at the physical nodes, when v holds the degrees of freedom u and update_ghosts! has filled its ghosts. A stencil reaching beyond a facet with ghosts reads the ghost there, at x₁ - h₂ or xₙ + hₙ; beyond a facet without ghosts it is truncated, as on a plain space.

out is a vector of the physical length, a GhostVector on the mesh of v (its physical part is written) or a grid function on the physical mesh; any other mesh, or an out sharing memory with v, refuses with an ArgumentError. The in-place forms allocate nothing on CpuSerial and CpuPolyester. A result composed with another operator is truncated at the facets, as on a plain space.

See also: Δₕ, update_ghosts!

_dot(::DeviceLocality, ::GpuPolicy, u, w::SeparableWeights, v) -> Number

The device counterpart of _dot(u, w::SeparableWeights, v) above: the weight tensor is materialized once via _separable_weights_full, then the product with u and v (reshaped to the grid's own dims, not copied) is one broadcasted reduction.

_dot_masked(::DeviceLocality, ::GpuPolicy, u, w::SeparableWeights, v, mask::BitVector) -> Number

The device counterpart of the BitVector-masked specialization above: mask (always host memory, gpena/Bramble.jl#298) is copied onto u's own device once, reshaped to dims alongside u and v, and the masked sum is one broadcasted reduction.

_dot_masked(::DeviceLocality, ::GpuPolicy, u, w::SeparableWeights, v, mask::MarkedIndicesUnion) -> Number

The multi-marker counterpart: mask is walked once, on the host, into a plain BitVector, then handled exactly as the BitVector method above.

_surface_sum(::DeviceLocality, ::GpuPolicy, Ωₕ, mask, u, v) -> Number

Device counterpart of _surface_sum above: a canonical face (or a union of them) is a set of at most 2D axis-aligned (D-1)-dimensional slices of the grid, so the CPU walk over every grid point – reading half_spacings and u/v one scalar at a time – becomes, per active face, one broadcasted reduction over that slice.

u and v are reshaped to the grid's own dims (a view, not a copy). For each axis d and side (min/max) the mask marks, selectdim takes the boundary slice directly – no mask array is built or walked here, unlike _dot_masked's device path, since a coordinate face is already a contiguous slice rather than an arbitrary marked set. The transverse surface measure for that slice is built once from host_half_spacings (never read on the device), then copied onto u's own device in one transfer via _transverse_weight_device; the reduction itself is sum(conj.(su) .* sv .* w), no @kernel.

Faces are summed independently rather than through one shared mask, exactly as the host _surface_weight_dir recursion does: a corner point on two active faces receives both faces' contributions, which is the intended, non-double-counting definition documented above _face_mask.

_transverse_weight_device(Ωₕ::AbstractMeshType{D}, u, ::Val{d}, ::Val{D}) -> AbstractArray

The (D-1)-dimensional transverse surface measure for the face normal to axis d, as a device array matching u's own backend.

Built the way _separable_weights_full builds the full weight tensor – one reshape per transverse axis so broadcasting multiplies them out to the full transverse shape – except every factor here comes from host_half_spacings rather than a SeparableWeights factor that is already device-resident, since this is a fresh geometric quantity with nothing cached to read. The host tensor is built once, from vectors of length O(n^{1/D}), then copied onto u's device in a single transfer – never read back scalar by scalar the way half_spacing would on the host path.