LevelSetEvolution

struct LevelSetEvolution{A,B} <: AbstractLevelSetEvolution

A wrapper to hold a level-set evolver and reinitialiser.

evolve!(s::LevelSetEvolution,φ,args...)

Evolve the level set function φ using the evolver in LevelSetEvolution.

get_evolver(s::LevelSetEvolution)

Return the level-set function evolver.

get_ls_space(s::LevelSetEvolution)

Return the finite element space used for the level-set function.

get_min_dof_spacing(s::LevelSetEvolution)

Return the minimum spacing of DOFs for the level-set function.

reinitialiser(s::LevelSetEvolution)

Return the level-set function reinitialiser.

reinit!(::LevelSetEvolution,φ,args...)

Reinitialise the level set function φ using the reinitialiser in LevelSetEvolution.

struct FiniteDifferenceEvolver{O}

A standard forward Euler in time finite difference method for solving the Hamilton-Jacobi evolution equation on order O finite elements in serial or parallel.

Based on the scheme by Osher and Fedkiw (link).

Hamilton-Jacobi evolution equation

$\frac{\partial\phi}{\partial t} + V(\boldsymbol{x})\lVert\boldsymbol{\nabla}\phi\rVert = 0,$

with $\phi(0,\boldsymbol{x})=\phi_0(\boldsymbol{x})$ and $\boldsymbol{x}\in D,~t\in(0,T)$.

Parameters

• stencil::Stencil: Spatial finite difference stencil for a single step HJ equation and reinitialisation equation.
• model: A CartesianDiscreteModel.
• space: FE space for level-set function
• perm: A permutation vector
• params: Tuple of additional params

FiniteDifferenceEvolver(stencil::Stencil,model,space,max_steps)

Create an instance of FiniteDifferenceEvolver given a stencil, model, FE space, and additional optional arguments. This automatically creates the DoF permutation to handle high-order finite elements.

Optional Arguments

• max_steps: number of timesteps
• correct_ls: Boolean for whether or not to ensure LS DOFs aren't zero this MUST be true for differentiation in unfitted methods.

struct FirstOrderStencil{D,T} <: Stencil end

A first order upwind difference scheme based on Osher and Fedkiw (link).

mutable struct CutFEMEvolver{V,M} <: Evolver

CutFEM method for level-set evolution based on method developed by

• Burman et al. (2018). DOI: 10.1016/j.cma.2017.09.005.
• Burman et al. (2017). DOI: 10.1016/j.cma.2016.12.021.
• Burman and Fernández (2009). DOI: 10.1016/j.cma.2009.02.011

This solves the tranport equation $\frac{\partial\phi(t,oldsymbol{x})}{\partial t}+\boldsymbol{\beta}\cdot\boldsymbol{\nabla}\phi(t,\boldsymbol{x})=0,$ with $boldsymbol{\beta}=\boldsymbol{n}v_h$, $\phi(0,\boldsymbol{x})=\phi_0(\boldsymbol{x}),$ and $\quad\boldsymbol{x}\in D,~t\in(0,T)$.

Parameters

• ode_solver::ODESolver: ODE solver
• Ωs::B: EmbeddedCollection holding updatable triangulation and measures from GridapEmbedded
• dΩ_bg::C: Measure for integration
• space::B: Level-set FE space
• assembler::Assembler: FE assembler
• params::D: Tuple of stabilisation parameter γg, mesh sizes h, and max steps max_steps, and background mesh skeleton parameters

CutFEMEvolver(V_φ::B,Ωs::EmbeddedCollection,dΩ_bg::A,h;
  max_steps=10,
  γg = 0.1,
  ode_ls = LUSolver(),
  ode_nl = ode_ls,
  ode_solver = MutableRungeKutta(ode_nl, ode_ls, 0.1, :DIRK_CrankNicolson_2_2),
  assembler=SparseMatrixAssembler(V_φ,V_φ)
) where {A,B}

Create an instance of CutFEMEvolver with the space for the level-set V_φ, the EmbeddedCollection Ωs for the triangulation and measures, the measure dΩ_bg for the background mesh, and the mesh size h. The mesh size h can either be a scalar or a CellField object.

The optional arguments are:

• correct_ls: Boolean for whether or not to ensure LS DOFs aren't zero this MUST be true for differentiation in unfitted methods.
• max_steps: Maximum number of steps for the ODE solver.
• γg: Stabilisation parameter for the continuous interior penalty term.
• ode_ls: Linear solver for the ODE solver.
• ode_nl: Non-linear solver for the ODE solver.
• ode_solver: ODE solver, default is MutableRungeKutta(ode_nl, ode_ls, 0.1, :DIRK_CrankNicolson_2_2).
• assembler: Assembler for the finite element space, default is SparseMatrixAssembler(V_φ,V_φ).

Note

• The stepsize dt = 0.1 in MutableRungeKutta is a place-holder and is updated using the γ passed to evolve!.

struct FiniteDifferenceReinitialiser{O}

A standard forward Euler in time finite difference method for solving the reinitialisation equation on order O finite elements in serial or parallel.

Based on the scheme by Osher and Fedkiw (link).

Reinitialisation equation

$\frac{\partial\phi}{\partial t} + \mathrm{sign}(\phi_0)(\lVert\boldsymbol{\nabla}\phi\rVert-1) = 0,$

with $\phi(0,\boldsymbol{x})=\phi_0(\boldsymbol{x})$ and $\boldsymbol{x}\in D,~t\in(0,T)$.

Parameters

• stencil::Stencil: Spatial finite difference stencil for a single step HJ equation and reinitialisation equation.
• model: A CartesianDiscreteModel.
• space: FE space for level-set function
• perm: A permutation vector
• params: Tuple of additional params

FiniteDifferenceReinitialiser(stencil::Stencil,model,space;γ_reinit,tol,max_steps)

Create an instance of FiniteDifferenceReinitialiser given a stencil, model, FE space, and additional optional arguments. This automatically creates the DoF permutation to handle high-order finite elements.

Optional Arguments

• γ_reinit: coeffient on the time step size.
• tol: stopping tolerance for reinitialiser
• max_steps: number of timesteps
• correct_ls: Boolean for whether or not to ensure LS DOFs aren't zero this MUST be true for differentiation in unfitted methods.

struct StabilisedReinitialiser{V,M} <: Reinitialiser

Stabilised FE method for level-set reinitialisation. Artificial viscosity approach (ArtificialViscosity) based on Mallon et al. (2023). DOI: 10.1016/j.cma.2025.118203. Interior jump penalty approach (InteriorPenalty) adapted from that work and replaces the artifical viscosity term with an interior jump penalty term.

Parameters

• nls::NonlinearSolver: Nonlinear solver for solving the reinitialisation equation
• stabilisation_method::A: A StabilisationMethod method for stabilising the problem
• Ωs::B: EmbeddedCollection holding updatable triangulation and measures from GridapEmbedded
• dΩ_bg::D: Background mesh measure
• space::C: FE space for level-set function
• assembler::Assembler: Assembler for LS FE space
• params::E: Tuple of Nitsche parameter γd and mesh size h

Note

• We expect the EmbeddedCollection Ωs to contain :dΓ. If this is not available we add it to the recipe list in Ωs and a warning will appear.

StabilisedReinitialiser(V_φ::C,Ωs::EmbeddedCollection,dΩ_bg::B,h;
  γd = 20.0,
  nls = NewtonSolver(LUSolver();maxiter=20,rtol=1.e-14,verbose=true),
  assembler=SparseMatrixAssembler(V_φ,V_φ),
  stabilisation_method::A = InteriorPenalty(V_φ)
) where {A,B,C}

Create an instance of StabilisedReinitialiser with the space for the level-set V_φ, the EmbeddedCollection Ωs for the triangulation and measures, the measure dΩ_bg for the background mesh, and the mesh size h. The mesh size h can either be a scalar or a CellField object.

The optional arguments are:

• correct_ls: Boolean for whether or not to ensure LS DOFs aren't zero this MUST be true for differentiation in unfitted methods.
• γd: Interface penalty parameter for the reinitialisation equation.
• nls: Nonlinear solver for solving the reinitialisation equation.
• assembler: Assembler for the finite element space.
• stabilisation_method: A StabilisationMethod method for stabilising the problem.

struct HeatReinitialiser <: Reinitialiser end

A heat method for reinitialisation of the level-set function as an approximate signed distance function. This is based on Feng and Crane (2024) [doi: 10.1145/3658220].

This works using in three steps:

• Step 1: Solve (I - tΔ)X = {Γ_n on Γ, 0 otherwise} for small t on M
• Step 2: Get Y = X/|X| (L2 projection to ensure we have an FEFunction)
• Step 3: Solve ΔΦ = ∇ ⋅ Y on M with n ⋅ ∇Φ = ∂M_n ⋅ Y on ∂M.

The resulting Φ is an approximation of the signed distance function. In the above M is our background domain and Γ is the interface defined by the zero isosurface of the level-set function.

Note, in Step 3 we add a penalty term to ensure that the resulting isosurface Φ=0 lies close to φ=0. This results in a multifield tangent space for yφtosdf.

The implementation is in terms of AffineFEStateMaps to enable AD through the maps using Zygote. The differentiable map φ_to_sdf can be obtained using the function get_lsf_to_sdf_map.

Parameters

• φ_to_x: AffineFEStateMap mapping from φ to the vector field x (Step 1)
• x_to_y: AffineFEStateMap mapping from the vector field x to the vector field y (Step 2)
• yφ_to_sdf: AffineFEStateMap mapping from the vector field y and φ to the sdf (Step 3)
• Ωs: An EmbeddedCollection for holding information regarding the geometry.

HeatReinitialiser(V_φ,model;
  t  = minimum(get_element_diameters(model))^2,
  γd = 10.0,
  boundary_tags = "boundary",
  V_xy = build_V_xy(model,V_φ),
  M_xyφ = MultiFieldFESpace([Vxy,V_φ]),
  V_sdf = build_V_sdf(model,V_φ),
  Ωs = build_Ωs(model,boundary_tags,V_φ)
)

Create an instance of HeatReinitialiser.

Optional parameters:

• t: Time step for the heat equation, default is the square of the minimum element diameter.
• γd: Penalty parameter to preserve level sets, default is 10.
• boundary_tags: Tag/s for the boundary, default is "boundary" (e.g., for CartesianDiscreteModel). For your own mesh (e.g., unstructured) covering M, you should provide the tags that consistute ∂M.

Advanced optional parameters:

• V_xy: TestFESpace for the vector field x, default is built using build_V_xy.
• M_xyφ: MultiFieldFESpace for the vector field x and φ, default is built using MultiFieldFESpace([V_xy,V_φ]). This is used in the third step to solve for the signed distance function.
• V_sdf: TestFESpace for the signed distance function, default is built using build_V_sdf.
• Ωs: An EmbeddedCollection for holding information regarding the geometry, default is built using build_Ωs
• φ_to_x_ls: Solver for the φtox map, default is LUSolver().
• x_to_y_ls: Solver for the xtoy map, default is LUSolver().
• yφ_to_sdf_ls: Solver for the yφtosdf map, default is LUSolver().
• φ_to_x_adjoint_ls: Adjoint solver for the φtox map, default is φ_to_x_ls.
• x_to_y_adjoint_ls: Adjoint solver for the xtoy map, default is x_to_y_ls.
• yφ_to_sdf_adjoint_ls: Adjoint solver for the yφtosdf map, default is yφ_to_sdf_ls.

hsq(h) = h^2
t = hsq ∘ get_element_diameter_field(model)

struct IdentityReinitialiser <: Reinitialiser end

A level-set function reinitialiser that does nothing.

abstract type Evolver

Your own level-set evolution method can be created by implementing concrete functionality for solve!, get_dof_spacing, and get_ls_space.

evolve!(::Evolver,φ::AbstractVector,vel::AbstractVector,γ,cache)

Evolve the level-set function φ using a velocity field vel and parameter γ according to an Evolution method. Reuse the supplied cache.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the evolver

evolve!(::Evolver,φ::AbstractVector,vel::AbstractVector,γ)

Evolve the level-set function φ using a velocity field vel and parameter γ according to an Evolution method.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the evolver

evolve!(::Evolver,φh,velh,γ)

Evolve the level-set function φh using a velocity field velh and parameter γ according to an Evolution method.

φh and velh should be FEFunctions.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the evolver

evolve!(::Evolver,φh,velh,γ,cache)

Evolve the level-set function φh using a velocity field velh and parameter γ according to an Evolution method. Reuse the supplied cache.

φh and velh should be FEFunctions.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the evolver

get_ls_space(m::Evolver)

Return the finite element space used for the level-set function.

get_min_dof_spacing(m::Evolver)

Return the minimum spacing of DOFs for the level-set function.

abstract type Reinitialiser

Your own level-set reinitialisation method can be created by implementing concrete functionality for solve!.

reinit!(::Reinitialiser,φ::AbstractVector,cache)

Reinitialise the level-set function φ according to an Reinitialiser method. Reuse the supplied cache.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the reinitialiser

reinit!(::Reinitialiser,φ::AbstractVector)

Reinitialise the level-set function φ according to an Reinitialiser method.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the reinitialiser

reinit!(::Reinitialiser,φh,cache)

Reinitialise the level-set function φh according to an Reinitialiser method. Reuse the supplied cache.

φh should be an FEFunction.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the reinitialiser

reinit!(::Reinitialiser,φh)

Reinitialise the level-set function φh according to an Reinitialiser method.

φh should be an FEFunction.

Returns

• φ: The updated level-set function as an AbstractVector
• cache: The cache for the reinitialiser

abstract type Stencil

Finite difference stencil for a single step of the Hamilton-Jacobi evolution equation and reinitialisation equation.

Your own spatial stencil can be implemented by extending the methods below.

evolve!(::Stencil,φ,vel,Δt,Δx,isperiodic,caches) -> φ

Single finite difference step of the Hamilton-Jacobi evoluation equation for a given Stencil.

reinit!(::Stencil,φ_new,φ,vel,Δt,Δx,isperiodic,caches) -> φ

Single finite difference step of the reinitialisation equation for a given Stencil.

allocate_caches(::Stencil,φ,vel)

Allocate caches for a given Stencil.

check_order(::Stencil,order)

Throw error if insufficient reference element order to implement stencil in parallel.