Reference

abstract type AbstractDomain end

Abstract type for the domain of a MeshField.

abstract type AbstractMesh{N,T}

An abstract mesh structure in dimension N with primite data of type T.

AbstractSignedDistanceFunction <: Function

Abstract type for signed distance functions.

Subtypes should be callable with sdf(x), and must implement update! to rebuild from a new level set.

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 ϕ.

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

Types used to specify boundary conditions.

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.

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)

struct CurvatureTerm{V} <: LevelSetTerm

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

struct DirichletBC{T} <: BoundaryCondition

A Dirichlet boundary condition taking values of f(x, t) at the boundary, where x is the spatial coordinate and t is the current time.

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.
• 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).

using LevelSetMethods
ExtrapolationBC(4)

# output

Degree 4 extrapolation

struct ForwardEuler

First-order explicit Forward Euler time integration scheme.

using LevelSetMethods
ForwardEuler()

# output

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

struct FullDomain <: AbstractDomain

Represents a field defined on the entire mesh.

const LevelSet{N, T, B}

Alias for MeshField on a CartesianGrid{N,T} with values stored as an Array{T,N} and a FullDomain. B is the type of the boundary conditions.

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

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 LevelSet for a full-grid discretization or a NarrowBandLevelSet for a narrow-band one.

Calling integrate!(eq, tf) will evolve the equation up to time tf, modifying current_state(eq) and current_time(eq) in place.

Boundary conditions are specified via the required bc keyword. 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.

Reinitialization is controlled by the reinit keyword, which accepts:

• nothing (default): no automatic reinitialization.
• an Int: reinitialization every reinit steps using NewtonReinitializer with default settings.
• a NewtonReinitializer: full control over algorithm parameters and frequency.

using LevelSetMethods, StaticArrays
grid = CartesianGrid((-1, -1), (1, 1), (50, 50))    # define the grid
ϕ = LevelSet(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
  ├─ reinit:   none
  ├─ state: MeshField on CartesianGrid in ℝ²
  │  ├─ domain:  [-1.0, 1.0] × [-1.0, 1.0]
  │  ├─ nodes:   50 × 50
  │  ├─ spacing: h = (0.04082, 0.04082)
  │  ├─ bc:     Degree 0 extrapolation (all)
  │  ├─ eltype:  Float64
  │  └─ values:  min = -0.2492,  max = 1.75
  ├─ log: SimulationLog (empty)
  ╰─

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.

MeshField(f::Function, m)

Create a MeshField by evaluating a function f on a mesh m.

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)
  ├─ eltype:  Float64
  └─ values:  min = -0.25,  max = 1.75

using LevelSetMethods
grid = CartesianGrid((-1, -1), (1, 1), (5, 5))
# vector-valued field
MeshField(x -> (x[1], x[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)
  └─ eltype:  Tuple{Float64, Float64}

struct MeshField{V,M,B,D}

A field described by its discrete values on a mesh.

• vals: the discrete values of the field.
• mesh: the underlying mesh (e.g. CartesianGrid).
• bcs: boundary conditions, used for indexing outside the mesh bounds.
• domain: the domain on which the field is defined (e.g. FullDomain).

Base.getindex of an MeshField is overloaded to handle indices that lie outside the CartesianIndices of its MeshField by using bcs.

MeshField(vals, mesh, bcs)

Construct a MeshField with explicit values, mesh, and boundary conditions. Defaults to FullDomain.

struct NarrowBandDomain{T} <: AbstractDomain

Domain for a narrow-band level set.

• halfwidth: half-width of the narrow band, typically on the order of a few grid spacings.

Active indices are the keys of the associated values dict and need not be stored separately.

const NarrowBandLevelSet{N, T, B}

Alias for MeshField on a CartesianGrid{N,T} with values stored as a Dict{CartesianIndex{N},T} and a NarrowBandDomain{T}. B is the type of the boundary conditions.

NarrowBandLevelSet(f, grid::CartesianGrid, halfwidth::Real; bc = nothing)

Construct a NarrowBandLevelSet by evaluating f at each node of grid and keeping only those where |f(x)| < halfwidth. No dense array is allocated.

NarrowBandLevelSet(f, grid::CartesianGrid; nlayers = 8, bc = nothing)

Construct a NarrowBandLevelSet with halfwidth automatically computed as nlayers * minimum(meshsize(grid)).

NarrowBandLevelSet(ϕ::LevelSet, halfwidth::Real; reinitialize = true)

Construct a NarrowBandLevelSet from a full-grid LevelSet. Active nodes are those where |ϕ[I]| < halfwidth. Boundary conditions are inherited from ϕ.

If reinitialize is true (the default), ϕ is first reinitialized to a signed distance function using NewtonReinitializer.

NarrowBandLevelSet(ϕ::LevelSet; nlayers = 3, reinitialize = true)

Construct a NarrowBandLevelSet with halfwidth automatically computed as nlayers * minimum(meshsize(ϕ)). nlayers sets the number of cell layers on each side of the interface included in the band.

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 NewtonReinitializer

Reinitializes a level set to a signed distance function using a Newton closest-point method. At each reinitialization step, a piecewise polynomial interpolant is built from the current level set, the interface is sampled, and the level set values are overwritten with the signed distances.

Keyword arguments

• reinit_freq: reinitialization frequency in time steps (default: 1)
• order: polynomial interpolation order (default: 3)
• upsample: interface sampling density per cell side (default: 8)
• maxiters: maximum Newton iterations (default: 20)
• xtol: tolerance on the KKT residual (default: 1e-8)
• ftol: tolerance on the function value (default: 1e-8)

using LevelSetMethods
NewtonReinitializer()

# output

NewtonReinitializer
  ├─ frequency: every step
  ├─ order:     3
  ├─ upsample:  8×
  ├─ maxiters:  20
  ├─ xtol:      1.0e-8
  └─ ftol:      1.0e-8

mutable struct NewtonSDF{I,Tr,P} <: AbstractSignedDistanceFunction

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(ϕ; order=3, kwargs...)

Construct a NewtonSDF from a level set by first creating a piecewise polynomial interpolant of the given order. Additional keyword arguments are forwarded to NewtonSDF(itp; ...). Works for both LevelSet and NarrowBandLevelSet, sampling only the relevant candidate cells in each case.

NewtonSDF(itp; upsample=8, maxiters=20, xtol=1e-8, ftol=1e-8)

Construct a NewtonSDF from a PiecewisePolynomialInterpolant.

The interface is sampled by projecting uniformly-spaced points in each cell onto the zero level set. A KDTree is built from these samples for fast nearest-neighbor queries.

Keyword arguments

• upsample: sampling density per cell side. Larger values means a denser sampling of the interface is used to build the KDTree, which in turn usually means a better initial guess for the Newton solver.
• maxiters: maximum Newton iterations
• xtol: tolerance on iterate updates for convergence of the Newton solver
• ftol: tolerance on the function value for convergence of the Newton solver

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.

using LevelSetMethods
PeriodicBC()

# output

Periodic

struct PiecewisePolynomialInterpolant{Φ, N, T}

A piecewise polynomial interpolant built over a MeshField.

The domain is partitioned into cells of the underlying mesh. On each cell a BernsteinPolynomial is fitted to the surrounding stencil of grid values, so that evaluating the interpolant at a point x amounts to:

1. locating the cell containing x,
2. (lazily) computing the local Bernstein coefficients, and
3. evaluating the resulting polynomial.

Construct via interpolate. The returned object itp supports:

• itp(x) — evaluate at point x
• make_interpolant(itp, I) — local BernsteinPolynomial for cell I
• gradient(itp, x) / hessian(itp, x) — via ForwardDiff
• cell_extrema(itp, I) — certified bounds via convex-hull property
• proven_empty(itp, I) — certified absence of interface or interior

The struct is mutable because the local coefficients and stencil values are cached in place (coeffs, vals) to avoid allocation on repeated evaluations.

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 (Mikula et al.) for advection problems.

using LevelSetMethods
SemiImplicitI2OE()

# output

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

mutable struct SimulationLog

Accumulates per-step timing and progress data for a LevelSetEquation. Persists and accumulates across multiple calls to integrate!. Use reset_log! to clear the history.

struct StepRecord

Record of a single time step during level-set integration.

• step: cumulative step number across all integrate! calls (1-indexed)
• t: simulation time at end of step
• wall_time: total wall-clock time for this step (seconds)
• reinit_time: time spent in reinitialization (0.0 if not performed)
• did_reinit: whether reinitialization was performed this step
• update_times: time per update_term!, summed over all RK stages (seconds)
• compute_times: time per _compute_term loop, summed over all RK stages (seconds)
• ϕ_min, ϕ_max: level-set extrema at end of step

abstract type TimeIntegrator end

Abstract type for time integrators. See subtypes(TimeIntegrator) for a list of available time integrators.

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

Copy the values from src to dest. The meshes, boundary conditions, and domains of the dest fields are not modified.

intersect!(ϕ1::LevelSet, ϕ2::LevelSet)

In-place intersection of two level sets: $ϕ_1 = \max(ϕ_1, ϕ_2)$.

intersect(ϕ1::LevelSet, ϕ2::LevelSet)

Return the intersection of two level sets: $\max(ϕ_1, ϕ_2)$.

setdiff!(ϕ1::LevelSet, ϕ2::LevelSet)

In-place set difference: $ϕ_1 = \max(ϕ_1, -ϕ_2)$.

setdiff(ϕ1::LevelSet, ϕ2::LevelSet)

Return the set difference: $\max(ϕ_1, -ϕ_2)$.

union!(ϕ1::LevelSet, ϕ2::LevelSet)

In-place union of two level sets: $ϕ_1 = \min(ϕ_1, ϕ_2)$.

union(ϕ1::LevelSet, ϕ2::LevelSet)

Return the union of two level sets: $\min(ϕ_1, ϕ_2)$.

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)

Upward 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)

Second-order 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.

D⁻⁻(ϕ, I, dim)

Second-order backward finite difference scheme for first order derivative at grid point I along dimension dim.

_apply_extrapolation_bc(ϕ, I, ::ExtrapolationBC{P}, ax, dim)

Return the extrapolated value of ϕ at out-of-bounds index I in dimension dim using degree-P Lagrange extrapolation (see _lagrange_extrap_weight). The boundary node and its P interior neighbors are used as stencil nodes.

_base_lookup(nb::NarrowBandLevelSet, I) -> value

Entry point for value lookup on a NarrowBandLevelSet at index I, which is assumed to be inside the grid (out-of-grid indices are handled by _getindexbc before reaching this function).

Tries the dict first; if I is not stored (i.e. it is inside the grid but outside the narrow band), falls back to _extrapolate_nb_rec to approximate the value from nearby band nodes. Throws an error if no path to stored values can be found.

_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.

_clear_buffer!(ϕ::MeshField)

Clear the active entries of a buffer before it is used as a write target in a time-integration step. For a NarrowBandLevelSet this empties the values dict so that stale entries from a previous band do not survive the src/dst swap. For a full-domain field this is a no-op.

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

Find the point on the zero level-set of 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.

_closest_point_on_interface(sdf, x)

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

If the first solve does not converge (e.g. because the closest point lies on a neighbouring polynomial patch), a single retry is attempted using the best iterate from the failed solve as a new seed on its own patch.

_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 │.

_ensure_boundary_conditions(ϕ)

Return ϕ unchanged if it already carries boundary conditions, otherwise wrap it with LinearExtrapolationBC on every face.

_extrapolate_nb_rec(nb::NarrowBandLevelSet, I, max_dim) -> value or nothing

Approximate the value at an in-grid index I that is not stored in the band dict, by bilinear/trilinear (degree-1) Lagrange extrapolation from nearby band values.

The algorithm processes dimensions 1 through max_dim in order. For each dimension, it searches outward from I (nearest first, both sides) for an anchor point in the dict and uses the two consecutive nodes starting there (anchor and anchor+step) as a linear stencil.

Stencil values are resolved by calling _extrapolate_nb_rec recursively with max_dim = dim - 1, so each stencil point can itself be extrapolated using lower dimensions. This produces a tensor-product extrapolation that handles indices outside the band in multiple dimensions simultaneously.

Returns nothing if no dimension yields a valid stencil.

Why degree 1 and no higher?

Outside the band, |ϕ| ≥ halfwidth. The only property we need from extrapolated ghost values is sign correctness (no spurious zeros). Linear extrapolation from a well-conditioned SDF preserves sign as long as halfwidth > Δx, which any reasonable band satisfies. Quadratic or higher extrapolation introduces polynomial oscillations that can flip the sign, creating false interface crossings that corrupt NewtonSDF and cause the band to migrate or collapse.

_lagrange_extrap_from(nb, I, dim, anchor, side, k, P) -> value or nothing

Attempt to evaluate a degree-P Lagrange extrapolant at I[dim] using P+1 consecutive stencil nodes along dimension dim:

anchor, anchor + side, anchor + 2*side, …, anchor + P*side

The stencil extends from anchor deeper into the band (away from I). The target I[dim] is at distance k from anchor in the opposite direction, corresponding to local coordinate ξ = -k relative to stencil nodes at ξ = 0, 1, …, P. This matches the convention of _lagrange_extrap_weight(j, k, P).

Each stencil value is resolved via _extrapolate_nb_rec(nb, Ij, dim - 1), using only dimensions lower than dim. Returns nothing if any stencil point falls outside the grid or cannot be resolved.

_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

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

_level_set_extrema(src) -> (ϕ_min, ϕ_max)

Compute the extrema of the level-set values. Handles both Array and Dict storage.

_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.

_project_to_interface(p, x_start, maxiters, ftol, safeguard_dist)

Use Newton's method to project a starting point onto the zero level set of p. Returns the converged point or nothing if Newton fails to converge or if the iterate moves more than safeguard_dist from x_start.

_push_record!(log, tc, t_step, reinit_time, did_reinit, update_times, compute_times, ϕ_min, ϕ_max)

Add a new StepRecord to log.

_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.

_show_fields(io, ϕ::MeshField; prefix="  ")

Print the fields of ϕ as indented tree lines: grid info (via _show_fields for CartesianGrid), boundary conditions, narrow-band info (if applicable), element type, and value range (for real-valued fields).

_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 → "²".

_timed_reinit!(ϕ, reinit, nsteps) -> (elapsed, did_reinit)

Run reinitialize! if due at this step, returning the elapsed time and whether it ran.

_timed_update_terms!(terms, ϕ, t, update_times)

Run update_term! for each term, accumulating the elapsed time per term into update_times.

active_cells(ϕ)

Return the cell indices that are active for interface sampling. A cell is active if all its corners are active indices. For a full-grid level set, this is all cells; for a narrow band, only cells where all corners are within the band. Dispatch point for NewtonSDF construction.

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.

boundary_conditions(eq::LevelSetEquation)

Return the boundary conditions of the equation.

cell_extrema(itp::PiecewisePolynomialInterpolant, I::CartesianIndex)

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

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.

circle(grid; center = (0, 0), radius = 1)

Create a 2D circle with the specified center and radius on a grid. Returns a LevelSet field.

complement!(ϕ::LevelSet)

In-place complement of a level set (negates the values).

complement(ϕ::LevelSet)

Return the complement of a level set (negates the values).

compute_cfl(terms, ϕ::LevelSet, 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(itp::PiecewisePolynomialInterpolant, x)

Compute the multi-index of the cell containing point x. If x is outside the domain, the index of the closest cell is returned (clamping to the boundary).

current_state(eq::LevelSetEquation)

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

current_time(eq::LevelSetEquation)

Return the current time of the simulation.

curvature(ϕ::LevelSet, I)

Compute the mean curvature of ϕ at I using κ = ∇ ⋅ (∇ϕ / |∇ϕ|). We use the formula κ = (Δϕ |∇ϕ|^2 - ∇ϕ^T Hϕ ∇ϕ) / |∇ϕ|^3 with first order finite differences. https://en.wikipedia.org/wiki/Meancurvature#Implicitformofmean_curvature

curvature(ϕ::LevelSet)

Compute the mean curvature of ϕ at I using κ = ∇ ⋅ (∇ϕ / |∇ϕ|). See curvature(ϕ::LevelSet, I) for more details.

using LevelSetMethods
N = 50
grid = CartesianGrid((-1, -1), (1, 1), (N, N))
ϕ = LevelSetMethods.star(grid)
using GLMakie
coeff = exp.(-40.0 * values(ϕ) .^ 2)
κ = curvature(ϕ) .* coeff
xs = LevelSetMethods.grid1d(grid, 1)
ys = LevelSetMethods.grid1d(grid, 2)
fig, ax, hm = heatmap(xs, ys, κ)
Colorbar(fig[:, end+1], hm)
contour!(xs, ys, values(ϕ); levels = [0.0])

dumbbell(grid; width = 1, height = 0.2, radius = 0.25, center = (0, 0))

Create a 2D dumbbell shape consisting of two circles connected by a rectangle. Returns a LevelSet field.

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

Export a surface mesh of the 3D interface (where ϕ = 0) to filename. Requires the MMG extension to be loaded.

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

Export a 3D volume mesh of the interior domain (where ϕ < 0) to filename. Requires the MMG extension to be loaded.

extend_along_normals!(F, ϕ::LevelSet;
                      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!(itp::PiecewisePolynomialInterpolant, base_idxs::CartesianIndex)

Fill the internal buffer of itp with the Bernstein coefficients for the cell at base_idxs.

get_sample_points(sdf::NewtonSDF)

Return the interface sample points used to build the KDTree of sdf.

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).

grad_norm(ϕ::LevelSet)

Compute the norm of the gradient of ϕ, i.e. |∇ϕ|, at all grid points.

gradient(p, x)

Compute the gradient of p at point x using ForwardDiff.

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

Compute the gradient vector $∇ϕ$ at grid index I using centered finite differences. Returns an SVector (or Vector) of derivatives.

gradient(ϕ::LevelSet)

Compute the gradient vector $∇ϕ$ for all grid points.

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.

hessian(p, x)

Compute the hessian of p at point x using ForwardDiff.

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

Compute the Hessian matrix $\mathbf{H}ϕ = ∇∇ϕ$ at grid index I using second-order finite differences. Returns a Symmetric matrix.

hessian(ϕ::LevelSet)

Compute the Hessian matrix for all grid points.

integrate!(ls::LevelSetEquation,tf,Δt=Inf)

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.

interpolate(ϕ::MeshField, order::Int = 3)

Create a piecewise polynomial interpolant of the given order for ϕ.

A deep copy of ϕ is made so that the interpolant is independent of future modifications to ϕ. If ϕ has no boundary conditions, ExtrapolationBC{2} is added automatically on all sides (this is necessary to evaluate the interpolant near the boundary, where the stencil may require out-of-bounds values).

The returned object itp behaves like a function and supports:

• itp(x): evaluate the interpolant at point x
• make_interpolant(itp, I): return the local interpolant for cell I
• gradient(itp, x): gradient at x (via make_interpolant + ForwardDiff)
• hessian(itp, x): hessian at x (via make_interpolant + ForwardDiff)
• cell_extrema(itp, I): lower and upper bounds of the interpolant in cell I

To avoid unnecessary copying, call update!(itp, ϕ) to update the interpolant's internal field with new values from ϕ while reusing the existing interpolation matrix and scratch buffers; see update! for details.

make_interpolant(itp::PiecewisePolynomialInterpolant, I::CartesianIndex)

Create a BernsteinPolynomial for the cell at multi-index I.

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.

nodeindices(g::CartesianGrid)

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

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

Compute the unit exterior normal vector $\mathbf{n} = \frac{∇ϕ}{\|∇ϕ\|}$ at grid index I.

normal(ϕ::LevelSet)

Compute the unit exterior normal vector $\mathbf{n} = \frac{∇ϕ}{\|∇ϕ\|}$ for all grid points.

using LevelSetMethods
N = 50
grid = CartesianGrid((-1, -1), (1, 1), (N, N))
ϕ = LevelSetMethods.star(grid)
using GLMakie
n = normal(ϕ)
xs = LevelSetMethods.grid1d(grid, 1)
ys = LevelSetMethods.grid1d(grid, 2)
coeff = exp.(-40.0 * values(ϕ) .^ 2)
us = getindex.(n, 1) .* coeff
vs = getindex.(n, 2) .* coeff
arrows(xs, ys, us, vs; arrowsize = 10 * vec(coeff), lengthscale = 2.0 / (N - 1))
contour!(xs, ys, values(ϕ); levels = [0.0])

perimeter(ϕ::LevelSet)

Compute the perimeter area of the level-set function.

Note: this function does not compute the perimeter on the borders of the domain.

using LevelSetMethods
R = 0.5
S0 = 2π * R
grid = CartesianGrid((-1, -1), (1, 1), (200, 200))
ϕ = LevelSetMethods.circle(grid; center = (0, 0), radius = R)
LevelSetMethods.perimeter(ϕ), S0

# output

(3.1426415491430366, 3.141592653589793)

proven_empty(itp::PiecewisePolynomialInterpolant, I::CartesianIndex; surface=false)

Return true if the cell I is guaranteed to not contain the interface (if surface=true) or to not contain any part of the interior (if surface=false).

Note that proven_empty being false does not mean the cell is non-empty, but rather that we can't guarantee emptiness based on the convex hull property of the Bernstein basis.

rebuild_band!(nb::NarrowBandLevelSet, sdf)

Rebuild the active node set from the signed distance function sdf using a breadth-first search seeded from all previously active nodes. The BFS expands axis-aligned neighbors and adds a node whenever |sdf(x)| < halfwidth, stopping a branch when a node falls outside the band.

rectangle(grid; center = (0, 0), width = (1, 1))

Create a rectangle (or N-dimensional box) with the specified center and width on a grid. Returns a LevelSet field.

reinitialize!(ϕ::LevelSet, r::NewtonReinitializer)

Reinitialize the level set ϕ to a signed distance function in place using the Newton closest-point method.

reinitialize!(eq::LevelSetEquation)

Reinitialize the current state of the level-set equation using its attached NewtonReinitializer.

reinitializer(eq::LevelSetEquation)

Return the NewtonReinitializer (if any) attached to the equation.

reset_log!(log::SimulationLog, t0 = log.t0)

Clear all records from log and optionally reset the initial time.

reset_log!(eq::LevelSetEquation)

Clear the simulation log of eq, resetting step count and all timing records.

set_makie_theme!()

Set the Makie theme to LevelSetMethods.makie_theme().

sphere(grid; center = (0, 0, 0), radius)

Create a 3D sphere with the specified center and radius on a grid. Returns a LevelSet field.

star(grid; radius = 1, deformation = 0.25, n = 5.0)

Create a 2D star shape defined in polar coordinates by $r = R(1 + d \cos(nθ))$. Returns a LevelSet field.

terms(eq::LevelSetEquation)

Return the tuple of LevelSetTerms of the equation.

time_integrator(eq::LevelSetEquation)

Return the TimeIntegrator of the equation.

update!(sdf::AbstractSignedDistanceFunction, ϕ::LevelSet)

Rebuild sdf in place from the new level set ϕ.

update!(sdf::NewtonSDF, ϕ)

Rebuild sdf in place from the new level set ϕ, reusing the existing interpolant buffers, upsample density, and solver tolerances.

update!(itp::PiecewisePolynomialInterpolant, ϕ)

Copy the values of ϕ into the interpolant's internal field and invalidate the cell cache. This is cheaper than calling interpolate again because it reuses the existing interpolation matrix and scratch buffers.

update_bc!(bc::BoundaryCondition, t)

Update the current time stored in a boundary condition. Only meaningful for DirichletBC; all other BC types are no-ops.

update_bcs!(ϕ::MeshField, t)

Update the current time in all DirichletBC boundary conditions of ϕ. Called automatically by the time-stepper at each stage.

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

Update the internal state of a LevelSetTerm before computing its contribution. This is called at each stage of the time integration.

value_gradient_hessian(p, x)

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

volume(ϕ::LevelSet)

Compute the volume of the level-set function.

using LevelSetMethods
R = 0.5
V0 = π * R^2
grid = CartesianGrid((-1, -1), (1, 1), (200, 200))
ϕ = LevelSetMethods.circle(grid; center = (0, 0), radius = R)
LevelSetMethods.volume(ϕ), V0

# output

(0.7854362890190668, 0.7853981633974483)

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.

zalesak_disk(grid; center = (0, 0), radius = 0.5, width = 0.25, height = 1)

Create a Zalesak disk (a circle with a rectangular slot cut out). Used for testing advection schemes. Returns a LevelSet field.

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

Compute a mesh of the LevelSet ϕ 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_surface_mesh(eq::LevelSetEquation, args...; kwargs...)

Call export_surface_mesh on current_state(eq).

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

Compute a mesh of the domains associated with LevelSet 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.

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

Call export_volume_mesh on current_state(eq).