Reference

abstract type AbstractMesh{N,T}

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

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.

See also berninterp.

abstract type BoundaryCondition

Types used to specify boundary conditions.

CartesianGrid(lc, hc, n)

Create a uniform cartesian grid with lower corner lc, upper corner hc and 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{2, Int64}([0, 0], [1, 1], (10, 4))

const CartesianMeshField{V,M<:CartesianGrid} = MeshField{V,M}

MeshField over a CartesianGrid.

struct CurvatureTerm{V,M} <: LevelSetTerm

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

struct DirichletBC{T} <: BoundaryCondition

A Dirichlet boundary condition taking values of f(x) at the boundary.

struct ExtrapolationBC{P} <: BoundaryCondition

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

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

using LevelSetMethods
ExtrapolationBC(5)

# output

ExtrapolationBC{5}()

LevelSet

Alias for MeshField with vals as an AbstractArray of Reals.

LevelSetEquation(; terms, levelset, boundary_conditions, t = 0, integrator = RK2(),
reinit = nothing)

Create a of a level-set equation of the form ϕₜ + sum(terms) = 0, where each t ∈ terms is a LevelSetTerm and levelset is the initial LevelSet.

Calling integrate!(ls, tf) will evolve the level-set equation up to time tf, modifying the current_state(eq) and current_time(eq) of the object eq in the process (and therefore the original levelset).

Boundary conditions can be specified in two ways. 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 of the level-set function is controlled by the reinit parameter, a Reinitializer object that specifies both the reinitialization algorithm parameters and the frequency (in time steps). By default, no automatic reinitialization is performed. Using this feature requires the ReinitializationExt to be loaded.

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 and curvature terms
bc = PeriodicBC()                                   # periodic boundary conditions
eq = LevelSetEquation(; terms, levelset = ϕ, bc)    # level-set equation

# output

Level-set equation given by

 	 ϕₜ + 𝐮 ⋅ ∇ ϕ = 0

Current time 0.0

abstract type LevelSetTerm

A typical term in a level-set evolution equation.

struct MeshField{V,M,B}

A field described by its discrete values on a mesh.

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

MeshField(f::Function, m)

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

NeumannBC = ExtrapolationBC{1}

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

NeumannGradientBC = ExtrapolationBC{2}

Homogeneous Neumann gradient boundary condition (∂²ϕ/∂n² = 0). Alias for ExtrapolationBC{2}: linear extrapolation into ghost cells.

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.

struct PiecewisePolynomialInterpolation{Φ, N, T}

A piecewise polynomial interpolant for a MeshField.

See interpolate.

struct RK2

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

struct RK3

Third order total variation dimishing Runge-Kutta scheme.

struct ReinitializationTerm <: LevelSetTerm

Level-set term representing sign(ϕ) (|∇ϕ| - 1). This LevelSetTerm should be used for reinitializing the level set into a signed distance function: for a sufficiently large number of time steps this term allows one to solve the Eikonal equation |∇ϕ| = 1.

There are two ways of constructing a ReinitializationTerm:

• using ReinitializationTerm(ϕ₀::LevelSet) precomputes the sign term on the initial level set ϕ₀, as in equation 7.5 of Osher and Fedkiw;
• using ReinitializationTerm() constructs a term that computes the sign term on-the-fly at each time step, as in equation 7.6 of Osher and Fedkiw.

Reinitializer(; upsample=8, maxiters=20, xtol=1e-8, ftol=1e-8, reinit_freq=1)

Configuration for Newton-based reinitialization to a signed distance function. The reinit_freq parameter specifies how often reinitialization is performed (in time steps).

struct SemiImplicitI2OE

Semi-implicit finite-volume scheme of the I2OE family (Mikula et al.) for advection problems.

abstract type TimeIntegrator end

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

D2(ϕ::CartesianMeshField,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⁰(ϕ::CartesianMeshField,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⁺⁺(ϕ::CartesianMeshField,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⁻⁺(ϕ::CartesianMeshField,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⁰(ϕ::CartesianMeshField,I,dim)

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

D⁺(ϕ::CartesianMeshField,I,dim)

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

D⁻(ϕ::CartesianMeshField,I,dim)

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

_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 P nodes total.

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

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

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

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

compute_index(itp::PiecewisePolynomialInterpolation, x)

Compute the multi-index of the cell containing point x.

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

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::PiecewisePolynomialInterpolation, base_idxs::CartesianIndex)

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

grad_norm(ϕ::LevelSet[, I])

Compute the norm of the gradient of ϕ at index I, i.e. |∇ϕ|, or for all grid points if I is not provided.

gradient(ϕ::LevelSet, I)

Return the gradient vector ∇ϕ of ϕ at I

hessian(ϕ::LevelSet, I)

Return the Hessian matrix Hϕ of ϕ at I

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{3} is added automatically on all sides.

The returned object itp behaves like a function and supports:

• itp(x): evaluate the interpolant at point x
• gradient(itp, x): gradient at x (via make_interpolant + ForwardDiff)
• make_interpolant(itp, I): return the local BernsteinPolynomial for cell I, which itself supports p(x), gradient(p, x), value_and_gradient(p, x), value_gradient_hessian(p, x)
• cell_extrema(itp, I): min/max of the interpolant over cell I
• proven_empty(itp, I): check whether cell I can be skipped

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

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

makie_theme()

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

normal(ϕ::LevelSet, I)

Compute the unit exterior normal vector of ϕ at I using n = ∇ϕ/|∇ϕ|

normal(ϕ::LevelSet)

Compute the unit exterior normal vector of ϕ using n = ∇ϕ/|∇ϕ|

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::PiecewisePolynomialInterpolation, 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).

reinitialize!(ϕ::LevelSet, reinitializer = Reinitializer())

Reinitializes the level set ϕ to a signed distance, modifying it in place.

The method works by first sampling the zero-level set of the interface, and then for each grid point, finding the closest point on the interface using a Newton-based method. The distance to the closest point is then used as the new value of the level set at that grid point, with the sign determined by the original level set value. See [4] for more details.

Arguments

• ϕ: The level set to reinitialize.
• reinitializer: Configuration for the reinitialization. Defaults to Reinitializer(). See Reinitializer for details.

set_makie_theme!()

Set the Makie theme to LevelSetMethods.makie_theme().

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

Called before computing the term at each stage of the time evolution.

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)

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