Grid spaces#

ScalarGridSpace(mesh::MType, weights::SpaceWeights{D, T, VT})
ScalarGridSpace{D, T, VT, MType}(mesh::MType, weights::SpaceWeights{D, T, VT})

Represents a function space for scalar fields defined on a mesh.

A ScalarGridSpace pairs the mesh with the precomputed weight vectors (SpaceWeights) its discrete inner products need.

Fields

• mesh::MType: the underlying mesh of the grid space.

• weights::SpaceWeights{D, T, VT}: precomputed inner product weight vectors, held by reference (SpaceWeights is a mutable struct with const fields).

Discrete inner products

The weights object stores vectors for different discrete \(L^2\) inner products on the space of grid functions. They are defined as follows:

- :innerₕ: The standard discrete \(L^2\) inner product, weighted by the cell measure \(|\square_k|\).

• 1D case:

• 2D case:

• 3D case:

Here, \(|\cdot|\) denotes the measure of the set (length, area, or volume). See cell_measure for details.

- :inner₊, :inner₊ₓ, :inner₊ᵧ, :inner₊₂: Modified discrete \(L^2\) inner products, weighted by a mix of forward/backward spacings (\(h_k\)) and cell widths (\(h_{k+1/2}\)).

• 1D case (:inner₊):

• 2D case (:inner₊ₓ, :inner₊ᵧ):

• 3D case (:inner₊ₓ, :inner₊ᵧ, :inner₊₂):

CompositeGridSpace(spaces::Tuple)
CompositeGridSpace{N}(spaces::Spaces) where {N, Spaces <: Tuple}
CompositeGridSpace{N, Spaces}(spaces::Spaces) where {N, Spaces <: Tuple}

A CompositeGridSpace represents a grid space formed by composing N individual sub-spaces. It is immutable and stack-allocatable, wrapping a tuple of spaces.

Fields

• spaces::Spaces: the tuple of constituent sub-spaces.

VectorGridSpace{N}

Type alias for a CompositeGridSpace{N}, representing a vector-valued function space.

gridspace(Ωₕ::AbstractMeshType{D}) -> ScalarGridSpace{D}
gridspace(Ωₕ::AbstractMeshType{D}; ghosts) -> GhostGridSpace{D}

Constructs a ScalarGridSpace defined on the mesh Ωₕ, precomputing the inner product weights listed in ScalarGridSpace.

With ghosts = :boundary (every boundary facet) or a tuple of facet names such as (:xmin, :ymax), it returns a GhostGridSpace instead: the physical nodes of Ωₕ and one ghost beyond each node of the named facets, with inner products that weigh the physical nodes, cells and facets only. The default ghosts = () is the space above.

Scratch memory is supplied explicitly by callers through in-place mutating operators (such as D₋ₓ!(vₕ, uₕ)), avoiding hidden internal vector buffers.

Weights are a snapshot: they read Ωₕ's grid at the moment of this call and are never refreshed afterward. If Ωₕ is later mutated in place (set_points!, change_points!, iterative_refinement!), the space this call returns keeps the old weights, and its innerₕ/inner₊*/every norm then throw naming the mismatch rather than silently computing against a mesh that no longer exists (gpena/Bramble.jl#221) – call gridspace(Ωₕ) again to get a space that reads the mutated mesh.

Examples

using Bramble
using Bramble: weights
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11))
ndofs(Wₕ) == 11 && sum(weights(Wₕ, Bramble.Innerh())) ≈ 1.0

# output
true

gridspace(Ωₕ::AbstractMeshType, ::Val{N}) where N -> CompositeGridSpace{N}
gridspace(Ωₕ::AbstractMeshType, N::Int) -> CompositeGridSpace{N}

Constructs a vector function space with N components on mesh Ωₕ. The underlying scalar space and its weights are computed once and shared across components.

N == 1 yields the ScalarGridSpace itself rather than a one-component composite, for both spellings. The element interface is uniform either way: uₕ(1) and components(uₕ) work on a scalar element.

The Val form is always type stable. The Int form is stable wherever N is a literal or otherwise constant-foldable, and has a non-concrete return type when N is only known at run time; prefer Val, or Wₕ^N for N ≤ 3, on hot paths.

vector_gridspace(Ωₕ::AbstractMeshType, [N = dim(Ωₕ)]) -> CompositeGridSpace

Convenience constructor for a vector grid space on mesh Ωₕ. If N is omitted, it defaults to the spatial dimension of the mesh (dim(Ωₕ)).

ndofs(Wₕ::AbstractSpaceType) -> Int
ndofs(Wₕ::AbstractSpaceType, ::Type{Tuple}) -> NTuple{N, Int}

Returns the total number of degrees of freedom (DOFs) in the function space Wₕ.

weights(Wₕ::ScalarGridSpace) -> SpaceWeights
weights(Wₕ::ScalarGridSpace, ::Innerh) -> AbstractVector
weights(Wₕ::ScalarGridSpace, ::Innerplus) -> NTuple{D, AbstractVector}
weights(Wₕ::ScalarGridSpace, ::Union{Innerh, Innerplus}, i::Int) -> AbstractVector
weights(Wₕ::ScalarGridSpace{D}, ::Val{S}) where {D, S} -> AbstractVector

Returns the precomputed weight vectors for discrete inner products.

The weights are diagonal matrices (stored as vectors) used in computing discrete \(L^2\) inner products. They represent cell measures or staggered grid spacings.

Every one of them is a SeparableWeights (gpena/Bramble.jl#115, #234): a lazy per-axis product, not a full-grid vector, whatever S is asked for. innerh and each innerplus[d] are built once, at gridspace construction time, and returned unchanged here; the rest are built fresh on each call from the same per-axis factors.

Methods

1. weights(Wₕ) - Returns the full SpaceWeights struct

2. weights(Wₕ, Innerh()) - Returns weights for standard \(L^2\) inner product (cell volumes)

3. weights(Wₕ, Innerplus()) - Returns tuple of weights for modified inner products (all directions)

4. weights(Wₕ, Innerplus(), i) - Returns weights for modified inner product in direction i

5. weights(Wₕ, Innerh(), i) - Same as weights(Wₕ, Innerh()); the cell measures do not depend on a direction, so i is accepted and ignored for interface symmetry

6. weights(Wₕ, Val(S)) - Returns the weights for the staggered set S ⊆ 1:D, whose entry at grid index I is \(\prod_{d \in S} h_d(I_d) \cdot \prod_{d \notin S} h_d(I_d + 1/2)\) (gpena/Bramble.jl#115, #234). S is an NTuple{K, Int} with K ≤ D and any axis order; weights(Wₕ, Val(())) is weights(Wₕ, Innerh()) and weights(Wₕ, Val((d,))) is weights(Wₕ, Innerplus(), d), returning the very same object rather than a recomputed copy. Every other valid S returns a freshly-built SeparableWeights. An S that is not a repeat-free subset of 1:D throws an ArgumentError, as the symbolic inner₊(p, q, Val(S)) does.

Examples

using Bramble: weights
Wₕ = gridspace(Ωₕ)

# Get all weights
w = weights(Wₕ)  # Returns SpaceWeights{D, T, VT}

# Get standard L² weights
w_h = weights(Wₕ, Innerh())  # SeparableWeights, entry I = cell volume at I

# Get modified inner product weights for x-direction
w_plus_x = weights(Wₕ, Innerplus(), 1)  # SeparableWeights for x-direction

# Use in inner product
result = dot(uₕ.data, w_h, vₕ.data)  # Weighted inner product

Defined for a ScalarGridSpace only, the same rule normₕ/norm₊ follow: a composite grid space's leaves can have different meshes and therefore different weights, so there is no single vector that could correctly answer for the whole composite. A CompositeGridSpace raises a MethodError; take a scalar component of it with components first.

Staleness (gpena/Bramble.jl#221)

Every method here funnels through the one-argument form, which checks Wₕ's stored SpaceWeights against mesh(Wₕ)'s current _mesh_version and throws naming the mismatch if an in-place mutator (set_points!, change_points!, iterative_refinement!) has run on the mesh since these weights were built – rather than innerₕ/inner₊*/every norm silently computing against weights for a mesh that no longer exists. Call gridspace again to get a space that reads the mutated mesh. On a CompositeGridSpace this is checked per leaf, the moment innerₕ/etc. recurse into it – there is no separate composite-level check to keep in step.

Form assembly checks once per walk rather than at every point: each walk entry (assemble, an assemble! refill, a matrix-free product, a linear form's assembly, once per leaf of a composite) checks the leaf it walks before its first point, with the same error, and the per-point stencil then reads the stored weights unchecked (gpena/Bramble.jl#437). A bilinear walk entry also checks its trial leaf, which may sit on another mesh (gpena/Bramble.jl#466).

See also: SpaceWeights, SeparableWeights, Innerh, Innerplus, innerₕ

spaces(Wₕ::AbstractSpaceType) -> Tuple

Returns the constituent subspace(s) of Wₕ as a tuple.

space(Wₕ::AbstractSpaceType) -> AbstractSpaceType

Returns the function space Wₕ itself.

space(uₕ::VectorElement) -> AbstractSpaceType

Returns the grid space associated with VectorElement uₕ.

space(sd::Semidiscretization)

Return the test space the semidiscretisation was built on.

ncomponents(Wₕ::AbstractSpaceType) -> Int
ncomponents(::Type{<:AbstractSpaceType}) -> Int

Returns the number of field components of the function space (e.g. 1 for scalar, D for vector).

On a CompositeGridSpace it counts the top-level sub-spaces: a nested composite counts its direct children. leaf_count counts the scalar leaves, and components(uₕ) lists them.

VectorElement(data::VT, space::S)
VectorElement{S, T, VT}(data::VT, space::S)

Represents a grid function (a vector) that belongs to a specific function space.

This is a wrapper that bundles the raw numerical data (the vector data) with its parent space. The space provides the essential context, such as the underlying mesh and associated operators. By subtyping AbstractVector, a VectorElement can be used just like a regular Julia vector in most operations.

Data with offset axes (for example a view over Base.IdentityUnitRange(2:11)) is refused with an ArgumentError: the element indexes its data from 1.

Fields

• data::VT: the raw vector data containing the degrees of freedom.

• space::S: the parent function space to which this vector belongs.

element(Wₕ::AbstractSpaceType) -> VectorElement
element(Wₕ::AbstractSpaceType, α::Number) -> VectorElement

Returns a VectorElement for grid space Wₕ with uninitialized components. If α is provided, the components are initialized to α.

element(Wₕ::AbstractSpaceType, ::Type{T}) -> VectorElement

Returns a VectorElement for grid space Wₕ holding coefficients of type T, with uninitialized components.

The coefficients of a grid function and the coordinates of the mesh under it are two different things, and this is where they part company. element(Wₕ) takes its type from the backend, which is the mesh's own; this takes whatever is asked for, and the container follows through similar, so a Vector backend gives a Vector{T} and a device array gives a device array of T.

The case this exists for is automatic differentiation: a ForwardDiff.Dual grid function over an ordinary Float64 mesh, so that the geometry is not differentiated along with the field. See Rₕ, which uses it to give back an element of whatever type the restricted function returns.

element(Wₕ::AbstractSpaceType, v::AbstractVector) -> VectorElement

Returns a VectorElement for a grid space Wₕ with the same coefficients as v.

parent(uₕ::VectorElement) -> AbstractVector

Returns the coefficient vector containing the degrees of freedom of VectorElement uₕ.

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.

components(uₕ::VectorElement) -> Tuple

Returns an NTuple of VectorElement views, one per leaf of uₕ's space, depth-first, the same leaf a matching-index uₕ(i) returns, regardless of nesting.

components(op::Union{TrialFunction, TestFunction}) -> Tuple
components(op::LazyOp, N::Integer) -> Tuple
components(op::LazyOp, space::AbstractSpaceType) -> Tuple

Returns an NTuple of components of the symbolic trial or test function, suitable for tuple destructuring: (u, v) = components(p).

component_range(Wₕ::CompositeGridSpace, i::Int) -> UnitRange{Int}

Returns the degree-of-freedom index range for the i-th leaf of composite space Wₕ, numbered depth-first (see leaf_spaces_offsets), the same numbering uₕ(i) and dirichlet_components use, so it agrees with them regardless of how deeply Wₕ nests.

component_ranges(Wₕ::CompositeGridSpace) -> NTuple{N, UnitRange{Int}}

N is the number of scalar leaves underneath Wₕ, counting through any nesting.

Returns the degree-of-freedom ranges for every leaf of Wₕ, depth-first; see component_range.

*(x, y...)

Multiplication operator.

Infix x*y*z*... calls this function with all arguments, i.e. *(x, y, z, ...), which by default then calls (x*y) * z * ... starting from the left.

Juxtaposition such as 2pi also calls *(2, pi). Note that this operation has higher precedence than a literal *. Note also that juxtaposition "0x..." (integer zero times a variable whose name starts with x) is forbidden as it clashes with unsigned integer literals: 0x01 isa UInt8.

Note that overflow is possible for most integer types, including the default Int, when multiplying large numbers.

Examples

julia> 2 * 7 * 8
112

julia> *(2, 7, 8)
112

julia> [2 0; 0 3] * [1, 10]  # matrix * vector
2-element Vector{Int64}:
  2
 30

julia> 1/2pi, 1/2*pi  # juxtaposition has higher precedence
(0.15915494309189535, 1.5707963267948966)

julia> x = [1, 2]; x'x  # adjoint vector * vector
5

Rₕ(Wₕ::AbstractSpaceType, f; markers = ()) -> VectorElement

Standard nodal restriction operator. Evaluates f at the grid points of mesh(Wₕ) and returns the result as a VectorElement.

Arguments

• Wₕ::AbstractSpaceType: grid space on which to restrict f.

• f: function of one grid point. It receives a scalar on a 1D mesh and an NTuple{D} on a D-dimensional one, never an SVector.

Keywords

• markers::NTuple{N,Symbol}: restrict evaluation to the named marked regions, leaving every other entry zero.

Examples

Rₕ(Wₕ, x -> sin(x))                # 1D: x is a Float64
Rₕ(Wₕ, x -> sin(x[1]) * x[2])      # 2D: x is a Tuple{Float64,Float64}

# Vector-valued spaces:
Rₕ(Vₕ, (f₁, f₂))                   # one function per component
Rₕ(Vₕ, x -> (f₁(x), f₂(x)))        # one function returning all components

Prefer x -> (f₁(x), f₂(x)) when components share computation, as it evaluates once per grid point on a space whose components share one mesh, whereas (f₁, f₂) always evaluates each component function separately. On a heterogeneous composite (components built over different meshes), the single-function form gives up that advantage, since there is no grid point shared by every component to evaluate it at only once.

See also: Rₕ!, avgₕ.

Rₕ!(uₕ::VectorElement, f; markers = ()) -> VectorElement

In-place version of the restriction operator Rₕ. Evaluates f at the grid points and writes the result into uₕ. Returns uₕ.

Arguments

• uₕ::VectorElement: pre-allocated element to write into.

• f: function of one grid point. It receives a scalar on a 1D mesh and an NTuple{D} on a D-dimensional one, never an SVector.

Keywords

• markers::NTuple{N,Symbol}: restrict evaluation to the named marked regions, leaving every other entry zero. Several markers act as a union.

Examples

Rₕ!(uₕ, x -> sin(x))                  # 1D: x is a Float64
Rₕ!(uₕ, x -> sin(x[1]) * cos(x[2]))   # 2D: x is a Tuple{Float64,Float64}

# only the points carrying the :left marker; the rest stay zero
Rₕ!(uₕ, x -> 1.0; markers = (:left,))

On a Metal (or other device) backend, f runs on the device, evaluated once per grid point inside a compiled kernel (gpena/Bramble.jl#94, #174). It must therefore be GPU-compilable: x -> sin(x[1]) is, but a closure boxing a captured value (a Ref, a mutable struct field, a global variable) or calling a function the device compiler cannot inline is not, and fails at kernel-compile time with a GPUCompiler.InvalidIRError naming the offending call – inherent to GPU execution, not specific to Rₕ!. A masked call (markers non-empty) runs the same kernel over just the marked points, after zeroing uₕ on the device (gpena/Bramble.jl#297).

For an N-component element either shape of f works and both give the same result; the single vector-valued function is evaluated once per grid point when every component shares the same mesh, whereas the tuple always evaluates each component function separately. On a heterogeneous composite, components built over different meshes, the single-function form is instead re-evaluated once per component, since there is no grid point shared by every component to evaluate it at only once:

Rₕ!(uₕ, (f₁, f₂))                     # one function per component
Rₕ!(uₕ, x -> (f₁(x), f₂(x)))          # one function returning all components

A tuple holds one function per leaf component (depth-first on a nested composite); any other length throws a DimensionMismatch.

See also: Rₕ, avgₕ!, element

avgₕ(Wₕ::AbstractSpaceType, f; quad_points = Val(AVG_QUAD_POINTS), markers = ()) -> VectorElement

Returns a VectorElement with the average of function f with respect to the cell_measure of mesh(Wₕ) around each grid point.

Each cell average is a tensor-product Gauss-Legendre rule with quad_points points per direction, exact for polynomials of degree 2 * quad_points - 1.

Arguments

• Wₕ::AbstractSpaceType: grid space on which to average f.

• f: function of one grid point. Receives coordinates as a scalar on 1D meshes or an NTuple{D} on D-dimensional meshes, never an SVector.

Keywords

• quad_points::Union{Integer, Val}: points per direction, per cell. Defaults to Val(AVG_QUAD_POINTS).

• markers::NTuple{N, Symbol}: restrict evaluation to the named marked regions, leaving every other entry zero.

Examples

avgₕ(Wₕ, x -> sin(x))
avgₕ(Wₕ, x -> sin(x[1]) * x[2]; quad_points = Val(4))

On a Metal (or other device) backend, f runs on the device inside the quadrature kernel, evaluated quad_points^D times per cell (gpena/Bramble.jl#94, #174, #693), so it must be GPU-compilable in the same sense Rₕ!'s docstring describes. A masked call runs the same kernel over just the marked cells, after zeroing the result on the device (gpena/Bramble.jl#297).

See also: avgₕ!, Rₕ.

avgₕ!(uₕ::VectorElement, f; quad_points = Val(AVG_QUAD_POINTS), markers = ()) -> VectorElement
avgₕ!(uₕ::VectorElement, f, quad_points::Val) -> VectorElement

In-place version of the averaging operator avgₕ. Returns uₕ.

Evaluates the tensor-product Gauss-Legendre cell average of f and writes the result into uₕ.

Arguments

• uₕ::VectorElement: pre-allocated element to write into.

• f: function of one grid point. Receives coordinates as a scalar on 1D meshes or an NTuple{D} on D-dimensional meshes.

Keywords

• quad_points::Union{Integer, Val}: points per direction, per cell. Defaults to Val(AVG_QUAD_POINTS). Using a Val allows compile-time specialization of the quadrature nodes and weights without boxing.

• markers::NTuple{N, Symbol}: restrict evaluation to the named marked regions, leaving every other entry zero.

On a composite element f may also be a tuple holding one function per leaf component (depth-first on a nested composite); any other length throws a DimensionMismatch.

See also: avgₕ, Rₕ!.

interpolate_at(uₕ::VectorElement, x; outside = :error)

The piecewise (multi)linear interpolant of uₕ at the physical point x, using uₕ's own mesh.

Locates the cell of mesh(space(uₕ)) containing x (locate_cell) and blends the \(2^D\) grid values at that cell's corners, weighted by x's relative position within it (the standard bilinear/trilinear construction), exact for any affine function of the coordinates and correct on a non-uniform mesh, since it reads the mesh's own point coordinates rather than assuming a fixed step.

outside names what happens when x falls outside the domain:

• :error (the default): throw an ArgumentError naming x and the domain extent.

• :clamp: clamp x to the domain first, then interpolate there.

• :extrapolate: blend with the boundary cell's own corner weights and slope, unclamped: a linear extrapolation rather than a hold.

• a Number (e.g. 0.0, NaN): return this value outright, bypassing the blend.

A point outside by no more than a few eps, relative to the domain extent, is always treated as exactly on the boundary, regardless of outside – this is what keeps πₕ between two meshes over the same nominal domain from tripping over floating-point endpoint noise at the rim.

uₕ's mesh is read as it is now, so once an in-place mutator (set_points!, change_points!, iterative_refinement!) has moved it, this throws an ArgumentError (the one weights throws) rather than locating x on a mesh uₕ's values were never computed on; call gridspace again and rebuild uₕ on it. The same check covers πₕ!, πₕ(Wₕ, src) and the symbolic source πₕ(uₕ), which all evaluate through here.

On a GhostGridSpace it interpolates the physical nodes on the physical mesh, the first npoints(Ωₕ) degrees of freedom: the ghost values never enter, and x is outside the domain of Ωₕ, not of the ghosts. A constrained space whose degrees of freedom omit fixed nodes or whose relations carry data has no grid function in them: interpolate_at throws an ArgumentError, and interpolate_at(Bramble.lift(Cₕ, u), x) is its interpolant.

This is the building block both πₕ!/πₕ (below, the numeric operator) and the one-argument, symbolic πₕ use: x -> interpolate_at(uₕ, x) is itself a valid source function, usable anywhere one is accepted, including directly as Rₕ's own argument: πₕ(Wₕ, src) is Rₕ applied to this one function, not a separate mechanism. Rₕ(Wₕ, f) restricts an arbitrary continuous f; when f happens to be another grid function's own interpolant, restricting it is interpolating it, which is why πₕ generalises Rₕ for the case the source is discrete rather than a closed-form function.

See also interpolation_matrix, whose outside accepts only :error, :clamp and :extrapolate – a fill value cannot be represented as a linear map (see the note atop this file).

πₕ!(dest::VectorElement, src::VectorElement; outside = :error) -> VectorElement

Fills dest with the piecewise (multi)linear interpolant of src, sampled at dest's own mesh points, providing the in-place numeric interpolation operator named after the Rₕ!/ avgₕ! convention.

interpolate_at(src, ·) is a genuine function of a physical point (evaluable anywhere, not only at src's own grid points), so this is equivalent to Rₕ!(dest, x -> interpolate_at(src, x; outside)). dest and src may be built over entirely different meshes. Rₕ! handles evaluating the interpolant at each of dest's grid points, following dest's backend execution_policy. outside is forwarded to interpolate_at unchanged; see its docstring for the policy this decides among.

Every call here re-locates, via locate_cell, which cell of src's mesh each of dest's points falls in. Interpolating repeatedly between the same two meshes (a time loop transferring a coefficient between two composite leaves, say) should build that once instead: see the interpolation_matrix-based method below, following the same "build the pattern once" shape allocate_system_matrix/assemble! already use.

When dest and src are both host-backed, this evaluates interpolate_at pointwise as described above. When either is device-backed, the pointwise path would compile a device kernel around the host closure x -> interpolate_at(src, x) – which cannot compile, since interpolate_at dispatches, calls locate_cell and reads src's own coefficients on the host. Instead, this assembles interpolation_matrix(space(dest), space(src); outside) once and runs it through the mul!-based method below, which uploads P and src next to dest as needed – so a fill value (outside::Number, meaningless as a linear map) throws the same ArgumentError interpolation_matrix itself throws, naming that function.

πₕ!(dest::VectorElement, P::AbstractMatrix, src::VectorElement) -> VectorElement

Fills dest via a precomputed interpolation_matrix P instead of re-locating each destination point's cell: parent(dest) .= P * parent(src), computed in place via mul!.

P must be interpolation_matrix(space(dest), space(src)) (or an equal-shape matrix built the same way) – a mismatched size throws the usual DimensionMismatch from mul!. interpolation_matrix always returns a host SparseMatrixCSC, regardless of dest's own backend: passing it straight to a device-backed dest/src would otherwise fall through SparseArrays' own generic sparse-times-vector method, which scalar-indexes the device vectors one entry at a time. So P and src are each brought to dest's own locality first – P uploaded with metal_sparse_csr if dest is device-backed and P is not already, src copied wholesale (never scalar-read) if its locality disagrees with dest's – and only then does mul! run, zero-allocation once P already lives next to dest and src (the loop in πₕ!'s own docstring above).

Repeated interpolation between the same two meshes should build P once and reuse it here every subsequent call, exactly as allocate_system_matrix/assemble! split the sparsity pattern (expensive, built once) from refilling values (cheap, every step):

using Bramble: interpolation_matrix
P = interpolation_matrix(space(dest), space(src))
for step in 1:nsteps
    Rₕ!(src, coefficient_at(step))
    πₕ!(dest, P, src)   # no locate_cell search, zero allocations
end

πₕ(Wₕ::ScalarGridSpace, src::VectorElement; outside = :error) -> VectorElement

Evaluates Rₕ(Wₕ, x -> interpolate_at(src, x; outside)); see πₕ!. Distinguished by multiple dispatch from the one-argument symbolic wrapper πₕ(uₕ), defined below. The element type is promoted from Wₕ's and src's own, so interpolating a Dual-valued src yields a Dual-valued result on an undifferentiated Wₕ. outside is forwarded to interpolate_at unchanged; a fill value joins that promotion and is converted to the result's element type, whichever grid points fall outside src's domain.

πₕ(uₕ::VectorElement; outside = :error) -> LazyOp

The interpolant of uₕ, as a symbolic source term; usable anywhere a source is, including inside another operator: innerₕ(D₋ₓ(πₕ(uₕ)), D₋ₓ(v)) differentiates the interpolated field the same way D₋ₓ differentiates any other source, innerₕ(Mₓ(πₕ(uₕ)), v) averages it, and so on. This enables a coupled form to evaluate a leaf's grid function on a different leaf's mesh.

Built as a SourceFunction whose function evaluates interpolate_at(uₕ, x; outside): a SourceFunction's own local_stencil evaluates its function at the current point of whichever mesh is being walked, so uₕ can originate from another leaf without special handling; the interpolation occurs once per point inside interpolate_at, where ordinary source function calls occur. outside is forwarded to interpolate_at unchanged, fill values included – this is a source (a function of uₕ's values), not a linear map, so unlike the operator πₕ(op) below it carries no such restriction.

Unlike a plain function source, a form does not sample this one once when it is built: every assemble/assemble! reads uₕ's current values, and throws an ArgumentError once uₕ's mesh has moved (see interpolate_at).

πₕ(op::LazyOp{D}; outside = :error) -> InterpolationNode

The interpolation operator onto whichever mesh the form integrates over, applied to the trial or the test function op. The space interpolated from is that function's own, and is not written here: it is bound during assembly, once the concrete leaf is known (_bind_interp_spaces), which is also what makes this work on a composite space, where the leaf a term names is only resolved block by block.

This is the bilinear counterpart of πₕ(uₕ): that one interpolates a grid function whose values are already known, and belongs on the source side of a linear form; this one interpolates the unknown, and so contributes matrix columns. innerₕ(πₕ(u), v) assembles Hᵥ · P, with P the same matrix interpolation_matrix builds: computed a row at a time during the sweep rather than as a matrix product, which allows assemble! to refill it with zero allocations.

Operators wrap it from the outside, acting on the mesh being integrated over: inner₊(D₋ₓ(πₕ(u)), D₋ₓ(v)) is D_x^⊤ H_+ D_x P. Writing an operator inside is a different operation and is refused, since it would difference on the source mesh instead.

innerₕ(u, πₕ(w)) is the mirror: the test side interpolates, so the trial mesh is the one integrated over and the entries name absolute rows instead of absolute columns. The assembled matrix is Pᵀ · H · (trial factor). A single term interpolating both sides is refused – one side has to stay native, since it is the side whose mesh carries the quadrature weight.

op must be a trial- or test-function leaf, plain or indexed. outside accepts only :error (the default), :clamp and :extrapolate – see InterpolationNode's own docstring for why a fill value is refused here.

On a GhostGridSpace the blend reads the physical nodes, as interpolate_at does, and a ghost point lies outside the domain. Where the product's weight is zero, as at every ghost point in innerₕ(πₕ(u), v), the value is discarded and outside is not consulted, so the default :error assembles there. It throws where a ghost value enters a row: inner_Γ(∂ₙ(πₕ(u)), v) and D₋ₓ(πₕ(u)) read πₕ at a ghost point from a physical node. Pass outside = :clamp or outside = :extrapolate there.

interpolation_matrix(Wdest::ScalarGridSpace, Wsrc::ScalarGridSpace; outside = :error) -> SparseMatrixCSC

The piecewise (multi)linear interpolant of πₕ/interpolate_at as a sparse matrix P rather than applied pointwise: P * parent(src) ≈ parent(πₕ(Wdest, src)) for any src::VectorElement over Wsrc. P is ndofs(Wdest) × ndofs(Wsrc), generally rectangular (since Wdest and Wsrc are built over different meshes), with at most \(2^D\) nonzero entries per row: the corner weights of the source cell locate_cell places that destination point in.

outside accepts only :error (the default), :clamp, and :extrapolate – not a fill value, unlike interpolate_at: P is a linear map, and a row that returns a constant regardless of src cannot be written as a weighted combination of src's own entries unless that constant is exactly zero. Passing a Number throws, naming this. Agrees entry-for-entry with pointwise interpolate_at under each of the three policies it does accept, except far outside the domain under :extrapolate: once a cell fraction t passes 1/eps, 1 - t rounds to -t and P's corner weights cancel, while interpolate_at blends differences there and keeps a constant constant.

Unlike D₋ₓ(Wₕ) and the other operator matrices, this is always a SparseMatrixCSC, regardless of either space's own backend matrix_type. Those matrices are built from shift (a fixed diagonal offset generalized via Kronecker products), but which source cell a destination point falls in has no such regular structure across two independent meshes: it is genuinely sparse and irregular, assembled directly from locate_cell rather than composed from a handful of shifts. Converting the result to another matrix type, where that is meaningful, is left to the caller.

GhostGridSpace{D, T, S, ET, DT, TD, GU} <: AbstractSpaceType{1}

A scalar grid space with ghost points beyond the named facets of its mesh, built by gridspace(Ωₕ; ghosts = :boundary) or gridspace(Ωₕ; ghosts = (:xmin, ...)).

Its degrees of freedom are the physical nodes of Ωₕ, in the mesh's linear order, then the ghosts in the order of GhostMesh. padded_space(Wₕ) is the ScalarGridSpace on the padded mesh, the one forms assemble on: its innerₕ, inner₊ and inner_Γ weigh the physical nodes, cells and facets only, and are 0 at every ghost and absent point, so they equal those of gridspace(Ωₕ) on the physical nodes. closure(Wₕ) is the pair (E, d) that maps the degrees of freedom u to the padded values E u + d; on a space whose relations have data depending on time (see constrained), d depends on time too.

The weights follow the version of the mesh they were built from: after change_points! of Ωₕ the ghost coordinates follow, and the padded space's weights throw as those of any stale space; call gridspace(Ωₕ; ghosts = ...) again.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 2.0)), (5, 7), (false, false))
Wₕ = gridspace(Ωₕ; ghosts = :boundary)
ndofs(Wₕ)                     # 35 physical nodes + 24 ghosts
E, d = Bramble.closure(Wₕ)    # E is 63 × 59, d = 0

constrained(Wₕ, relations) -> GhostGridSpace

The space Wₕ, built by gridspace(Ωₕ; ghosts = ...), with the ghosts of some facets eliminated through boundary relations. relations is a function of the unknown u that returns one pair facet => (operator => data) per related facet, or a tuple of them:

Cₕ = constrained(Wₕ, u -> (:xmin => ∂ₙ(u) => 0.0,
    :xmax => 2.0 * u + 0.5 * ∂ₙ(u) => x -> exp(x[1])))

The operator is either u itself, a Dirichlet relation, or a facet operator that reads the facet's ghost: ∂ₙ (conormal ∂ₙ(u; K) included), ∂ₙₙ, and their linear combinations with u. The data is a number, a function of the point x -> …, or a function of the point and the time (x, t) -> … (told apart by the methods it has). At each node of the facet, a ghost relation is solved for the node's ghost; at a 2D corner where two conormal relations read each other's ghost, the two are solved together. A Dirichlet relation :F => u => g removes the nodes of F from the unknowns and fixes them to g; its facet needs no ghosts, and if it has them it needs one ghost relation besides, which solves for them. A ghosted facet without any relation keeps its ghosts as unknowns.

The degrees of freedom of the result are the physical nodes no Dirichlet relation fixes, in the mesh's linear order, then the ghosts that remain unknowns, in the order of Wₕ. Its closure (E, d) maps them to the padded values E u + d, where d carries the data. The load of a pair whose trial space has d ≠ 0 reads the bilinear form, so assemble(a, l) builds the two together, its matrix Eᵀ A E symmetric when A is; assemble(l) alone refuses a test space with d ≠ 0, which is exact when both sides are one space. Dirichlet conditions can instead replace rows, assemble(a, l; dirichlet, symmetrize), on a space whose ghosts are all eliminated and none of whose nodes are fixed by u => g: its unknowns are then the physical nodes, as on gridspace(Ωₕ).

When some data depend on time, E stays the same and only d depends on time, d(t). lift(Cₕ, u, t) gives E u + d(t); lift(Cₕ, u), lift(Cₕ) and assemble(a, l) refuse with an ArgumentError, as there is no single d, while form and assemble(a), which read E only, work as on a static space. A time-dependent problem on such a space is posed through semidiscretize. constrained evaluates the data at t = 0, so they must be defined there, and closure(Cₕ) gives d(0) as its d. On a composite of ghost spaces, constrained(C₁ × C₂, relations) relates the components together, a relation reading several of them (a traction) included: see constrained(::CompositeGridSpace, relations).

Dirichlet data in time that meet at a corner are checked to agree at t = 0 only; the corner then takes the first facet's datum at every t. Under CpuPolyester a datum that closes over an array crosses each task whole and allocates on every evaluation; one that closes over numbers only, or over nothing, allocates nothing.

Refused with an ArgumentError: a space without ghosts, an already constrained space, a facet with two Dirichlet or two ghost relations, a ghost relation on a facet that carries no ghosts, a ghosted facet with u => g and no ghost relation, a relation that does not read its facet's ghost (β = 0, or a jumpₙ alone), two Dirichlet relations disagreeing at a shared node, and a singular block at a corner (a singular K).

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0)), 11, false)
Wₕ = gridspace(Ωₕ; ghosts = :boundary)
Cₕ = constrained(Wₕ, u -> (:xmin => ∂ₙ(u) => 1.0, :xmax => u + ∂ₙ(u) => 2.0))
a = form(Cₕ, Cₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
l = form(Cₕ, v -> innerₕ(x -> 1.0, v))
A, F = assemble(a, l)

Dₕ = constrained(gridspace(Ωₕ; ghosts = (:xmax,)), u -> (:xmin => u => 0.0,
    :xmax => ∂ₙ(u) => x -> x[1]))  # 10 unknowns: the node at xmin is fixed

Tₕ = constrained(Wₕ, u -> (:xmin => ∂ₙ(u) => (x, t) -> sin(t), :xmax => u + ∂ₙ(u) => 2.0))
Bramble.lift(Tₕ, zeros(Bramble.ndofs(Tₕ)), 0.5)   # E u + d(0.5)

See also: gridspace, ∂ₙ, ∂ₙₙ

constrained(Vₕ::CompositeGridSpace, relations) -> ConstrainedCompositeSpace

A composite of ghost spaces (C₁ × C₂, each built by gridspace(Ωₕ; ghosts = ...) with the same ghosts) constrained as a whole. relations(u) returns the pairs facet => operator => data of constrained on a ghost space, written of the components u(i). A Dirichlet relation names its component: :ymin => u(2) => 0.0 fixes the nodes of :ymin in component 2, at most one per component and facet. A ghost relation may read one component, :xmin => ∂ₙ(u(1)) => g, or several, as the traction of component i, Σⱼ ∂ₙ(u(j); K = Kᵢⱼ), does; its tangential terms are the conormal's. A facet's ghost relations are not keyed by component. Relations of one component each take at most one per component, and a component without one keeps its ghosts on that facet as unknowns, as on the product; a facet with a relation reading several components has one ghost relation per component. Each is attached to the ghosts of one component (the one it reads alone, if any), and the relations of a node, or of a corner, are solved together, so a traction closes jointly with its partner.

The relations of all the components are solved together on the padded composite, and the result stores the closure (E, d) they give (Bramble.closure), whose E couples the components wherever a relation does. Its degrees of freedom are component 1's, then component 2's, …, each in the order constrained gives on the component alone, so with relations of one component each it equals constrained(C₁, r₁) × constrained(C₂, r₂). Forms on it, assemble(a, l) and lift work as on that product.

A displacement facet beside a traction facet carries ghosts, which the traction's tangential difference reads at their corner: close them by ∂ₙₙ(u(i)) => 0 besides u(i) => g. At a corner of two facets that each have a relation reading several components (two tractions), the tangential difference of each facet's relations is one-sided into the domain, second order over the three nodes along the facet, so each facet's relations read its own ghosts alone and the corner closes facet by facet: solved together, isotropic tractions would give a singular block, their shear relations being one equation there (the stress is symmetric). Elsewhere the difference stays centred.

A form on it reads the ghosts of a facet with a relation reading several components through second differences as the closure gives them, and, in the left argument of a product, through a first difference along one axis of a first difference along another (under sums and scales only, such as Dcₓ(Dcᵧ(u(2))) or Dcₓ(2.0 * Dcᵧ(u(2)))) as the quadratic through the three nodes inward along the facet's normal, extrapolated to the ghost: the difference across the facet is then one-sided, second order on the non-uniform spacing. The right argument reads the closure's ghosts; an average or shift between the two differences has no facet rule and reads past a corner of the padded grid, so a form with one is refused with an ArgumentError. Read through the traction's tangential difference, the ghosts made the collocated elasticity operator unstable on cells longer across the facet than along it. Relations of one component each leave every form as on the product.

Refused with an ArgumentError, besides what constrained and the product refuse: a relation naming no component, two ghost relations of one component on a facet, a facet with a relation reading several components but not one ghost relation per component, relations whose block at a node is singular, a corner of two tractions with fewer than three nodes along a facet, a facet with a relation reading several components and fewer than three nodes across it (refused by form), data depending on time (constrain the components apart and take the product), the dirichlet keyword of assemble, dirichlet_bc! and symmetrize! (give u(i) => g instead), and GhostVector and update_ghosts!.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (7, 9), (false, false))
Sₕ = gridspace(Ωₕ; ghosts = (:xmin, :xmax))
Vₕ = constrained(Sₕ × Sₕ, u -> (:xmin => ∂ₙ(u(1)) => 1.0, :xmax => ∂ₙ(u(1)) => 0.0,
    :xmin => ∂ₙ(u(2)) => 0.0, :xmax => u(2) + ∂ₙ(u(2)) => 0.0))
E, d = Bramble.closure(Vₕ)

# isotropic traction on :xmin, displacement on the other facets
λ, μ = 2.0, 0.5
K(i, j) = [λ * (i == d) * (j == e) + μ * ((i == j) * (d == e) + (i == e) * (j == d))
           for d in 1:2, e in 1:2]
t(u, i) = ∂ₙ(u(1); K = K(i, 1)) + ∂ₙ(u(2); K = K(i, 2))
disp(F) = u -> (F => u(1) => 0.0, F => u(2) => 0.0, F => ∂ₙₙ(u(1)) => 0.0,
    F => ∂ₙₙ(u(2)) => 0.0)
Wₕ = gridspace(Ωₕ; ghosts = :boundary)
Tₕ = constrained(Wₕ × Wₕ, u -> (:xmin => t(u, 1) => 1.0, :xmin => t(u, 2) => 0.0,
    disp(:ymin)(u)..., disp(:xmax)(u)..., disp(:ymax)(u)...))

lift(Wₕ::GhostGridSpace) -> VectorElement
lift(Wₕ::GhostGridSpace, u::AbstractVector) -> VectorElement
lift(Wₕ::GhostGridSpace, u::AbstractVector, t::Number) -> VectorElement

Grid functions on the padded space of Wₕ (Bramble.padded_space(Wₕ)) from its closure (E, d) (Bramble.closure(Wₕ)).

lift(Wₕ) is the data d: on a space built by constrained, the ghost values the relations' data give when every degree of freedom is 0, and 0 everywhere else. lift(Wₕ, u) is E u + d, the grid function of the degrees of freedom u (a solution of a system assembled on Wₕ) with its eliminated ghosts filled in by the relations. The runtime operators read the ghosts of the result, and both are data in forms on Wₕ, so a linear form can name, say, the full Δₕ of the lifting at the facet nodes. Absent padded points (corners beyond two ghosted facets) hold 0. After change_points! of the mesh, lift throws as the space's weights do: build the space again.

lift(Wₕ, u, t) is E u + d(t), for a space whose relation data depend on time, (x, t) -> … (see constrained): d(t) comes from the data at t, and E is the same at every t. On a static space it is lift(Wₕ, u). It evaluates the data into a buffer of its own, so it may be called concurrently on one space, and a time of a non-float type, such as a ForwardDiff.Dual, gives values of that type. On a space with time-dependent data, lift(Wₕ) and lift(Wₕ, u) refuse with an ArgumentError, as there is no single d: give the time, lift(Wₕ, u, t), or pose the problem in time through semidiscretize.

lift(Wₕ, u) runs in the bands of the mesh's execution policy, each padded point summed by one band, so it gives the same values on every CPU policy. Bramble.lift!(uₕ, Wₕ, u) and Bramble.lift!(uₕ, Wₕ, u, t) refill a grid function lift returned, allocating nothing.

lift is public, not exported: write Bramble.lift or using Bramble: lift.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0)), 9, false)
Cₕ = constrained(gridspace(Ωₕ; ghosts = :boundary),
                 u -> (:xmin => ∂ₙ(u) => 1.0, :xmax => u + ∂ₙ(u) => 2.0))
a = form(Cₕ, Cₕ, (u, v) -> innerₕ(-1.0 * Bramble.D̃ₓ(Bramble.D₋ₓ(u)), v))
l = form(Cₕ, v -> innerₕ(x -> 1.0, v))
A, F = assemble(a, l)
uₕ = Bramble.lift(Cₕ, A \ F)          # 11 padded values, ghosts included
Δₕ(uₕ)                                 # reads the ghosts at the facet nodes

Dₜ = constrained(gridspace(Ωₕ; ghosts = :boundary),
                 u -> (:xmin => ∂ₙ(u) => (x, t) -> sin(t), :xmax => ∂ₙ(u) => 0.0))
vₕ = Bramble.lift(Dₜ, zeros(Bramble.ndofs(Dₜ)), 0.5)   # the ghosts of the data at t = 0.5

lift(Vₕ::CompositeGridSpace) -> VectorElement
lift(Vₕ::CompositeGridSpace, u::AbstractVector) -> VectorElement
lift(Vₕ::CompositeGridSpace, u::AbstractVector, t::Number) -> VectorElement

lift on a composite of ghost spaces (C₁ × C₂): the grid function on Bramble.padded_space(Vₕ) whose component i is lift(Cᵢ, uᵢ) (d, E u + d or E u + d(t)) for the slice uᵢ of u at component_range(Vₕ, i). Each component is checked, and refused, as lift on it is. Bramble.lift!(uₕ, Vₕ, u) and Bramble.lift!(uₕ, Vₕ, u, t) refill a grid function lift returned.

On a composite constrained as a whole, constrained(C₁ × C₂, relations), they give d and E u + d from its stored closure over the whole padded composite, as a relation reading several components (a traction) makes a component's ghosts depend on the others' degrees of freedom; its data do not depend on time, so t changes nothing.

lift!(uₕ::VectorElement, Wₕ::GhostGridSpace, u::AbstractVector) -> uₕ
lift!(uₕ::VectorElement, Wₕ::GhostGridSpace, u::AbstractVector, t::Number[, g]) -> uₕ

Write E u + d into the grid function uₕ, the in-place form of lift(Wₕ, u): uₕ lives on the padded space of Wₕ, as the result of lift does, so a solve loop lifts each new solution without allocating. The values are those of lift(Wₕ, u), bit for bit, on every CPU policy. With t, it writes E u + d(t), the values of lift(Wₕ, u, t); on a space whose relation data depend on time, the form without t refuses, as lift(Wₕ, u) does. The data at t go into a buffer the space stores, so lift!(uₕ, Wₕ, u, t) is not safe to call concurrently on one space; each concurrent caller passes its own buffer g = Bramble._affine_buffer(Wₕ) instead.

lift! refuses a grid function on another space, or of an element type that cannot hold the values (a real uₕ for a complex u), with an ArgumentError. It is public, not exported: write Bramble.lift!.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0)), 9, false)
Cₕ = constrained(gridspace(Ωₕ; ghosts = :boundary),
                 u -> (:xmin => ∂ₙ(u) => 1.0, :xmax => u + ∂ₙ(u) => 2.0))
a = form(Cₕ, Cₕ, (u, v) -> innerₕ(-1.0 * Bramble.D̃ₓ(Bramble.D₋ₓ(u)), v))
l = form(Cₕ, v -> innerₕ(x -> 1.0, v))
A, F = assemble(a, l)
uₕ = Bramble.lift(Cₕ, A \ F)
assemble!(A, F, a, l)                  # the same system, refilled in place
Bramble.lift!(uₕ, Cₕ, A \ F)

GhostVector{T, VP, VG, G} <: AbstractVector{T}

The values of a grid function on a ghost space: the physical values, the very vector the solver holds, and the values at the ghosts beyond the named facets, kept in a buffer of their own.

GhostVector(Wₕ, u) takes the vector u of the values at the physical nodes of the mesh of Wₕ, in the mesh's linear order, and holds u itself as its physical part, not a copy: a time loop that writes its state into u sees it in the GhostVector. The ghost buffer has one entry per ghost of mesh(Wₕ), in the ghost order of GhostMesh (the facets xmin, xmax, ymin, …, each facet's nodes in the mesh's column-major order), and starts at 0; update_ghosts! fills it from the relations of a constrained space. GhostVector(Wₕ) starts from zeros, its physical part a vector of the backend of Wₕ. Wₕ is any space with ghosts, as built by gridspace(Ωₕ; ghosts = ...) or constrained.

A GhostVector is an AbstractVector of length np + ng: the np physical values, then the ng ghosts, the order of the degrees of freedom of the unconstrained ghost space gridspace(Ωₕ; ghosts = ...). Bramble.physical(v) and Bramble.ghost_values(v) return the two parts. It carries the mesh with ghosts, mesh(Wₕ), and update_ghosts! refuses a space on another mesh.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 2.0)), (5, 7), (false, false))
Cₕ = constrained(gridspace(Ωₕ; ghosts = :boundary),
                 u -> (:xmin => ∂ₙ(u) => 1.0, :xmax => u + ∂ₙ(u) => 0.0,
                       :ymin => ∂ₙ(u) => 0.0, :ymax => ∂ₙ(u) => 0.0))
u = rand(35)                    # the values at the 5 × 7 physical nodes
v = GhostVector(Cₕ, u)          # Bramble.physical(v) === u
update_ghosts!(v, Cₕ)           # the 24 ghosts the relations give
length(v)                       # 59

update_ghosts!(v::GhostVector, Cₕ::GhostGridSpace) -> v
update_ghosts!(v::GhostVector, Cₕ::GhostGridSpace, t::Number[, g]) -> v

Write into the ghost buffer of v the ghosts that the relations of Cₕ eliminate, from the values in v: each is its row of the closure (E, d) of Cₕ (see lift), applied to the physical values of v and to the ghosts Cₕ keeps as unknowns, which are read from the buffer. These are the ghost values of lift(Cₕ, u) for the degrees of freedom u that v holds, the conormal corner blocks included; a point a Dirichlet relation fixes contributes its datum, whatever v holds there. The ghosts of a facet without a relation, and the physical part, are left as they are. On gridspace(Ωₕ; ghosts = ...) itself, with no relations, nothing is written.

The update reads the eliminated ghosts' rows only, stored with the space, so it costs a sweep over facet-sized data, allocates nothing, and runs in the bands of the mesh's execution policy, each ghost summed by one band: every CPU policy gives the same values.

With t, the relation data are taken at time t (data (x, t) -> …, see constrained): the ghosts of lift(Cₕ, u, t). On a static space the time is ignored; on a space whose data depend on time, the form without t refuses with an ArgumentError. The data at t go into a buffer the space stores, so update_ghosts!(v, Cₕ, t) is not safe to call concurrently on one space; each concurrent caller passes its own buffer g = Bramble._affine_buffer(Cₕ).

update_ghosts! refuses with an ArgumentError a GhostVector on the mesh of another space, one whose parts were resized, one whose element type cannot hold the ghost values (a real v on a complex space), and a space on a mesh changed in place since it was built.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0)), 9, false)
Cₕ = constrained(gridspace(Ωₕ; ghosts = :boundary),
                 u -> (:xmin => ∂ₙ(u) => (x, t) -> sin(t), :xmax => u + ∂ₙ(u) => 0.0))
v = GhostVector(Cₕ, collect(points(Ωₕ)))
update_ghosts!(v, Cₕ, 0.5)
Bramble.ghost_values(v)         # the ghosts at xmin and xmax at t = 0.5

physical(v::GhostVector) -> AbstractVector

The values of v at the physical nodes, in the mesh's linear order: the vector v was built from, not a copy. physical is public, not exported: write Bramble.physical.

ghost_values(v::GhostVector) -> AbstractVector

The values of v at the ghosts, one per ghost of its mesh in the ghost order of GhostMesh: the buffer update_ghosts! writes, not a copy. ghost_values is public, not exported: write Bramble.ghost_values.

assemble!(A, F, a::GhostBilinearForm, l::GhostLinearForm;
    dirichlet = nothing, dirichlet_components = nothing, symmetrize = false)

Refill A and F with the matrix and the load that A, F = assemble(a, l) returns, and return A, F. A keeps the pattern and F the element type that assemble(a, l) gave, so call that once and refill its two outputs. The call allocates nothing warm.

Both forms are assembled together because, when the closure of the trial space has an affine part d ≠ 0, the load is E_vᵀ(F − A d) and reads the matrix. assemble!(b, l) alone refuses a test space with d ≠ 0, which is exact when both sides are one space: on a pair whose trial space has data and whose test space has none, it gives E_vᵀ F without the lift. dirichlet, dirichlet_components and symmetrize act as in assemble(a, l).

It throws an ArgumentError when the trial space of a or the test space of l is time-dependent (pose the problem with semidiscretize), for symmetrize = true without dirichlet, for A that is not a SparseMatrixCSC with the reduced pattern of a, for a dirichlet the rows of the space do not allow, and for an F whose element type cannot hold the load's. A DimensionMismatch follows an F of the wrong length and forms whose test spaces differ: not the same space, and meshes with other points or closures (E, d) that are not equal.

polygon_spaces(Ωₕ::MeshnD{2}, label::Symbol, g) -> (U, W)

The trial space U and the test space W of the scheme of Ferreira and Grigorieff (2006) on the region marked label on Ωₕ, a mesh of a domain with the marker label => polygon(vertices) or label => levelset(φ), with the Dirichlet datum g(x) on the region's boundary nodes. For a level-set marker the region is the polygon that mesh inscribes in the curve \(φ = 0\) (its vertices are the boundary nodes), and g is read at those vertices, not projected onto the curve: a vertex lies off the curve by \(O(h^2)\).

Both are GhostGridSpaces over Ωₕ, whose degrees of freedom are the interior nodes of the polygon (the label label), in the mesh's linear order. lift(U, u) is the grid function on all of Ωₕ: u at the interior nodes, g at the boundary nodes (Symbol(label, :_boundary)), the extension (3.2) of the paper at each exterior diagonal neighbour p = k + (sx, sy) of an interior node k, u_p = u_{k+(sx,0)} + u_{k+(0,sy)} − u_k, and 0 at every other node. W selects the interior nodes. Assemble with form(U, W, f) and form(W, f):

U, W = Bramble.polygon_spaces(Ωₕ, :L, g)
A, F = assemble(form(U, W, (u, v) -> innerₕ(-1.0 * Δₕ(u), v)),
                form(W, v -> innerₕ(f, v)))
uₕ = Bramble.lift(U, A \ F)

A form may read, in the equation of an interior node, its 5-point neighbours and its diagonal neighbours (Δₕ, Dcₓ(Dcᵧ(u)), first differences): an exterior node beyond them has no value and the form throws. An exterior node two interior nodes extend differently (a V-notch) is substituted per reading equation in the reduction.

Rₕ(U, f) is the element of U holding f at the interior nodes, and normₕ, snorm₁ₕ, norm₁ₕ and norminf of an element of U measure the grid function on Ω̄ₕ that is the element at the interior nodes and 0 at the boundary nodes, so that the error of a solve on the region is

e = element(U, A \ F) - Rₕ(U, u)
normₕ(e), snorm₁ₕ(e), norminf(e)

W refuses Rₕ and the norms.

U × U and U^Val(n) of the trial spaces, tested by W × W, assemble per component, each with its own g. The components must be built on the same mesh Ωₕ. constrained on such a composite throws, as g fixes each component's boundary.

Throws

• ArgumentError: if Ωₕ has no polygon or level-set marker label; if the grid does not fit the region (a cell holding interior and exterior nodes other than across a diagonal between two boundary nodes, as after change_points! to lines the closure would not give); if the region has no interior node. form(U, U, f) and form(U, f) throw, as U is no test space, and so do dirichlet = in assemble and constrained on a composite of these spaces: g fixes the boundary.