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

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.

LevelSet

Alias for MeshField with vals as an AbstractArray of Reals.

LevelSetEquation(; terms, levelset, boundary_conditions, t = 0, integrator = RK2(),
reinit_freq = 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 can be performed every reinit_freq time steps. By default, no 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.

struct NeumannBC <: BoundaryCondition

Homogenous Neumann boundary condition, i.e. ∂x ϕ = 0.

struct NeumannGradientBC <: BoundaryCondition

Homogenous Neumann gradient boundary condition, i.e. ∂xx ϕ = 0.

struct NormalMotionTerm{V,M} <: 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.

struct PeriodicBC <: BoundaryCondition

Singleton type representing periodic boundary conditions.

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.

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.

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

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.

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)

reinitialize!(ϕ::LevelSet; upsample=2, maxiters=10, xtol=1e-8)

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.

Keyword Arguments

• upsample: number of samples to take in each cell when sampling the interface. Higher values yield better initial guesses for the closest point search, but increase the computational cost.
• maxiters: maximum number of iterations to use in the Newton's method to find the closest point on the interface.
• xtol: convergence tolerance for the Newton's method. The iterations stop when the change in position is less than xtol.

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