Reference

abstract type AbstractMeshField{N,T,V}

A mutable, node-centered discrete field on a CartesianGrid: it attaches a value to every node of the grid and can be both read and written by Cartesian multi-index.

The type parameters expose the field's basic shape: N is the spatial dimension, T the grid's coordinate type, and V the type of the stored values (e.g. Float64 for a level set, SVector{N,Float64} for a velocity field).

Interface

A concrete subtype must implement:

• mesh(ϕ): the underlying CartesianGrid; every geometric quantity (ndims, axes, meshsize, …) derives from it.
• ϕ[I] (Base.getindex): the value at node I. This must also return a value for indices just outside the grid, so that finite-difference stencils work near the boundary.
• ϕ[I] = v (Base.setindex!): write v at node I.

The following optional methods can be implemented for efficiency:

• active_nodeindices(ϕ): node indices the field is actually defined on; defaults to nodeindices(ϕ) (every grid node).
• active_cellindices(ϕ): likewise for cells; defaults to cellindices(ϕ).
• copy(ϕ): clone the field; defaults to deepcopy. The time integrators allocate stage buffers with it, so overriding it to share the immutable mesh/bcs is cheaper.
• copy!(dest, src): overwrite dest's values (and active set) with src's; a generic fallback copies through values. The time integrators recycle stage buffers with it.

See MeshField (dense) and NarrowBandMeshField (sparse band) for examples.

AdvectionTerm(𝐮[, scheme = WENO5(), update_func = nothing])

Advection term representing 𝐮 ⋅ ∇ϕ. Available schemes are Upwind and WENO5.

If passed, update_func will be called as update_func(𝐮, ϕ, t) before computing the term at each stage of the time evolution. This can be used to update the velocity field 𝐮 depending not only on t, but also on the current level set ϕ.

struct BernsteinPolynomial{N, T}

A multidimensional Bernstein polynomial with coefficients of type T in N dimensions, defined on a hyperrectangle [low_corner[i], high_corner[i]].

BernsteinPolynomial(c::AbstractArray, lc, hc)

Create a multidimensional Bernstein polynomial with coefficients c defined on the hyperrectangle [lc[1], hc[1]] × … × [lc[N], hc[N]].

Calling p(x) evaluates the polynomial at the point x = (x[1], …, x[N]) by the formula

\[p(x_1,\dots,x_D)=\sum_{i_j=0}^{d_j}c_{i_1\dots i_D}\prod_{j=1}^D\binom{d_j}{i_j}(x_j-l_j)^{i_j}(r_j-x_j)^{d_j-i_j}\]

where $l_j = lc[j]$ and $r_j = hc[j]$ are the lower and upper bounds of the hyperrectangle, respectively, and $d_j = size(c)[j] - 1$ is the degree of the polynomial in dimension j.

BernsteinPolynomials can be differentiated using ForwardDiff; see gradient and hessian.

abstract type BoundaryCondition

Abstract supertype for boundary conditions. A boundary condition determines the value of a field at ghost indices that lie outside the mesh, closing finite-difference stencils near the domain boundary.

A concrete subtype bc must implement:

• bc_stencil(bc, I, ax, dim): return an iterable of (weight, index) pairs whose weighted sum Σ wⱼ · ϕ[Iⱼ] gives the ghost value at the out-of-bounds index I along dimension dim (with ax the in-bounds range). dim is part of the signature because an index can be out of bounds in several dimensions at once (e.g. a corner ghost), and each returned index is resolved along the remaining dimensions by the caller.

bc_stencil is a pure function of indices: it never touches the field, so the field-reading and accumulation live entirely in the indexing code (see _getindexbc).

struct CartesianCell{N, T}

A cell of a CartesianGrid: the axis-aligned hypercube bounded by nodes at lc (lower corner) and hc (upper corner). Obtain via getcell(grid, I) where I is a cell index.

CartesianGrid(lc, hc, n)

Create a uniform cartesian grid with lower corner lc, upper corner hc and n nodes in each direction.

Node coordinates are accessed via getnode; cell geometry via getcell. Use nodeindices to iterate over all node indices and cellindices to iterate over all cell indices.

Examples

using LevelSetMethods
a = (0, 0)
b = (1, 1)
n = (10, 4)
grid = CartesianGrid(a, b, n)

# output

CartesianGrid in ℝ²
  ├─ domain:  [0.0, 1.0] × [0.0, 1.0]
  ├─ nodes:   10 × 4
  └─ spacing: h = (0.1111, 0.3333)

CartesianGrid(lc, hc; meshsize)

Create a uniform cartesian grid spanning [lc, hc] with node spacing approximately equal to meshsize in each dimension. meshsize may be a scalar (applied to every dimension) or a collection with one entry per dimension.

The domain [lc, hc] is honored exactly; the number of cells in each dimension is rounded up, so the realized spacing (see meshsize) is never coarser than meshsize but may be finer.

Examples

using LevelSetMethods
CartesianGrid((0, 0), (1, 1); meshsize = 0.3)

# output

CartesianGrid in ℝ²
  ├─ domain:  [0.0, 1.0] × [0.0, 1.0]
  ├─ nodes:   5 × 5
  └─ spacing: h = (0.25, 0.25)

struct CurvatureTerm{V} <: LevelSetTerm

Level-set curvature term representing bκ|∇ϕ|, where κ = ∇ ⋅ (∇ϕ/|∇ϕ|) is the curvature.

struct EikonalReinitializationTerm <: LevelSetTerm

A LevelSetTerm representing sign(ϕ)(|∇ϕ| - 1), which drives the level set toward a signed distance function by solving the Eikonal equation |∇ϕ| = 1 via pseudo-time marching.

There are two constructors:

• EikonalReinitializationTerm(ϕ₀): freezes the sign from the initial level set ϕ₀ (equation 7.5 of Osher & Fedkiw). Recommended when the interface may drift. ϕ₀ may be a full-grid or narrow-band field; the frozen sign lives on the same active set.
• EikonalReinitializationTerm(): recomputes the sign from the current ϕ at each step (equation 7.6 of Osher & Fedkiw).

struct ExtrapolationBC{P} <: BoundaryCondition

Degree-P one-sided polynomial extrapolation boundary condition. Uses the P+1 nearest interior cells to construct a degree-P polynomial that is extrapolated into the ghost region.

ExtrapolationBC{0} (aliased as NeumannBC) gives constant extension (∂ϕ/∂n = 0 at the boundary face). ExtrapolationBC{1} (aliased as LinearExtrapolationBC) gives linear extrapolation (∂²ϕ/∂n² = 0).

struct ForwardEuler

First-order explicit Forward Euler time integration scheme.

using LevelSetMethods
ForwardEuler()

# output

ForwardEuler (1st order explicit)
  └─ cfl: 0.5

mutable struct InterpolatedField{F, T, B}

A continuous field obtained by equipping a discrete AbstractMeshField with piecewise polynomial interpolation. Evaluating cf(x) returns the value of the local BernsteinPolynomial on the cell containing x; gradient and hessian differentiate the same local patch.

The Bernstein basis is chosen for its convex-hull property: the interpolant on a cell is bounded by the extrema of its coefficients, which proven_empty exploits to discard cells with no interface or interior.

Construct via InterpolatedField(ϕ, order).

mutable struct InterpolationScratch{N, T}

Per-task working memory and per-cell memo for evaluating an InterpolatedField. It holds the scratch arrays written on every cell fill (coeffs, vals, temp1, temp2) together with the memo key (Ic, gen) identifying which cell the current coeffs are for. It carries no interpolation matrix (that is shared by the field) and no reference to the underlying field.

Construct via InterpolationScratch(N, order, T).

InterpolationScratch(N::Int, order::Int, T::Type)

Allocate the per-task scratch for N-dimensional fields with element type T and polynomial interpolation of degree order.

LevelSetEquation(; terms, ic, bc, t = 0, integrator = RK2())

Create a level-set equation of the form ϕₜ + sum(terms) = 0, where each t ∈ terms is a LevelSetTerm and ic is the initial condition — either a MeshField for a full-grid discretization or a NarrowBandMeshField for a narrow-band discretization.

Calling integrate!(eq, tf) will evolve the equation up to time tf, modifying current_state(eq) and current_time(eq) in place. The initial condition ic is copied into the equation, so the field you pass in is never modified; reach the evolving state through current_state.

Boundary conditions come from the bc keyword, or from ic if it already carries them (passing both warns and bc wins); it is an error if neither supplies them. If a single BoundaryCondition is provided, it will be applied uniformly to all boundaries of the domain. To apply different boundary conditions to each boundary, pass a tuple of the form (bc_x, bc_y, ...) with as many elements as dimensions in the domain. If bc_x is a BoundaryCondition, it will be applied to both boundaries in the x direction. If bc_x is a tuple of two BoundaryConditions, the first will be applied to the left boundary and the second to the right boundary. The same logic applies to the other dimensions.

The optional parameter t specifies the initial time of the simulation, and integrator is the TimeIntegrator used to evolve the level-set equation.

using LevelSetMethods, StaticArrays
grid = CartesianGrid((-1, -1), (1, 1), (50, 50))    # define the grid
ϕ = MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid)    # initial shape
𝐮 = MeshField(x -> SVector(1, 0), grid)             # advection velocity
terms = (AdvectionTerm(𝐮),)                          # advection term
bc = NeumannBC()                                     # zero-gradient boundary conditions
eq = LevelSetEquation(; terms, ic = ϕ, bc)          # level-set equation

# output

LevelSetEquation
  ├─ equation: ϕₜ + 𝐮 ⋅ ∇ ϕ = 0
  ├─ time:     0.0
  ├─ integrator: RK2 (2nd order TVD Runge-Kutta, Heun's method)
  │  └─ cfl: 0.5
  ├─ state: MeshField on CartesianGrid in ℝ²
  │  ├─ domain:  [-1.0, 1.0] × [-1.0, 1.0]
  │  ├─ nodes:   50 × 50
  │  ├─ spacing: h = (0.04082, 0.04082)
  │  ├─ bc:     Neumann (all)
  │  ├─ valtype: Float64
  │  └─ values:  min = -0.2492,  max = 1.75
  ╰─

abstract type LevelSetTerm

A typical term in a level-set evolution equation.

LinearExtrapolationBC = ExtrapolationBC{1}

Alias for ExtrapolationBC{1}: linear extrapolation into ghost cells. Corresponds to ∂²ϕ/∂n² = 0 at the boundary face.

struct MeshField{N,T,V,B} <: AbstractMeshField{N,T,V}

A node-centered field defined on the entire mesh, described by its discrete values at each node.

• vals: dense Array{V,N} of values at each node.
• mesh: the underlying CartesianGrid{N,T}.
• bcs: boundary conditions, used for indexing outside the mesh bounds.

Base.getindex(ϕ, I) returns the value of the field at node index I. If I lies outside the mesh bounds, bcs are applied to determine the value.

Use nodeindices and cellindices to iterate over active node and cell indices respectively.

MeshField(vals, grid::CartesianGrid; bc=nothing)

Construct a node-centered MeshField with explicit values on a grid.

• bc: boundary conditions (normalized via _normalize_bc).

MeshField(f::Function, grid; kwargs...)

Create a node-centered MeshField by evaluating f at each node of grid. All keyword arguments are passed to the vals-based constructor.

Examples

using LevelSetMethods, StaticArrays
grid = CartesianGrid((-1, -1), (1, 1), (5, 5))
# scalar field without boundary conditions
MeshField(x -> x[1]^2 + x[2]^2 - 0.5^2, grid)

# output

MeshField on CartesianGrid in ℝ²
  ├─ domain:  [-1.0, 1.0] × [-1.0, 1.0]
  ├─ nodes:   5 × 5
  ├─ spacing: h = (0.5, 0.5)
  ├─ valtype: Float64
  └─ values:  min = -0.25,  max = 1.75

struct NarrowBandMeshField{N,T,V,B} <: AbstractMeshField{N,T,V}

A node-centered field defined on a narrow band: a sparse set of active nodes, with values stored only there. The fields are:

• vals: sparse Dict{CartesianIndex{N},V} mapping active node indices to values.
• mesh: the underlying CartesianGrid{N,T}.
• bcs: boundary conditions, used for indexing outside the mesh bounds.
• nlayers: halo depth of the band — the number of node layers grown around the cut cells by update_band!, and the depth it rebuilds to during evolution.

Use active_nodeindices and active_cellindices to iterate over active node and cell indices respectively.

Indexing is defined on the band and the surrounding stencil halo; an in-grid node farther out throws. At present, periodic boundary conditions are not supported.

NarrowBandMeshField(vals, grid; bc=nothing, nlayers=3)

Construct a NarrowBandMeshField from a pre-built vals dict whose keys are the active nodes. Useful for an auxiliary field on an adopted band (e.g. a velocity), where the active set is filled from another field rather than derived from vals.

NarrowBandMeshField(f::Function, grid::CartesianGrid; bc=nothing, nlayers=3)

Construct a narrow-band field by evaluating f at each node of grid, then restricting to a band of nlayers around the interface. Equivalent to NarrowBandMeshField(MeshField(f, grid; bc); nlayers).

NarrowBandMeshField(ϕ::MeshField; nlayers = 3)

Construct a topological narrow-band field from a full-grid MeshField. The narrow-band field is defined on the band and a halo of nlayers nodes around the cut-cells.

const NeumannBC = ExtrapolationBC{0}

Homogeneous Neumann boundary condition (∂ϕ/∂n = 0 at the boundary face). Alias for ExtrapolationBC{0}: ghost cells take the value of the nearest boundary node (constant extension).

struct NewtonSDF{Φ,Tr,P,T}

A signed distance function to the zero level set of an underlying level set function, computed using a Newton-based closest point method.

Evaluating sdf(x) returns the signed distance from point x to the interface. An optional second argument sdf(x, s) can be used to supply the sign directly (e.g. sign(ϕ(x))) when it is already known, avoiding an extra interpolant evaluation.

NewtonSDF(ϕ::AbstractMeshField; order=3, upsample=2, maxiters=10, xtol=nothing, ftol=nothing)

Construct a NewtonSDF from an [AbstractMeshField]. The interface is sampled by projecting uniformly-spaced points in each cell onto the zero level set, building a KDTree for fast nearest-neighbour queries. Given a query point, the closest interface point is found by seeding a Newton-Lagrange solve from the nearest sample.

Keyword arguments

• order: polynomial interpolation order
• upsample: sampling density per cell side
• maxiters: maximum Newton iterations
• xtol: tolerance on the KKT residual; defaults to sqrt(eps(T)) with T = valtype(ϕ)
• ftol: tolerance on the function value; defaults to sqrt(eps(T)) with T = valtype(ϕ)

struct NormalMotionTerm{V,F} <: LevelSetTerm

Level-set advection term representing v |∇ϕ|. This LevelSetTerm should be used for internally generated velocity fields; for externally generated velocities you may use AdvectionTerm instead.

If passed, update_func will be called as update_func(v, ϕ, t) before computing the term at each stage of the time evolution.

struct PeriodicBC <: BoundaryCondition

Singleton type representing periodic boundary conditions. Ghost cells are filled by wrapping around to the opposite side of the domain.

struct RK2

Second order total variation dimishing Runge-Kutta scheme, also known as Heun's predictor-corrector method.

using LevelSetMethods
RK2()

# output

RK2 (2nd order TVD Runge-Kutta, Heun's method)
  └─ cfl: 0.5

struct RK3

Third order total variation dimishing Runge-Kutta scheme.

using LevelSetMethods
RK3()

# output

RK3 (3rd order TVD Runge-Kutta)
  └─ cfl: 0.5

struct SemiImplicitI2OE

Semi-implicit finite-volume scheme of the I2OE family ([3]) for advection problems.

using LevelSetMethods
SemiImplicitI2OE()

# output

SemiImplicitI2OE (semi-implicit advection, Mikula et al.)
  └─ cfl: 2.0

struct SymmetryBC <: BoundaryCondition

Symmetry-plane (reflective) boundary condition: the field is mirrored across the boundary node, ϕ[b - k] = ϕ[b + k]. Like NeumannBC it enforces ∂ϕ/∂n = 0, but by reflection rather than flat extension, so it also preserves curvature symmetrically and makes the interface meet the boundary perpendicularly. This is the correct condition on a symmetry axis, e.g. for axisymmetric simulations.

abstract type TimeIntegrator end

Abstract type for time integrators. See subtypes(TimeIntegrator) for a list of available time integrators. Every concrete integrator stores a cfl safety factor.

Upwind <: SpatialScheme

First-order upwind discretization of the spatial derivatives in a level-set term. Cheap and robust, but introduces significant numerical diffusion; prefer WENO5 unless cost is the overriding concern.

WENO5 <: SpatialScheme

Fifth-order WENO (Weighted Essentially Non-Oscillatory) discretization of the spatial derivatives in a level-set term. More expensive than Upwind, but much more accurate in smooth regions while suppressing oscillations near steep gradients.

Base.copy!(dest::AbstractMeshField, src::AbstractMeshField)

Make dest hold the same values as src, leaving its mesh and boundary conditions untouched. For a sparse NarrowBandMeshField this also replaces dest's active set with src's (the Dict copy! syncs keys), so the two agree on which nodes are stored — the time integrator relies on this to keep a reused stage buffer's band in sync with the evolving field.

Base.copy(ϕ::AbstractMeshField)

Clone ϕ. Generic fallback (deepcopy); concrete fields override it to share their immutable mesh and boundary conditions and only copy the values.

copy(ϕ::MeshField)

Return a MeshField with an independent copy of the values of ϕ, sharing its (immutable) mesh and boundary conditions.

copy(nb::NarrowBandMeshField)

Return a NarrowBandMeshField with an independent copy of the band values of nb, sharing its (immutable) mesh and boundary conditions.

intersect!(ϕ₁::MeshField, ϕ₂::MeshField)

In-place intersection of two level sets: ϕ₁ = max(ϕ₁, ϕ₂).

intersect(ϕ₁::MeshField, ϕ₂::MeshField)

Intersection of two level sets: max(ϕ₁, ϕ₂).

map(f, ϕ::MeshField)

Return a MeshField on the same mesh and boundary conditions as ϕ, with each value replaced by f applied to it. The element type may change with f.

map(f, nb::NarrowBandMeshField)

Return a NarrowBandMeshField with the same mesh, boundary conditions and band as nb, with each stored value replaced by f applied to it. The element type may change with f.

setdiff!(ϕ₁::MeshField, ϕ₂::MeshField)

In-place set difference: ϕ₁ = max(ϕ₁, -ϕ₂).

setdiff(ϕ₁::MeshField, ϕ₂::MeshField)

Set difference of two level sets: max(ϕ₁, -ϕ₂).

union!(ϕ₁::MeshField, ϕ₂::MeshField)

In-place union of two level sets: ϕ₁ = min(ϕ₁, ϕ₂).

union(ϕ₁::MeshField, ϕ₂::MeshField)

Union of two level sets: min(ϕ₁, ϕ₂).

D2(ϕ,I,dims)

Finite difference scheme for second order derivative at grid point I along the dimensions dims.

If dims[1] == dims[2], it is more efficient to call D2⁰(ϕ,I,dims[1]).

D2⁰(ϕ,I,dim)

Centered finite difference scheme for second order derivative at grid point I along dimension dim. E.g. if dim=1, this approximates ∂ₓₓ.

D2⁺⁺(ϕ,I,dim)

Forward finite difference scheme for second order derivative at grid point I along dimension dim. E.g. if dim=1, this approximates ∂ₓₓ.

D2⁻⁻(ϕ,I,dim)

Backward finite difference scheme for second order derivative at grid point I along dimension dim. E.g. if dim=1, this approximates ∂ₓₓ.

D⁰(ϕ,I,dim)

Centered finite difference scheme for first order derivative at grid point I along dimension dim.

D⁺(ϕ,I,dim)

Forward finite difference scheme for first order derivative at grid point I along dimension dim.

D⁻(ϕ,I,dim)

Backward finite difference scheme for first order derivative at grid point I along dimension dim.

_add_boundary_conditions(ϕ::MeshField, bc)

Return a new MeshField with bc as boundary conditions. All of the underlying data is aliased (shared) with the original MeshField.

_add_boundary_conditions(nb::NarrowBandMeshField, bc)

Return a new NarrowBandMeshField with bc as boundary conditions. All underlying data is aliased with the original.

_bc_str(bcs)

Format a boundary-conditions tuple (as returned by boundary_conditions) into a compact human-readable string. Each element of bcs is a (left, right) pair for one spatial dimension.

_closest_point(p::BernsteinPolynomial, xq, x0, maxiters, xtol, ftol, safeguard_dist) -> (x_closest, converged)

Find the point on the zero level-set of the single-cell patch p closest to xq, starting from x0. Uses a Newton-Lagrange solver on the KKT conditions of min ||x - xq||² s.t. p(x) = 0. The patch is fixed: unlike _project_to_interface, the iterate is not re-dispatched to neighbouring cells.

_closest_point_from_seed(sdf, x, seed)

Run the Newton-Lagrange closest-point solver on the Bernstein patch of the cell containing seed (a point already on or near the interface). Returns (closest_point, converged, ∇φ(closest_point)), where the gradient gives the (unnormalized, outward) interface normal at the closest point.

_closest_point_on_interface(sdf, x)

Find the point on the interface closest to x by seeding a local Newton-Lagrange solve from the nearest interface sample. Returns (closest_point, converged, ∇φ(closest_point)).

_domain_str(g::CartesianGrid) -> String

Format the domain as "[lo, hi] × [lo, hi] × …".

_embed_show(io, label, obj; indent="  ")

Print the text/plain representation of obj indented under label as a tree branch. The first line becomes ├─ label: <header>, and any remaining lines are indented with │.

_extrapolate_to_ghost(nb::NarrowBandMeshField, I) -> value

Extrapolate a value at the out-of-band index I so that finite-difference/WENO stencils reaching just past the band edge get a consistent value. It is the linear (affine) Lagrange extrapolant anchored at the nearest band node.

_getindexbc(ϕ::AbstractMeshField, I, Val(dim)) -> value

Recursive resolutions ϕ at an out-of-grid index I via boundary conditions. The recursion peels one dimension at a time from dim down to 0; each level that is in range is skipped, and each out-of-range level is resolved by accumulating the boundary condition's bc_stencil, recursing on the remaining dimensions so corner ghosts (out of bounds in several dimensions) compose correctly. At level 0 every component is back in range, so the value is read through the field's own in-grid getindex (ϕ[I]).

dim is carried as a Val so that it is a compile-time constant: the per-dimension lookups axes(ϕ)[dim] and boundary_conditions(ϕ)[dim] then resolve to concrete types, keeping the boundary-condition dispatch static even when the conditions differ between dimensions.

_i2oe_relation(bc, Ioff, I, ax, dim, T) -> (α, β, idx, γ)

Express the ghost value of the SemiImplicitI2OE scheme at the out-of-bounds neighbor Ioff of node I (along dimension dim, with in-bounds range ax and element type T) as the affine relation α·u[I] + β·u[idx] + γ. Optional interface method of BoundaryCondition; the fallback rejects conditions the integrator cannot handle.

_ingrid_value(nb::NarrowBandMeshField, I) -> value

Value at the in-grid index I. Returns the stored value if I is a band node; otherwise extrapolates it from the nearest band node (see _extrapolate_to_ghost).

_interpolation_matrix(order, T)

Build the 1D interpolation matrix mapping the stencil_order + 1 equispaced nodal values on a super-cell to the order + 1 coefficients of the Bernstein basis on the central cell.

The matrix depends only on order and the element type T — not on the cell or the spatial dimension — so it is built once per InterpolatedField and shared (read-only) by every evaluating task.

_itp_buffer(cf::InterpolatedField)

Return the calling task's private InterpolationScratch for cf, lazily created on first use and cached per task by the OncePerTask buffer. No lock is needed and concurrent tasks never share scratch memory.

_lagrange_extrap_weight(j, k, P)

Lagrange weight for the j-th interior node (0-indexed) when extrapolating to a ghost point at distance k outside the boundary, using a degree-P polynomial fitted to P+1 nodes.

Nodes are at positions 0, 1, …, P (relative to the boundary); the ghost is at position -k. The weight is

TODO: maybe write this as proper latex math?

wⱼ = ∏_{m=0, m≠j}^{P}  (-k - m) / (j - m)

_normalize_bc(bc, dim)

Normalize the bc argument into a dim-tuple of (left, right) pairs, one per spatial dimension, as expected by _add_boundary_conditions.

• A single BoundaryCondition is applied to all sides.
• A length-dim collection applies each entry to both sides of the corresponding dimension.
• A length-dim collection of 2-tuples applies each entry as (left, right) for that dimension.

_normalize_bc_pair(bc, dim) -> (left, right)

Normalize the boundary condition(s) for a single dimension dim into a (left, right) pair. A single BoundaryCondition applies to both sides; a length-2 collection is taken as (left, right). Periodicity must hold on both sides of a dimension or neither.

_project_to_interface(itp::InterpolatedField, x_start, maxiters, ftol, safeguard_dist)

Use Newton's method to project a starting point onto the zero level set of itp. Returns the converged point or nothing if Newton fails to converge or if the iterate moves more than safeguard_dist from x_start. The full InterpolatedField (not a single-cell patch) is used so the projection follows the iterate across cells.

_sample_cell!(pts, lk, grid, itp, I, upsample, maxiters, ftol, safeguard_dist)

Project uniformly-spaced sample points in cell I onto the interface, appending to the shared pts (under lock lk) the converged projections that land back inside I. Cells proven to contain no interface are skipped.

_sample_interface(grid, itp, cells, upsample, maxiters, ftol)

Project uniformly-spaced sample points in each candidate cell onto the interface. Returns all converged projections; cells proven to be empty are skipped. Cells are processed in parallel, appending to a single shared vector guarded by a lock.

_show_fields(io, g::CartesianGrid; prefix="  ", last=true)

Print domain, nodes, and spacing of g as indented tree lines. When last=true the spacing line uses └─ (terminal); otherwise ├─ (continuing).

_superscript(n::Int) -> String

Convert an integer to its Unicode superscript representation, e.g. 2 → "²".

active_cellindices(ϕ::AbstractMeshField)

Return the cell indices ϕ is defined on. Defaults to cellindices (every grid cell); sparse fields override it with their active set.

active_cellindices(nb::NarrowBandMeshField)

Return the set of active cell indices of nb: cells whose corners are all in-band nodes.

active_nodeindices(ϕ::AbstractMeshField)

Return the node indices ϕ is defined on. Defaults to nodeindices (every grid node); sparse fields such as NarrowBandMeshField override it with their active set.

active_nodeindices(nb::NarrowBandMeshField)

Return the set of active (in-band) node indices of nb.

bc_stencil(bc::BoundaryCondition, I, ax, dim) -> NTuple of (weight, index)

Express the ghost value at the out-of-bounds index I (along dimension dim, with in-bounds range ax) as a weighted sum Σ wⱼ · ϕ[Iⱼ], returned as a tuple of (weight, index) pairs. Each returned index is in range along dim but may still be out of bounds in other dimensions; the caller (_getindexbc) resolves those by recursion. Primary interface method of BoundaryCondition; pure function of indices, with no field access.

boundary_conditions(eq::LevelSetEquation)

Return the boundary conditions of the equation.

cell_extrema(cf::InterpolatedField, I::CartesianIndex)

Compute the minimum and maximum values of the local Bernstein interpolant in cell I.

cellindices(ϕ::AbstractMeshField)

Return all cell indices of the underlying mesh. See active_cellindices for the subset a field is actually defined on.

cellindices(g::CartesianGrid)

Return a CartesianIndices ranging over all cell indices of g. Cell I is the hypercube bounded by nodes I and I + 1 in each dimension. Cells are indexed 1:n[d]-1 in each dimension d.

complement!(ϕ::MeshField)

In-place complement of a level set: ϕ = -ϕ.

complement(ϕ::MeshField)

Complement of a level set: -ϕ.

compute_cfl(terms, ϕ::MeshField, t)

Compute the maximum stable time-step $Δt$ for the given terms and level set ϕ at time t, based on the Courant-Friedrichs-Lewy (CFL) condition.

compute_index(ϕ::AbstractMeshField, x)

Return the multi-index of the cell containing the point x, delegating to compute_index(g::CartesianGrid, x) on the underlying mesh.

compute_index(g::CartesianGrid, x)

Return the multi-index of the cell of g containing the point x. If x lies outside the grid, the index of the closest cell is returned (clamped to the boundary).

current_state(eq::LevelSetEquation)

Return the current state of the level-set equation (a MeshField).

current_time(eq::LevelSetEquation)

Return the current time of the simulation.

curvature(ϕ::AbstractMeshField, I::CartesianIndex)

Mean curvature κ = ∇ ⋅ (∇ϕ / ‖∇ϕ‖) of the level set of ϕ at grid index I, computed from the gradient and hessian (centered finite differences) via

κ = (Δϕ ‖∇ϕ‖² - ∇ϕᵀ Hϕ ∇ϕ) / ‖∇ϕ‖³ .

Returns zero where ∇ϕ vanishes, since the curvature is undefined there. See the Wikipedia article on the implicit form of the mean curvature.

export_surface_mesh(eq::LevelSetEquation, filename; kwargs...)

Export a surface mesh of the 3D interface (where ϕ = 0) to filename.

export_volume_mesh(eq::LevelSetEquation, filename; kwargs...)

Export a 3D volume mesh of the interior domain (where ϕ < 0) to filename.

extend_along_normals!(F, ϕ::MeshField;
                      nb_iters = 50,
                      cfl = 0.45,
                      frozen = nothing,
                      interface_band = 1.5,
                      min_norm = 1e-14)

Extend a scalar speed field F away from the interface of ϕ by solving in pseudo-time ∂τF + sign(ϕ) n⋅∇F = 0, with n = ∇ϕ / |∇ϕ|.

The equation is discretized with first-order upwind derivatives. The update preserves frozen nodes (Dirichlet constraint). If frozen is not provided, a mask is built from the interface band abs(ϕ) <= interface_band * Δ, where Δ = minimum(meshsize(ϕ)).

F can be an AbstractArray (same size as ϕ) or a MeshField on the same mesh.

Reference from [Peng et al. 1999]

fill_coefficients!(ϕ::AbstractMeshField, mat, scratch::InterpolationScratch, base_idxs::CartesianIndex)

Fill scratch.coeffs with the Bernstein coefficients for the cell at base_idxs, using the shared interpolation matrix mat and the surrounding stencil of grid values of ϕ.

get_sample_points(sdf::NewtonSDF)

Return the interface sample points, lying on the surface, used to build the KDTree of sdf.

getcell(ϕ::AbstractMeshField, I)

Return the cell at index I, delegating to getcell on the underlying mesh.

getcell(g::CartesianGrid, I::CartesianIndex)

Return the CartesianCell with lower corner at node I and upper corner at node I+1. I must be a valid cell index, i.e. I ∈ cellindices(g).

getnode(ϕ::AbstractMeshField, I)

Return the coordinates of node I, delegating to getnode on the underlying mesh.

getnode(g::CartesianGrid, I::CartesianIndex)
getnode(g::CartesianGrid, i1, i2, ...)

Return the coordinates of the node at multi-index I as an SVector. I must be a valid node index (i.e. I ∈ nodeindices(g)).

gradient(cf::InterpolatedField, x)

Evaluate the spatial gradient vector of the interpolated field at the point x.

gradient(p::BernsteinPolynomial, x)

Compute the gradient of p at point x using ForwardDiff.

gradient(ϕ::AbstractMeshField, I::CartesianIndex)

Gradient ∇ϕ at grid index I from centered finite differences, returned as an SVector.

grid1d(g::CartesianGrid[, dim::Integer])

Return a LinRange of the coordinates along the given dimension dim. If dim is not provided, return a tuple of LinRanges for all dimensions.

hausdorff_distance(sdf₁::NewtonSDF, sdf₂::NewtonSDF) -> Float64

Approximate the Hausdorff distance between the zero level sets Γ₁ and Γ₂ of sdf₁ and sdf₂. The symmetric Hausdorff distance is

d_H(Γ₁, Γ₂) = max(max_{p ∈ Γ₁} d(p, Γ₂),  max_{p ∈ Γ₂} d(p, Γ₁))

Each one-sided maximum is approximated by iterating over the interface sample points of one NewtonSDF and finding the closest point on the other interface via the Newton-Lagrange solver.

hessian(cf::InterpolatedField, x)

Evaluate the spatial Hessian matrix of the interpolated field at the point x.

hessian(p::BernsteinPolynomial, x)

Compute the hessian of p at point x using ForwardDiff.

hessian(ϕ::AbstractMeshField, I::CartesianIndex)

Hessian Hϕ = ∇∇ϕ at grid index I from second-order centered finite differences, returned as a Symmetric SMatrix.

integrate!(ls::LevelSetEquation, tf, Δt = Inf; prehook = identity, posthook = identity)

Integrate the LevelSetEquation ls up to time tf, mutating the levelset and current_time of the object ls in the process.

An optional parameter Δt can be passed to specify a maximum time-step allowed for the integration. Note that the internal time-steps taken to evolve the level-set up to tf may be smaller than Δt due to stability reasons related to the terms and integrator employed.

prehook and posthook are functions called once per accepted time step (as opposed to once per Runge-Kutta stage): prehook(ls) runs at the start of the step, before the state is advanced, and posthook(ls) runs after the step has been committed. Both receive ls with current_state(ls) and current_time(ls) synced to that point, and may mutate current_state(ls) (mutations are carried into the step). Their return values are ignored; the default identity is a no-op.

Reinitialization of the level-set to a signed distance function is not built in by default, but can be easily added by creating a prehook: pass prehook = eq -> reinitialize!(current_state(eq)) (see reinitialize!) to reinitialize before every step, or gate the call on any criterion (elapsed time, |∇ϕ| drift, a step counter closed over by the hook). The integrator never reinitializes on its own — keeping a signed distance function is entirely the caller's choice.

isrealvalued(ϕ::AbstractMeshField)

Return true if the values of ϕ are real numbers. Real-valued fields are the ones that may be interpreted as a level set (with their zero contour representing an interface).

local_interpolant(cf::InterpolatedField, x)

Return the local BernsteinPolynomial representing the field in the cell containing the point x.

make_interpolant(cf::InterpolatedField, I::CartesianIndex)

Return a BernsteinPolynomial for the cell at multi-index I, lazily computing (and caching) its Bernstein coefficients from the surrounding stencil of grid values.

makie_theme()

Return a Makie theme for plots of level-set functions.

mesh(eq::LevelSetEquation)

Return the underlying CartesianGrid of the equation.

meshsize(g::CartesianGrid[, dim::Integer])

Return the spacing between grid nodes along the given dimension dim. If dim is not provided, return a SVector of spacings for all dimensions.

nlayers(nb::NarrowBandMeshField)

Return the halo depth of the band (see NarrowBandMeshField).

nodeindices(g::CartesianGrid)

Return a CartesianIndices ranging over all node indices of g. Nodes are indexed 1:n[d] in each dimension d.

nodeindices(ϕ::AbstractMeshField)

Return all node indices of the underlying mesh. See active_nodeindices for the subset a field is actually defined on.

normal(ϕ::AbstractMeshField, I::CartesianIndex)

Unit exterior normal n = ∇ϕ / ‖∇ϕ‖ of the level set of ϕ at grid index I.

perimeter(ϕ::MeshField)

Measure of the interface {x : ϕ(x) = 0} described by the zero level set of ϕ. In N dimensions this is the (N-1)-dimensional measure — a perimeter (arc length) in 2D, a surface area in 3D — approximated as ∫ δ(ϕ) ‖∇ϕ‖ dx with a smoothed Dirac delta δ. Contributions from the domain border are neglected.

using LevelSetMethods
R = 0.5
S₀ = 2π * R
grid = CartesianGrid((-1, -1), (1, 1), (200, 200))
ϕ = MeshField(x -> sqrt(sum(abs2, x)) - R, grid)
LevelSetMethods.perimeter(ϕ), S₀

# output

(3.1426415491430384, 3.141592653589793)

perimeter(nb::NarrowBandMeshField)

Measure of the interface {x : ϕ(x) = 0} described by the zero level set, matching perimeter(ϕ::MeshField). The smoothed Dirac delta is supported only within δmin of the interface — entirely inside the band — so the sum runs over the active nodes alone.

proven_empty(cf::InterpolatedField, I::CartesianIndex; surface=false)

Return true if cell I is guaranteed to contain no interface (when surface=true) or no interior (when surface=false), based on the convex-hull property of the Bernstein basis.

quadrature(mf::AbstractMeshField; interpolation_order, quadrature_order, surface=false)
quadrature(ϕ::InterpolatedField; quadrature_order, surface=false)

Generate a quadrature for the implicit domain defined by the level set. If surface=true, generate a quadrature for the interface ϕ=0; otherwise for the interior ϕ < 0.

The first form wraps mf in an InterpolatedField of degree interpolation_order and quadratures over it; pass an InterpolatedField directly (second form) to reuse an existing interpolant and skip interpolation_order. For a high-order rate, keep quadrature_order ≥ interpolation_order so the rule does not cap the interpolant's accuracy.

Returns a Dict mapping each cut cell's CartesianIndex to its ImplicitIntegration.Quadrature (which exposes coords and weights); provably empty cells are omitted.

reinitialize!(ϕ::AbstractMeshField; order=3, upsample=2, maxiters=20, xtol=nothing, ftol=nothing)

Reinitialize the level set ϕ to a signed distance function in place using the Newton closest-point method: a NewtonSDF is built from ϕ and each active grid node is overwritten with its signed distance to the interface (a Newton-Lagrange solve seeded from the nearest interface sample). The keyword arguments are forwarded to NewtonSDF.

Reinitialization is not built into the time integrator; drive it from a prehook, e.g. prehook = eq -> reinitialize!(current_state(eq); order = 5). See integrate!.

set_makie_theme!()

Set the Makie theme to LevelSetMethods.makie_theme().

terms(eq::LevelSetEquation)

Return the tuple of LevelSetTerms of the equation.

time_integrator(eq::LevelSetEquation)

Return the TimeIntegrator of the equation.

update_band!(ϕ::AbstractMeshField; nlayers)

Rebuild the active node set topologically: every node within nlayers axis-aligned steps of a cut cell (one whose corner values straddle zero). ϕ must be real-valued (a level set). Values are filled via ϕ's own indexing (stored value or affine extrapolation).

nlayers defaults to the band's stored halo depth (see NarrowBandMeshField).

A full-grid MeshField already holds every node, so there is no band to rebuild and this is a no-op — making it safe to call regardless of the state type. integrate! calls it automatically after each step.

update_term!(term::LevelSetTerm, ϕ, t)

Update the internal state of a LevelSetTerm before computing its contribution. Called before the CFL estimate and at each stage of the time integration.

value_and_gradient(cf::InterpolatedField, x)

Evaluate the value and spatial gradient vector of the interpolated field at the point x.

value_and_gradient(p::BernsteinPolynomial, x)

Fused computation of the value and gradient of p at point x using ForwardDiff.

value_gradient_hessian(cf::InterpolatedField, x)

Evaluate the value, spatial gradient vector, and Hessian matrix of the interpolated field at the point x.

value_gradient_hessian(p::BernsteinPolynomial, x)

Fused computation of the value, gradient, and hessian of p at point x using ForwardDiff.

volume(ϕ::MeshField)

Measure of the region {x : ϕ(x) ≤ 0} enclosed by the zero level set of ϕ. In N dimensions this is the N-dimensional measure — a length in 1D, an area in 2D, a volume in 3D — approximated as ∫ H(-ϕ) dx with a smoothed Heaviside H.

using LevelSetMethods
R = 0.5
V₀ = π * R^2
grid = CartesianGrid((-1, -1), (1, 1), (200, 200))
ϕ = MeshField(x -> sqrt(sum(abs2, x)) - R, grid)
LevelSetMethods.volume(ϕ), V₀

# output

(0.7854362890190668, 0.7853981633974483)

volume(nb::NarrowBandMeshField)

Measure of the region {x : ϕ(x) ≤ 0} enclosed by the zero level set, matching volume(ϕ::MeshField) but computed from the band alone, without materialising values on the full grid.

The smoothed Heaviside transition is supported within δmin of the interface — entirely inside the band — so every off-band node contributes exactly 0 (outside) or 1 (inside), and only its sign is needed. Those signs are recovered by a scanline sweep: the interface never crosses outside the band, so along each grid line the off-band interior shows up as same-sign gaps between consecutive band nodes (and the tails beyond the outermost ones), counted by arithmetic without ever visiting an interior node. A line with no band node carries no crossing and is classified by the sign of its nearest band node.

weno5⁺(ϕ, I, dim)

Fifth-order WENO (Weighted Essentially Non-Oscillatory) reconstruction of the derivative at grid point I along dimension dim, using a right-biased stencil.

weno5⁻(ϕ, I, dim)

Fifth-order WENO (Weighted Essentially Non-Oscillatory) reconstruction of the derivative at grid point I along dimension dim, using a left-biased stencil.

export_surface_mesh(eq::LevelSetEquation, args...; kwargs...)

Call export_surface_mesh on current_state(eq).

export_surface_mesh(ϕ::MeshField, output::String;
    hgrad = nothing, hmin = nothing, hmax = nothing, hausd = nothing)

Compute a mesh of the MeshField ϕ zero contour using MMGs_O3.

hgrad control the growth ratio between two adjacent edges

hmin and hmax control the edge sizes to be (respectively) greater than the hmin parameter and lower than the hmax one

hausd control the maximal distance between the piecewise linear representation of the boundary and the reconstructed ideal boundary

export_volume_mesh(eq::LevelSetEquation, output; kwargs...)

Call export_volume_mesh on current_state(eq).

export_volume_mesh(ϕ::MeshField, output::String;
    hgrad = nothing, hmin = nothing, hmax = nothing, hausd = nothing)

Compute a mesh of the domains associated with MeshField eq using either MMG2dO3 or MMG3dO3.

hgrad control the growth ratio between two adjacent edges.

hmin and hmax control the edge sizes to be (respectively) greater than the hmin parameter and lower than the hmax one.

hausd control the maximal distance between the piecewise linear representation of the boundary and the reconstructed ideal boundary.

For more information, see the official MMG documentation.