LevelSetEvolution

In GridapTopOpt, the a level-set function is evolved and reinitialised as a signed distance function using an Evolver and Reinitialiser implementations, respectively. These are wrapped in an object LevelSetEvolution for the purpose of optimisation.

List of Evolver types:

MethodAmbient mesh typeCachedPreferred Method
FiniteDifferenceEvolverCartesian
CutFEMEvolverUnstructured✗*

*: Caching is disabled due to a bug. This is will be improved in future.

List of Reinitialiser types:

MethodAmbient mesh typeCachedDifferentiablePreferred Method
FiniteDifferenceReinitialiserCartesian
StabilisedReinitialiserUnstructured
HeatReinitialiserUnstructured
Note

The finite difference methods are based on the upwinding scheme of Osher and Fedkiw (link). These are extremely efficient as they require no linear solves. They are, however, currently restricted to meshes built from CartesianDiscreteModel. Note as well that the FiniteDifferenceReinitialiser is can result in slight movement of the zero isosurface.

LevelSetEvolution

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

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

source
GridapTopOpt.reinit!Method
reinit!(::LevelSetEvolution,φ,args...)

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

source

Level-Set Evolvers

FiniteDifferenceEvolver

GridapTopOpt.FiniteDifferenceEvolverType
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
source
GridapTopOpt.FiniteDifferenceEvolverMethod
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.
source

CutFEMEvolver

GridapTopOpt.CutFEMEvolverType
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
Warning

Caching for the CutFEMEvolver method is currently disabled. This will be re-enabled in the future."

source
GridapTopOpt.CutFEMEvolverMethod
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!.
source

Level-Set Reinitialisers

FiniteDifferenceReinitialiser

GridapTopOpt.FiniteDifferenceReinitialiserType
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
source
GridapTopOpt.FiniteDifferenceReinitialiserMethod
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.
source

StabilisedReinitialiser

GridapTopOpt.StabilisedReinitialiserType
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.
source
GridapTopOpt.StabilisedReinitialiserMethod
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.
source

HeatReinitialiser

GridapTopOpt.HeatReinitialiserType
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.
source
GridapTopOpt.HeatReinitialiserMethod
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.
Note

When using an unstructured mesh with significant changes in mesh sizes, you should consider testing different values of t, e.g., minimum, maximum, or mean element diameter. It is also possible to set t to be the square of the element size field from get_element_diameter_field: E.g.,

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

We have found that the latter yields better results for some problems.

source

IdentityReinitialiser

Your own custom Evolver

GridapTopOpt.EvolverType
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.

source
GridapTopOpt.evolve!Method
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
source
GridapTopOpt.evolve!Method
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
source
GridapTopOpt.evolve!Method
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
source
GridapTopOpt.evolve!Method
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
source

Your own custom Reinitialiser

GridapTopOpt.ReinitialiserType
abstract type Reinitialiser

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

source
GridapTopOpt.reinit!Method
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
source
GridapTopOpt.reinit!Method
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
source
GridapTopOpt.reinit!Method
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
source
GridapTopOpt.reinit!Method
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
source

Custom Stencil for finite difference methods

GridapTopOpt.StencilType
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.

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

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

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

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

source
GridapTopOpt.check_orderFunction
check_order(::Stencil,order)

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

source