LevelSetMethods

Installation

To install the library, run the following command on a Julia REPL:

using Pkg; Pkg.add("LevelSetMethods")

This will install the latest tagged version of the package and its dependencies.

For visualization, you may also want to install one of the Makie backends (we suggest GLMakie for 3D plots and animations).

Overview

This package defines a LevelSetEquation type that can be used to solve partial differential equations of the form

\[\phi_t + \underbrace{\boldsymbol{u} \cdot \nabla \phi}_{\substack{\text{advection} \\ \text{term}}} + \underbrace{v |\nabla \phi|}_{\substack{\text{normal} \\ \text{term}}} + \underbrace{b \kappa |\nabla \phi|}_{\substack{\text{curvature} \\ \text{term}}} + \underbrace{\text{sign}(\phi)(|\nabla \phi| - 1)}_{\substack{\text{reinitialization}\\ \text{term}}} = 0\]

where

  • $\phi : \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}$ is the level set function
  • $\boldsymbol{u} :\mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}^d$ is a given (external) velocity field
  • $v : \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}$ is a normal speed
  • $b : \mathbb{R}^d \times \mathbb{R}^+ \to \mathbb{R}$ is a function that multiplies the curvature $\kappa = \nabla \cdot (\nabla \phi / |\nabla \phi|)$

Here is how it looks in practice to create a simple LevelSetEquation:

using LevelSetMethods
grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
# ϕ    = LevelSet(x -> sqrt(2*x[1]^2 + x[2]^2) - 1/2, grid) # a disk
ϕ    = LevelSetMethods.dumbbell(grid) # a predefined shape
𝐮    = (x,t) -> (-x[2], x[1])
eq   = LevelSetEquation(;
  terms = (AdvectionTerm(𝐮),),
  levelset = ϕ,
  bc = PeriodicBC()
)
Level-set equation given by

 	 ϕₜ + 𝐮 ⋅ ∇ ϕ = 0

Current time 0.0

You can easily plot the current state of your level set equation using the plot function from Makie:

using GLMakie # loads the MakieExt from LevelSetMethods
LevelSetMethods.set_makie_theme!() # optional theme customization
plot(ϕ)
Example block output

To step it in time, we can use the integrate! function:

integrate!(eq, 1)
Level-set equation given by

 	 ϕₜ + 𝐮 ⋅ ∇ ϕ = 0

Current time 1.0

This will advance the solution up to t = 1, modifying ϕ in the process:

plot(ϕ)
Example block output

Creating an animation can be achieved by calling integrate! in a loop and saving the results to a file:

using GLMakie
theme = LevelSetMethods.makie_theme()
anim = with_theme(theme) do
    obs = Observable(eq)
    fig = Figure()
    ax = Axis(fig[1, 1])
    plot!(ax, obs)
    framerate = 30
    t0 = current_time(eq)
    tf = t0 + π
    timestamps = range(t0, tf; step = 1 / framerate)
    record(fig, joinpath(@__DIR__, "ls_intro.gif"), timestamps) do t_
        integrate!(eq, t_)
        return obs[] = eq
    end
end
"/home/runner/work/LevelSetMethods.jl/LevelSetMethods.jl/docs/build/ls_intro.gif"

Here is what the .gif file looks like:

Dumbbell

For more interesting applications and advanced usage, see the examples section!

Other resources

There is an almost one-to-one correspondance between each of the LevelSetTerms described above and individual chapters of the book by Osher and Fedwick on level set methods [1], so users interested in digging deeper into the theory/algorithms are encourage to consult that refenrence. We also drew some inspiration from the great Matlab library ToolboxLS by Ian Mitchell [2].

Going further

As illustrated above, the LevelSetEquation type is the main structure of this package. Becoming familiar with its fields and methods is a good starting point to use the package:

LevelSetMethods.LevelSetEquationType
LevelSetEquation(; terms, levelset, boundary_conditions, t = 0, integrator = RK3())

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.

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
source

To learn more about the package, you should also check out the following sections:

  • The section on terms for a detailed description of each term and their corresponding customizations
  • The section on time integrators for a description of the available time integrators and how to use them
  • The section on boundary conditions for a description of the available boundary conditions and how to use them

Finally, the examples section contains a list of examples that demonstrate some hopefully cool applications.

Bibliography

[1]
[2]
[3]
G. Allaire, F. Jouve and A.-M. Toader. Structural optimization using sensitivity analysis and a level-set method. Journal of computational physics 194, 363–393 (2004).