# A predictor-corrector scheme for approximating signed distances using finite element methods

Amina El Bachari, Johann Rannou, Vladislav Yastrebov, Pierre Kerfriden, Susanne Claus

► To cite this version:

Amina El Bachari, Johann Rannou, Vladislav Yastrebov, Pierre Kerfriden, Susanne Claus. A predictor-corrector scheme for approximating signed distances using finite element methods. 2025. [hal-05122269]

HAL Id: hal-05122269

https://hal.science/hal-05122269v1

Preprint submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A predictor-corrector scheme for approximating signed distances using finite element methods

Amina El Bachari¹,², Johann Rannou³, Vladislav Yastrebov², Pierre Kerfriden², and Susanne Claus¹

¹DTIS, ONERA, Université Paris-Saclay, 91120, Palaiseau, France

²Centre de Matériaux, MINES Paris-PSL, CNRS UMR 7633, Satory, France

³DMAS, ONERA, Université Paris-Saclay, 92320, Châtillon, France

June 20, 2025

## Abstract

In this article, we introduce a finite element method designed for the robust computation of approximate signed distance functions to arbitrary boundaries in two and three dimensions. Our method employs a novel prediction-correction approach, involving first the solution of a linear diffusion-based prediction problem, followed by a nonlinear minimization-based correction problem associated with the Eikonal equation. The prediction step efficiently generates a suitable initial guess, significantly facilitating convergence of the nonlinear correction step.

A key strength of our approach is its ability to handle complex interfaces and initial level set functions with arbitrary steep or flat regions, a notable challenge for existing techniques. Through several representative examples, including classical geometries and more complex shapes such as star domains and three-dimensional tori, we demonstrate the accuracy, efficiency, and robustness of the method, validating its broad applicability for reinitializing diverse level set functions.

## 1 Introduction

The computation of signed distance functions (SDFs) plays a critical role in numerous computational methods, particularly within computational mechanics, computer graphics, and numerical simulations involving interface capturing and tracking [21, 24]. Accurate and efficient methods to compute SDFs are crucial for applications such as shape optimization [2, 13], multi-phase flow [27, 16, 11] or image processing [28, 18]. The signed distance function is typically defined as the solution of the Eikonal equation, a nonlinear partial differential equation whose numerical solution poses significant computational challenges [25].

Traditional methods to compute SDFs generally fall into two main categories: geometric algorithms such as fast marching methods [9, 24, 23, 17] and fast sweeping methods [29, 32, 22] which can be highly efficient but are challenging to parallelize [26], extend to unstructured grids [8], [15] or to higher-order geometric representations. PDE-based methods, which solve variants of the Eikonal equation directly using finite element or finite difference discretizations, offer greater flexibility for complex domains but require careful treatment to ensure stability and accuracy [21].

In this article, we propose a PDE-based finite element approach characterized by a novel prediction-correction strategy. Our method begins with a linear diffusion-based prediction
step, efficiently providing a robust initial approximation of the SDF. This is subsequently refined through a nonlinear correction step formulated as a minimization problem derived from the Eikonal equation. This minimization based formulation of the Eikonal equation has first been used in [18] for the reinitialization of a level set function. As the original method struggles in regions with a gradient norm of the level set function close to zero several modifications of the original problem have been proposed in [6, 31] which modify the equations for small norms of the gradient of the signed distance functions. However, these changes are not consistent with the original problem of solving the Eikonal equation. In this article, we propose to leave the minimization based problem unchanged. Instead, we focus on formulating a prediction problem to provide an ideal initial guess to the minimization based formulation of the Eikonal problem. We will demonstrate that this yields a robust predictor-corrector method that can handle initial functions with arbitrary steep and flat regions for complex two- and three-dimensional geometries.

We demonstrate the use of our method in the context of meshes that fit to the boundary of the complex geometry to which we want to compute the signed distance function as well as the case in which the boundary is immersed in a background finite element grid. In the case of immersed geometries, we reconstruct the interface and enforce Dirichlet conditions weakly on the unfitted interface using cut integration [7]. The implementation of our method is done in FEniCSx [5] and CutFEMx [10].

We validate our method through extensive numerical experimentation, including classical benchmarks such as circles and tori, as well as challenging test cases featuring noisy data and non-smooth boundaries. Our results highlight the accuracy, computational efficiency, and robustness of the approach, demonstrating its potential as a versatile and reliable tool for reinitializing level set functions in diverse applications within computational mechanics and engineering.

The remainder of this article is organized as follows. In Section 2, we introduce the mathematical formulation of the signed distance function and corresponding signed distance problem. In Section 3 we discuss the minimization-based Eikonal problem and introduce its finite element formulation. Section 4 presents the prediction-correction algorithm and its finite element discretization. In Section 5, we validate our predictor-corrector scheme on a range of numerical experiments in both two- and three-dimensions.

# 2 Level set method

## 2.1 Level sets and signed distance functions.

Let's consider fixed domain $D \subset \mathbb{R}^n$ containing implicitly defined domain $\Omega$. A level set function, $\phi: D \to \mathbb{R}$, is used to represent $\Omega \subset D$ such that:

$$ \left\{ \begin{array}{ll} \phi(x) < 0 & \text{if } x \in \Omega, \\ \phi(x) = 0 & \text{if } x \in \Gamma, \\ \phi(x) > 0 & \text{if } x \in D \setminus \{\Omega \cup \Gamma\} \end{array} \right. \qquad (1) $$

with $\Gamma = \{x \in D \text{ such that } x \in \partial\Omega \setminus \partial D\}$.

A level set function with signed distance property with respect to $\phi(x) = 0$ is defined as:

$$ \phi(x) = \begin{cases} -d(x, \Gamma) & \text{if } x \in \Omega, \\ d(x, \Gamma) & \text{if } x \in D \setminus \{\Omega \cup \Gamma\}, \end{cases} \qquad (2) $$

where $d$ is the euclidean distance function defined as:

$$ d(x, \Gamma) = \inf_{y \in \Gamma} d(x, y). \qquad (3) $$
If the domain $\Omega \subset \mathbb{R}^n$ has a piecewise smooth boundary, then the signed distance function to domain boundary $\Gamma$ is differentiable almost everywhere, and the equivalent signed distance property (SDP) holds

$$ |\nabla \phi| = 1, \text{ almost every where in } D, \qquad (4) $$

where $|\cdot|$ denotes the Euclidean norm, defined by

$$ |x| = \left( \sum_{i=1}^{n} x_i^2 \right)^{1/2}, $$

with $x = (x_1, x_2, \dots, x_n) \in \mathbb{R}^n$.

The normal to the contour lines of locally smooth level set function is given by

$$ \mathbf{n}_{\phi} = \frac{\nabla \phi}{|\nabla \phi|}, \qquad (5) $$

which can be simplified as $\mathbf{n}_{\phi} = \nabla \phi$, almost every where, for a level set function with signed distance property.

## 2.2 Need for approximate redistancing

In general, the level set method is employed in the context of moving boundary problems. In this framework, an advection equation is used to advect the domain boundary $\Gamma$ in the normal direction $\theta = v \mathbf{n}_{\phi}$. Regardless of the method used for the advection of the level set, the signed distance property $|\nabla \phi| = 1$ is not preserved. Nevertheless, maintaining the signed distance property is essential, as it greatly facilitates the following aspects:

(i) to compute normals on $\Gamma$, $\mathbf{n}_{\phi} = \nabla \phi / |\nabla \phi|$, which becomes inaccurate when $|\nabla \phi| \to 0$;

(ii) for the level set advection, because the term $\theta \cdot \nabla \phi = v |\nabla \phi|$ becomes difficult to stabilize for steep $\phi$;

(iii) to reconstruct the interface by linear interpolation, which becomes inaccurate for flat $\phi$ as too small change in $\phi$ values yields a large difference in determining the location of $\Gamma$.

The use of the level set method for numerous moving boundary problem (shape optimization, fluid simulation, propagation phenomena) has led to active research on reinitialization methods to address the diverse needs associated with each application.

The need for redistancing extends beyond moving boundary value problems. For instance, solving partial differential equations on domains derived from the semantic segmentation of an image can be made more efficient by redefining the domain's contour using the zero level-set, rather than relying on pixel or voxel boundaries (see for example [12]). In this context, the level-set is not introduced as a distance function but rather as a function that is typically constant inside and outside of the domain, and piece-wise polynomial in a band of interface elements. Consequently, the resulting level-set does not possess signed distance properties. Another example arises in mesh generation, where signed distance functions are used to represent the geometry and to guide the construction of the adaptive meshes (see [19]). More broadly, signed distance functions are widely employed in computer graphics (see [14]).
## 2.3 Signed Distance Problem

Starting from a level set function $\phi_0$ that does not have the signed distance property, we seek to determine the corresponding signed distance function denoted $\phi$ that is satisfying the following properties:

$$ \left\{ \begin{array}{ll} |\nabla \phi| = 1 & \text{in } D, \\ \phi(x) = 0 & \text{on } \Gamma, \\ \operatorname{sgn}(\phi) = \operatorname{sgn}(\phi_0) & \text{in } D \setminus \Gamma. \end{array} \right. \qquad (6) $$

To solve this problem, reinitialization methods are used. A wide range of methods are proposed in the literature; in the following section, we will present the used minimization based method.

# 3 Minimization based reinitialization

## 3.1 Least squares minimization

Our redistancing method is inspired by Basting and Kuzmin [6] who formulated the Eikonal equation, $|\nabla \phi| = 1$, as a least squares minimization problem

$$ \min \left( \int_D \frac{1}{2} (|\nabla \phi| - 1)^2 dx \right). \qquad (7) $$

Setting

$$ \mathcal{F}(\phi) = \left( \int_D \frac{1}{2} (|\nabla \phi| - 1)^2 dx \right) \qquad (8) $$

and taking the Gâteaux derivative of $\mathcal{F}(\phi)$ in direction $\psi$, we obtain

$$ d\mathcal{F}(\phi; \psi) = \int_D \frac{1}{2} 2(|\nabla \phi| - 1) \lim_{\tau \to 0} \left( \frac{|\nabla(\phi + \tau\psi)| - |\nabla \phi|}{\tau} \right) dx \qquad (9) $$

$$ = \int_D (|\nabla \phi| - 1) \frac{\nabla \phi}{|\nabla \phi|} \nabla \psi dx \qquad (10) $$

$$ = \int_D \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) \nabla \psi dx \qquad (11) $$

And the minimization problem becomes

$$ \int_D \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) \nabla \psi dx = 0. \qquad (12) $$

We can derive the strong form by integration by parts

$$ - \int_D \nabla \cdot \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) \psi dx + \int_{\partial D} \left( \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) \cdot \mathbf{n} \right) \psi ds = 0, \qquad (13) $$

and setting

$$ \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) \cdot \mathbf{n} = 0 \text{ on } \partial D, \qquad (14) $$

which yields

$$ \left\{ \begin{array}{ll} -\nabla \cdot \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) = 0 & \text{in } D, \\ \left( \nabla \phi - \frac{\nabla \phi}{|\nabla \phi|} \right) \cdot \mathbf{n} = 0 & \text{on } \partial D. \end{array} \right. \qquad (15) $$
Here, **n** denotes the outward normal to the boundary $\partial D$. To obtain the distance to the boundary defined by $\Gamma$, we impose the following Dirichlet condition

$$ \phi = 0 \text{ on } \Gamma. \qquad (16) $$

The solution to this non-linear problem is the unsigned distance to boundary $\Gamma$. This problem is very stiff and therefore very difficult to solve and does not have any notion of a sign of $\phi$. These two challenges can be overcome by reformulating (15) as an iterative scheme with an initial signed solution $\phi^0$ and $\phi^0 = 0$ on $\Gamma$. There are several options to formulate an iterative scheme. One possibility is to add an artificial time term $\frac{\partial \phi}{\partial \tau}$ to (15) and/or to use a fixed point iteration scheme. In this article, we use the most straightforward iterative scheme given by

$$ \left\{ \begin{array}{ll} -\nabla \cdot \nabla \phi^{n+1} & = -\nabla \cdot \frac{\nabla \phi^n}{|\nabla \phi^n|} & \text{in } D \\ \nabla \phi^{n+1} \cdot \mathbf{n} & = \frac{\nabla \phi^n}{|\nabla \phi^n|} \cdot \mathbf{n} & \text{on } \partial D, \end{array} \right. \qquad (17) $$

which is a fixed point Picard iteration scheme that linearizes the problem. Here, $n$ denotes the $n$-th iterate.

All these iterative schemes have one major drawback:

$$ \text{If } |\nabla \phi^n| \to 0, \text{ then } \frac{1}{|\nabla \phi^n|} \to \infty. \qquad (18) $$

If we do not pose any conditions on our initial function $\phi^0$ then none of the iterative schemes for problem (15) will converge for flat areas of $\phi^0$. One way to prevent this is to reformulate the original problem as suggested by [6, 31, 1] and detailed in the next Section. In [30], it is proposed to solve the Laplace equation in a neighborhood of the interface. The conditions on this neighborhood are chosen so as to preserve the sign of the initial level set. Nevertheless, this method cannot be extended to the entire domain $D$. In this article, however, we propose to address this shortcoming of the method by carefully designing a predictor problem which we solve to obtain an optimal initial guess $\phi^0$. Before discussing our method, we will briefly review the reformulations of [6, 1] and show why we felt the need to design a predictor problem instead.

## 3.2 Modified minimization and non-linear diffusion formulation

In order to prevent $\frac{1}{|\nabla \phi|} \to \infty$, [6] and [1] propose to use another function, $p(|\nabla \phi|)$, in

$$ \mathcal{F}(p, \phi) = \int_{\Omega} p(|\nabla \phi|) dx. \qquad (19) $$

For which the Gâteaux derivative is given by

$$ d\mathcal{F}(\phi; \psi) = \int_D \frac{p'(|\nabla \phi|)}{|\nabla \phi|} \nabla \phi \nabla \psi dx. \qquad (20) $$

Looking for the zero of this Gâteaux derivative can be interpreted as a non-linear diffusion problem

$$ \int_D \tilde{d}(|\nabla \phi|) \nabla \phi \nabla \psi dx = 0, \qquad (21) $$

where the non-linear diffusion term is given by

$$ \tilde{d}(|\nabla \phi|) := \frac{p'(|\nabla \phi|)}{|\nabla \phi|}. \qquad (22) $$
The strong form for this non-linear diffusion problem becomes

$$ \begin{cases} -\nabla \cdot (\tilde{d}(|\nabla \phi|) \nabla \phi) = 0 & \text{in } D, \\ (\tilde{d}(|\nabla \phi|) \nabla \phi) \cdot n = 0 & \text{on } \partial D. \end{cases} \qquad (23) $$

Basting et al. [6] introduced the following double well potential inspired function for $|\nabla \phi| \le 1$

$$ p(|\nabla \phi|) := \begin{cases} \frac{1}{2} (|\nabla \phi| - 1)^2 & \text{if } |\nabla \phi| > 1, \\ \frac{1}{2} |\nabla \phi|^2 (|\nabla \phi| - 1)^2 & \text{if } |\nabla \phi| \le 1. \end{cases} \qquad (24) $$

and keeping the original function for $|\nabla \phi| > 1$. The diffusion term for this function becomes

$$ \tilde{d}(|\nabla \phi|) := \begin{cases} 1 - \frac{1}{|\nabla \phi|} & \text{if } |\nabla \phi| > 1, \\ 1 - (3|\nabla \phi| - 2|\nabla \phi|^2) & \text{if } |\nabla \phi| \le 1. \end{cases} \qquad (25) $$

Adams et al. [1] pointed out that for $|\nabla \phi| < 0.5$ this diffusion term becomes positive meaning that an already flat functional profile is pushed towards $|\nabla \phi| = 0$ in any iterative scheme. In the original formulation the diffusion term is positive for $|\nabla \phi| > 1$, which causes the function profile to flatten and it is negative for $|\nabla \phi| < 1$ which causes the function profile to steepen through the iterations. Adams et al. [1] suggested the following modification of the function

$$ p(|\nabla \phi|) := \begin{cases} \frac{1}{2} \left(1 - \frac{1}{|\nabla \phi|}\right) & \text{if } |\nabla \phi| > 1, \\ \frac{(|\nabla \phi|)^3}{3} - \frac{(|\nabla \phi|)^2}{2} + \frac{1}{6} & \text{if } |\nabla \phi| \le 1. \end{cases} \qquad (26) $$

and the corresponding diffusion term becomes

$$ \tilde{d}(|\nabla \phi|) := \begin{cases} 1 - \frac{1}{|\nabla \phi|} & \text{if } |\nabla \phi| > 1, \\ 1 - (2 - |\nabla \phi|) & \text{if } |\nabla \phi| \le 1. \end{cases} \qquad (27) $$

To solve these non-linear diffusion problems, [6], [1] solve this problem by Picard iteration, which for the strong form can be written as

$$ \begin{cases} -\nabla \cdot \nabla \phi^{n+1} = -\nabla \cdot \left( \tilde{d}^* (|\nabla \phi^n|) \nabla \phi^n \right) & \text{in } D, \\ \nabla \phi^{n+1} \cdot \mathbf{n} = \left( \tilde{d}^* (|\nabla \phi^n|) \nabla \phi^n \right) \cdot \mathbf{n} & \text{on } \partial D, \end{cases} \qquad (28) $$

Here,

$$ \tilde{d}^* (|\nabla \phi|) = \frac{1}{|\nabla \phi|} \qquad (29) $$

in the original problem (15) and for $|\nabla \phi| > 1$ for (25) and (27). And for $|\nabla \phi| \le 1$, (25) yields

$$ \tilde{d}^* (|\nabla \phi|) = (3 |\nabla \phi| - 2 |\nabla \phi|^2) \qquad (30) $$

and (27) yields

$$ \tilde{d}^* (|\nabla \phi|) = (2 - |\nabla \phi|). \qquad (31) $$
## 3.3 Finite element formulation for the non-linear diffusion problem

Let $\phi^0: D \to \mathbb{R}$ be a function who decomposes our domain $D$ into two parts: a negative part $\Omega^-$ and a positive part $\Omega^+$ given by the sign of the function $\phi^0$, and an interface between $\Omega^-$ and $\Omega^+$, denoted by $\Gamma$, given by the zero contour line of $\phi^0$. This means

$$ \phi^0(x) > 0 \quad \text{for } x \in \Omega^+, \qquad (32) $$

$$ \phi^0(x) = 0 \quad \text{for } x \in \Gamma, \qquad (33) $$

$$ \phi^0(x) < 0 \quad \text{for } x \in \Omega^-. \qquad (34) $$

We do not make any further assumptions on $\phi^0$.

Let us now presume that $D$ is decomposed into mesh elements whose faces coincide with $\Gamma$. We denote this mesh by $K_h$ and the faces that coincide with $\Gamma$ by $\Gamma_h$ as well as the discretized domain by $D_h = \Omega_h^+ \cup \Omega_h^-$. In later sections in this article, we will discuss how this scheme can be extended to meshes that do not coincide with $\Gamma_h$.

Furthermore let $V_h^0$ denote the continuous piecewise linear finite element space which includes the zero Dirichlet boundary condition on $\phi$ along $\Gamma_h$, i.e.

$$ V_h^0 := \{v \in H^1(D) \cap C^0(D) \mid v|_{K_h} \in \mathbb{P}_1(K_h) \text{ for all } K_h \in \mathcal{K}_h, v|_{F_h} = 0 \text{ for all } F_h \in \Gamma_h\}. \quad (35) $$

Our finite element formulation to solve the non-linear diffusion problem (28) then reads: For $n = 0, \dots, N-1$, find $\phi_h^{n+1} \in V_h^0$ such that for all $v_h \in V_h^0$

$$ \int_{D_h} \nabla \phi_h^{n+1} \nabla v_h \, dx = \int_{D_h} (\tilde{d}^* (|\nabla \phi_h^n|) \nabla \phi_h^n) \nabla v_h \, dx. \qquad (36) $$

Note that the Neumann boundary condition from (28) is implicitly enforced through first integrating by parts and then setting the boundary term to zero. And the Dirichlet condition is strongly enforced through the finite element space. The number of iterations, $N$, is either fixed or we are iterating until

$$ \| (|\nabla \phi^{n+1}| - 1) \|_{L^2(D)} < \varepsilon \qquad (37) $$

$$ \left| \| (|\nabla \phi^{n+1}| - 1) \|_{L^2(D)} - \| (|\nabla \phi^n| - 1) \|_{L^2(D)} \right| < \varepsilon \qquad (38) $$

with a small positive real value for $\varepsilon$.

### 3.3.1 Test case for non-linear diffusion problem

Let us assume an initial function $\phi^0$ in a domain $D = [0, 1]$ defined as

$$ \phi^0(x) = \begin{cases} 0.6x - 0.25, & \text{if } x \in [0, \frac{1}{3}] \\ 0.3x - 0.15, & \text{if } x \in [\frac{1}{3}, \frac{2}{3}] \\ 0.05, & \text{if } x \in [\frac{2}{3}, 1]. \end{cases} \qquad (39) $$

displayed in Figure 1a. This initial function features three different slopes and will be used to compare the results obtained for the different potentials proposed by Kuzmin and Basting [6], Adams et al. [1], and for our prediction-correction method introduced in this article. The results shown in Figure 1 are after 10 fixed point iterations.

First, the simulation using the function (24) is conducted. In the subdomain $[\frac{1}{3}, 1]$, where the level set is initialized with a function whose gradient is less than 0.5, the solution iterates towards a function with a zero gradient norm, as shown in Figure 1b.
Secondly, the simulation using the function (26) is performed. The level set function converges to a signed distance function for any initial condition with a non-zero gradient norm, as shown in Figure 1c. However, for regions where the gradient norm is close to zero, the method does not converge to the expected signed distance function but is left unchanged. As a result, regions with a gradient norm very close to zero tend to retain their value.

Note that in these two methods, the diffusion function is changed to avoid an infinite diffusion for small gradient norms, however, the changes are not consistent with the original problem and we are not sure of the nature of the problem that we are solving any more. The idea of minimizing the distance between the gradient norm and 1 is lost and this means that they do not always perform in the ways needed.

This is why in this article, we propose another alternative. Instead of changing the non-linear diffusion problem, we are designing a predictor problem such that we obtain a good initial guess for $\phi^0$ that removes all steep or flat areas from the function to be redistanced. We then use the original unmodified method (28)-(29) to minimize the distance between $|\nabla\phi|$ and 1. The predictor-corrector scheme, proposed in this article, is able to cope with flat areas and the signed distance can be recovered for the given test case as shown in Figure 1d. We will now proceed to introducing our predictor-corrector scheme.

<figure><img src="image_2.png" /><figcaption>Figure 1: Level set after 10 fixed point iterations, the horizontal line shows the problem domain :a Piecewise-affine initialization of the level set with specified slopes: $a_1 = 0.6$, $a_2 = 0.3$ and $a_3 = 0$; b Solution with (25) ; c Solution with (27) ; d Solution with our method.</figcaption></figure>

# 4 Prediction-Correction method

In the first Subsection, the prediction-correction method is explained and detailed. The weak imposition of Dirichlet boundary conditions using the Nitsche's method is discussed, and a simple 1D example of the prediction problem is provided. In the second Subsection, the finite element discretization is explained, and the prediction-correction scheme is discretized, with specific details provided for the interface integration.

## 4.1 Predictor

Let us assume we have an initial level set function, $\phi^0$ of arbitrary shape that decomposes our domain $D$ into a positive and negative part

$$ \Omega_+ := \{x \in D : \phi^0(x) > 0\}, \quad (40) $$

$$ \Omega_- := \{x \in D : \phi^0(x) < 0\}, \quad (41) $$

$$ \Gamma := \{x \in D : \phi^0(x) = 0\}. \quad (42) $$
Furthermore, let sgn denote the sign function

$$ \operatorname{sgn}(x) = \begin{cases} 1 & \text{for } x \in \Omega_{-}, \\ -1 & \text{for } x \in \Omega_{+}, \\ 0 & \text{for } x \in \Gamma. \end{cases} \qquad (43) $$

We then formulate the following predictor problem: Find $\phi_p: D \to \mathbb{R}$, such that

$$ \left\{ \begin{array}{ll} -\Delta\phi_p = \operatorname{sgn}(x) & \text{in } D, \\ \phi_p = 0 & \text{on } \Gamma, \\ \nabla\phi_p \cdot \mathbf{n} = \operatorname{sgn}(x) & \text{on } \partial D. \end{array} \right. \qquad (44) $$

Here, we have carefully chosen sgn for the right hand side of a linear diffusion problem and as a Neumann boundary condition on $\partial D$.

The choice of sgn for the linear diffusion problem ensures the maintenance of the sign of the original function and has the same values of $|\nabla\phi_p|$ close to zero as the signed distance function, i.e. at equidistant points. This also means that there is no strong kink in the solution at $\Gamma$ as there would be if we were choosing the RHS as 1. Why this is a good first approximation for a signed distance will become even clearer in the 1D example, we give below.

For the choice of the Neumann boundary condition, let us revisit the Neumann boundary condition of the problem we wish to solve (15), i.e.

$$ \nabla\phi \cdot \mathbf{n} = \frac{\nabla\phi}{|\nabla\phi|} \cdot \mathbf{n} \qquad (45) $$

With $\mathbf{n}_\phi$ the normal to the contour lines of function $\phi$, the following holds

$$ \frac{\nabla\phi}{|\nabla\phi|} \cdot \mathbf{n} = \mathbf{n}_\phi \cdot \mathbf{n} \qquad (46) $$

Here, $\mathbf{n}_\phi$ points from $\Omega_-$ towards $\Omega_+$ and

$$ 0 \le \mathbf{n}_\phi \cdot \mathbf{n} \le 1 \text{ for } x \in \Omega_+, \qquad (47) $$

$$ -1 \le \mathbf{n}_\phi \cdot \mathbf{n} \le 0 \text{ for } x \in \Omega_-. \qquad (48) $$

We assume the "steepest" case

$$ \frac{\nabla\phi}{|\nabla\phi|} \cdot \mathbf{n} \approx \operatorname{sgn}(x) \qquad (49) $$

which yields

$$ \nabla\phi \cdot \mathbf{n} = \operatorname{sgn}(x). \qquad (50) $$

In addition to giving us a first approximation for the Neumann condition of the problem we want to solve, we also ensure that $\phi_p$ is not flattened near the boundary of $D$ (see the following 1D 4.1.1 example for more details).

### 4.1.1 Illustration of prediction problem in 1D

In one dimension the prediction problem becomes: Find a continuous function $\phi_p: [-A, A] \to \mathbb{R}$ such that

$$ -\frac{\partial^2 \phi_p(x)}{\partial x^2} = \operatorname{sgn}(x) \quad \forall x \in [-A, A], \qquad (51) $$

$$ \phi_p(x) = 0 \text{ for } x \in \{-T; T\}, \qquad (52) $$

$$ \frac{\partial \phi_p(x)}{\partial x} = -1 \text{ at } x = -A, \qquad (53) $$

$$ \frac{\partial \phi_p(x)}{\partial x} = 1 \text{ at } x = A, \qquad (54) $$
with $T \ge 0$, $T \in ]-A, A[$ and sgn defined as the following step function:

$$ \operatorname{sgn}(x) := \begin{cases} 1 & \text{if } x \in [-A, -T) \cup (T, A], \\ -1 & \text{if } x \in (-T, T), \\ 0 & \text{if } x \in \{-T, T\}. \end{cases} \qquad (55) $$

By integrating twice per interval, $[-A, -T]$, $[-T, T]$, $[T, A]$, and using the boundary conditions (52), (53) and (54), we obtain:

$$ \phi_p(x) = \begin{cases} -\frac{1}{2}x^2 - (A+1)x + \frac{1}{2}T^2 - (A+1)T & \text{if } x \in [-A, -T] \\ \frac{1}{2}x^2 - \frac{1}{2}T^2 & \text{if } x \in [-T, T] \\ -\frac{1}{2}x^2 + (A+1)x + \frac{1}{2}T^2 - (A+1)T & \text{if } x \in [T, A] \end{cases} \qquad (56) $$

Figure (2) shows a plot of this 1D solution (56). The solution has the same sign as the initial function and non-zero gradient norm except in the equidistant point $x = 0$, which is the point that is discontinuous in the gradient in the signed distance function.

Note that,

$$ \frac{\partial \phi_p}{\partial x} = -x - (A + 1) \quad \text{for } x \in [-A, -T] \qquad (57) $$

$$ \frac{\partial \phi_p}{\partial x} = x \quad \text{for } x \in [-T, T] \qquad (58) $$

$$ \frac{\partial \phi_p}{\partial x} = -x + (A + 1) \quad \text{for } x \in [T, A] \qquad (59) $$

Thus

$$ \frac{\nabla \phi_p}{|\nabla \phi_p|} = \begin{cases} -1, & x \in [-A, -T) \cup (T, A] \\ 1, & x \in (-T, T) \end{cases} \qquad (60) $$

with $\frac{\nabla \phi_p}{|\nabla \phi_p|}$ undefined at $x = 0$, exactly like the sign function. This now means in the corrector problem, for which we take (36)

$$ \int_{D_h} \nabla \phi_h \nabla v_h dx = \int_{D_h} \frac{\nabla \phi_p}{|\nabla \phi_p|} \nabla v_h dx. \qquad (61) $$

with $\nabla v_h = 1$ for piecewise linear shape functions and 1D

$$ \int_{D_h} \frac{\partial \phi_h}{\partial x} dx = \int_{D_h} \frac{\nabla \phi_p}{|\nabla \phi_p|} dx. \qquad (62) $$

which after integrating becomes

$$ \phi(x) = \frac{\nabla \phi_p}{|\nabla \phi_p|} x - T \quad \text{as} \quad \phi(-T) = 0. \qquad (63) $$

This is the signed distance function

$$ \phi(x) = \begin{cases} -x - T, & x \in [-A, 0) \\ x - T, & x \in (0, A], \end{cases} \qquad (64) $$

and hence the exact solution. This shows that by solving the predictor problem followed by one fixed point iteration, we obtain the exact solution for our simple test case. In practice, more iterations are needed, as the boundary of $\partial D$ is usually not aligned with $n_\phi$ and we introduce numerical errors for the points with $|\nabla \phi| \approx 0$. Nevertheless it motivates our choice
for this predictor problem to be a very good first initial guess for the Picard iteration scheme (36).

<figure><img src="image_3.png" /><figcaption>Figure 2: 1D solution of the prediction-correction method, with $T = 0.5$ and $A = 1$.</figcaption></figure>

## 4.2 Brief discussion on the gradient norm close to zero

As mentioned above the minimization based (re)-distancing method suffers from one major drawback, which is the division by the norm of the gradient. This norm of the gradient will inevitably become close to zero for a signed distance function in points that are equidistant to multiple points on the boundary. For example, for a circular boundary this is the point in the centre of the circle. For other geometries these might be points along equidistant lines or surfaces. The signed distance function is not differentiable at these points and piece-wise polynomial approximations around these points can yield a gradient norm close to zero. Our predictor problem also features these "smoothened" gradient norms close to zero around equidistant points. We show below that setting a very small minimal value for the gradient norm in these cases works well in the iteration scheme. Concretely, we set $\max(|\nabla\phi|, \epsilon)$ in the denominator with $\epsilon > 0$ very small. Alternatively, a narrow band approach could be used.
<figure><img src="image_4.png" /><figcaption>Figure 3: Example of pathological case during reinitialization procedure: a for 1D; for 2D</figcaption></figure>

This method allows to focus the redistancing to elements close to the boundary. In most applications, this is sufficient as the signed distance function (SDF) property is only required in a neighborhood close to the $\Gamma$. However, for any application, in which interfaces merge, equidistant points will occur even in the narrow band. This needs to be considered in any redistancing scheme.

In the next Section we detailed the finite element discretization of the proposed prediction-correction scheme.

## 4.3 Finite element discretization

In this section, we will first introduce the finite element discretization for a mesh that conforms to the interface $\Gamma$. Then, we will extend this formulation in case the mesh is not aligned with $\Gamma$.

### 4.3.1 Interface-fitted mesh

Let $D_h$ denote the discretized domain $D$ into a mesh conforming to $\Gamma_h$. Using the space $V_h^0$ previously defined by (35), the finite element problem to solve reads as follows.

**Predictor** Find $\phi_{p,h} \in V_h^0$ such that for all $v_h \in V_h^0$

$$ \int_{D_h} \nabla \phi_{p,h} \cdot \nabla v_h \, dx = \int_{D_h} \operatorname{sgn}(x) v_h \, dx + \int_{\partial D_h} \operatorname{sgn}(x) v_h \, ds. \qquad (65) $$

**Corrector** With $\phi_h^0 := \phi_{p,h}$ the initial value of the iterative scheme, the finite element element formulation to solve corrector problem becomes: For $n = 0, ..., N-1$, find $\phi_h^{n+1} \in V_h^0$ such that for all $v_h \in V_h^0$

$$ \int_{D_h} \nabla \phi_h^{n+1} \nabla v_h \, dx = \int_{D_h} \frac{\nabla \phi_h^n}{|\nabla \phi_h^n|} \nabla v_h \, dx. \qquad (66) $$
We now extend this formulation to a finite element formulation for meshes that do not conform to the interface $\Gamma$. This choice is motivated by the fact that re-distancing is often applied in the context of moving boundary problems in which re-meshing can become challenging and costly.

### 4.3.2 Unfitted finite element formulation

Let us now presume that the faces of the triangular or tetrahedral mesh $K_h$ don't coincide with $\Gamma_h$. A new finite element space is constructed, which does not incorporate the boundary conditions on the interface. Within the continuous Galerkin framework, this space is defined as

$$V_h := \{v \in H^1(D) \cap C^0(D) \mid v|_{K_h} \in [\mathbb{P}_1(K_h) \text{ for all } K_h \in \mathcal{K}_h]\}, \quad (67)$$

where $K$ denotes the elements of the mesh.

Now that the Dirichlet boundary conditions are no longer included in the finite element space, we enforce them weakly through additional terms in our finite element formulation.

**Weak enforcement of boundary conditions** Three main methods exist to enforce boundary conditions weakly: i) the penalty method, ii) Lagrange multiplier, iii) Nitsche's method. In [6], the penalty method is employed. While it is relatively straightforward to formulate, it requires careful selection of a penalty parameter suitable for the numerical problem. If the parameter is chosen too low, the solution to the problem is not unique. Although the algorithm converges well to a reinitialized level set, insufficient penalization on the interface can lead to significant interface displacement. Conversely, choosing a penalty parameter that is too high results in excessive constraints, generating matrix ill-conditioning and subsequent boundary oscillations. This phenomenon is known as locking.

In [1], the authors utilize the Lagrangian method to impose Dirichlet conditions on $\Gamma_h$. This method introduces a new unknown parameter, $\lambda$. At each iteration, the Lagrange multiplier $\lambda$ is updated by solving a PDE for the constraint. The choice of the Lagrange multiplier space is not always straightforward as it has to fulfill inf-sup stability conditions.

In this work, Nitsche's method [20] is used. Similar to the penalty method, it requires the selection of a penalty parameter. The choice of this parameter leads to difficulties akin to those previously mentioned for the penalization method. However, unlike the penalty method, Nitsche's method offers a stable and effective numerical formulation for moderate penalty parameters.

To integrate over $\Gamma_h$, which now intersects elements in an arbitrary fashion, we use the cut integration scheme introduced in the following section. 4.3.2.

**Sub-integration for interface integral** First, we reconstruct the interface in order to integrate over $\Gamma$ crossing element $K$. A simple linear interpolation of level set values is used in this work to calculate the points of intersection between the edges of the element $K$ and the interface and a piecewise linear interface, $\Gamma_K$, is reconstructed by connecting these points on the edges (see Figure (4)). Then, to integrate over the reconstructed interface part $\Gamma_K$, we generate a quadrature rule by mapping standard integration rules onto the piecewise linear approximation of $\Gamma_K$ inside the reference element. This generates a new quadrature rule with new quadrature points $\tilde{p}_i$ and new weights $\tilde{w}_i$ for each interface part. The mapping to generate the quadrature rule is mixed dimensional. An example of the approximation of the $\Gamma_K$ integral as part of $K$ is given in Figure (4).
<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of mapping corresponding to sub-integration over $\Gamma_K$.</figcaption></figure>

We can now formulate our predictor-corrector scheme for interfaces $\Gamma_h$ that freely intersect through a background mesh.

**Predictor** Using the Nitsche method, the discretized weak formulation of the predictor problem reads: Find $\phi_{p,h} \in V_h$, such that for all $v_h \in V_h$

$$a(\phi_{p,h}, v_h) = l(v_h), \qquad (68)$$

where

$$a(\phi_{p,h}, v_h) = \int_{D_h} \nabla\phi_{p,h} \cdot \nabla v_h \, dx - \int_{\Gamma_h} \nabla\phi_{p,h} \cdot \mathbf{n}_\phi v_h \, ds - \int_{\Gamma_h} \nabla v_h \cdot \mathbf{n}_\phi \phi_{p,h} \, ds \\ + \gamma_D \int_{\Gamma_h} h^{-1}\phi_{p,h} v_h \, ds \qquad (69)$$

$$l(v_h) = \int_{D_h} \operatorname{sgn}(x) \cdot v_h \, dx + \int_{\partial D_h} \operatorname{sgn}(x) v_h \, ds. \qquad (70)$$

The integration term over the discretization of $\Gamma$, denoted $\Gamma_h$, is done by the cut integration scheme introduce in the previous sub-Section 4.3.2. The rest of the terms are integrated with a standard integration scheme. A cut integration on the volume is not required as the solution to our elliptic problem is continuous at $\Gamma$ despite the discontinuous RHS term sgn($x$).

**Corrector** Using Nitsche's method, and $\phi_h^0 := \phi_{p,h}$ as the initial value of the iterative scheme, the finite element element formulation to solve corrector problem reads: For $n = 0, ..., N-1$, find $\phi_h^{n+1} \in V_h$ such that for all $v_h \in V_h$

$$a(\phi_h^{n+1}, v_h) = l(v_h, \phi_h^n), \qquad (71)$$

where

$$a(\phi_h^{n+1}, v_h) = \int_{D_h} \nabla\phi_h^{n+1} \cdot \nabla v_h \, dx - \int_{\Gamma_h} \nabla\phi_h^{n+1} \cdot \mathbf{n}_\phi v_h \, ds - \int_{\Gamma_h} \nabla v_h \cdot \mathbf{n}_\phi \phi_h^{n+1} \, ds \\ + \gamma_D \int_{\Gamma_h} h^{-1}\phi_h^{n+1} v_h \, ds \qquad (72)$$

$$l(v_h, \phi_h^n) = \int_{D_h} \frac{\nabla\phi_h^n}{\max(|\nabla\phi_h^n|, \epsilon)} \cdot \nabla v_h \, dx \qquad (73)$$

with $\epsilon$ being a very small positive scalar.
# 5 Numerical Results

In this section, we evaluate our prediction-correction method for four functions with increasing numerical complexity: (i) a smooth circular interface function, which enables a thorough assessment of the method (including convergence order, influence of the predictor problem, and analysis of the Nitsche parameter); (ii) a piecewise constant function; (iii) a complex example involving a sharply varying interface; and (iv) a three-dimensional case, demonstrating the ease with which our method can be extended to 3D problems. To analyze the convergence towards the exact solution, we denote the difference between the numerical approximation, $\phi_h$, and the exact solution, $\phi_{ex}$, as

$$e(\phi_h(x)) = \phi_h(x) - \phi_{ex}(x), \quad x \in D_h. \qquad (74)$$

We evaluate this error in the $L^2$ norm

$$\|e\|_{L^2(D_h)}^2 = \frac{\int_{D_h} e^2(\phi_h(x)) dx}{\int_{D_h} dx}. \qquad (75)$$

Furthermore, to evaluate the method's ability to converge towards the solution of the Eikonal equation (4), we introduce the Eikonal error, defined as

$$e_{\text{eik}}(\phi_h(x)) = 1 - |\nabla \phi_h(x)|, \quad x \in D_h. \qquad (76)$$

We compute this error using the $L^2$ norm

$$\|e_{\text{eik}}\|_{L^2(D_h)}^2 = \frac{\int_D e_{\text{eik}}^2(\phi_h(x)) dx}{\int_{D_h} dx}. \qquad (77)$$

Finally to assess the method's precision in preserving the initial interface $\Gamma_h^0$ we define the following norm with respect to the solution of the prediction problem $\phi_h^0$

$$\|\phi_h\|_{L^2(\Gamma)}^2 = \int_{\phi_h^0=0} (\phi_h(x) - \phi_h^0(x))^2 ds \qquad (78)$$

$$= \int_{\phi_h^0=0} (\phi_h(x))^2 ds. \qquad (79)$$

## 5.1 Classical case: circular level set

This first numerical example is used to study the order of convergence of the prediction-correction method for the three norms introduced above. We start with an initial quadratic function on the domain $D = [0, L]^2$, with $L = 1$, defined as

$$\phi_{\text{init}}(x, y) = (x - x_0)^2 + (y - y_0)^2 - r^2 \qquad (80)$$

with $(x_0, y_0)$ being the coordinate of the center of the circle and $r$ being the radius. We set $(x_0, y_0) = (0.5, 0.5)$ and $r = 0.25$. As shown in Figure (5a), this function is a parabole, and not a signed distance function. The exact solution of the reinitialized level set, denoted $\phi_{ex}$, is known as a cone centered at $(x_0, y_0)$ and such that the line of level zero coincides with that of $\phi_{init}$. The exact solution $\phi_{ex}$ is given by

$$\phi_{ex}(x, y) = \sqrt{(x - 0.5)^2 + (y - 0.5)^2} - r. \qquad (81)$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Circular level set for reinitialization with prediction-correction method with $h = L/40$ and $\gamma_D = 10^4$: a level set warped before reinitialization; b level set warped after prediction problem; c predicted level set warped after 1 iteration of correction; d level set warped after 10 iterations of correction.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: View of level set function along centre line: a initial function; b solution of predictor step; c solution of correction step (after 10 iterations).</figcaption></figure>

**Convergence study** We fix the Nitsche's parameter to $10^4$ and we study the errors defined above for different mesh sizes. By observing the evolution of the $L^2$-errors of $e$ and $e_{eik}$ as a function of iterations (see Figures (7a) and (7b)), we note that, regardless of the mesh size, these errors fall below 0.01 and 0.05, respectively, after just 3 iterations. The level set function after one iteration of the correction problem (see Figure (5c)) and after ten iterations (see (5d)) differ only very slightly. In practice, this means that only one iteration of the correction problem gives a sufficiently accurate approximation to the Eikonal equation for any application in which a strict verification of the Eikonal equation is not required, and an approximation within these error margins is acceptable (for example, in shape optimization with the level set method). The cross-section displayed in Figure (6) shows that even at the center of the domain, near the singularity, the level set shows no oscillation or locking effects.
<figure><img src="image_8.png" /><figcaption>Figure 7: Error evolution of circle example (prediction error is not plotted) with $\gamma_D = 10^4$ for different mesh sizes: a Error $\|e\|_{L^2(D)}$; b Error $\|e_{\text{eik}}\|_{L^2(D)}$.</figcaption></figure>

**Convergence order** In this section, we calculate the convergence order of the prediction-correction scheme for the three errors defined in the beginning of this section. We set $\gamma_D = 10^4$ for the convergence study. The solution is considered to have converged if

$$ \left\| \|e_{\text{eik}}(\phi^{n+1})\|_{L^2(D)} - \|e_{\text{eik}}(\phi^n)\|_{L^2(D)} \right\| < 10^{-8}. \qquad (82) $$

We study convergence with mesh refinement for the series of triangular meshes with size $h = L/60, L/80, L/120, L/160$.

The order of convergence is computed between the errors obtained with $h = L/60$ and $h = L/160$. Numerically, we observe in Figures (8a), (8b) and (8c) that the convergence order of our scheme on the interface is 1.48, the order of the Eikonal error is 0.89 and the order for the error $e$ in the whole domain is 1.39. These convergence rates are below optimal which would be order 2 for the $L^2$-errors and 1 for the Eikonal error. We have two main sources of error: (i) equidistant points with norm of the gradient close to zero (theoretically not differentiable), (ii) the splitting error in the corrector scheme from the Picard fixed point linearization. To test if the point in the centre of the circle is causing the reduction in the convergence order, we solve the circular level set problem only in a narrow band around the interface. The errors are denoted as "P.C. annulus" in Figures (8a), (8b) and (8c). We see that while the $L^2$ error does not significantly improve, we recover an optimal order of convergence for the Eikonal error.
<figure><img src="image_9.png" /><figcaption>Figure 8: The convergence order of the P.C. method with respect to the following errors: a Error $||e||_{L^2(D)}$; b Error $||e_{\text{eik}}||_{L^2(D)}$; c Error $||\phi||_{L^2(\Gamma)}$.</figcaption></figure>

**Comparison with Elliptic method:** In this section, we compare our method to the Elliptic one. The elliptic method can be used for this example because the norm of the gradient is not close to zero in large areas. Table 1 shows the errors at convergence, which remain unchanged regardless of whether the predictor problem is solved or not. This demonstrates that our method does not degrade the solution provided by conventional reinitialization approaches such as the elliptic based reinitialization.

The downside of using our method in such cases is that it requires a greater number of iterations to reach convergence compared to the classical minimization based reinitialization. In this example, the number of iterations increases by a factor of approximately 1.25.

However, the elliptic method fails to converge if the initial level set function contains large flat areas of very small values for the norm of the gradient. In these cases, solving the predictor problem becomes essential. We will demonstrate this in the next section where we take an initial flat function.
<table><caption>Table 1: Comparison between our method and the minimization method, for $N = 80$, $\gamma_D = 10000$ for minimization scheme and $\gamma_D = 20000$ for predictor-corrector scheme.</caption><thead><tr><th></th><th>Number of iteration, n</th><th>||e<sub>eik</sub>||<sub>L<sup>2</sup>(D)</sub></th><th>||e||<sub>L<sup>2</sup>(D)</sub></th><th>||e||<sub>L<sup>2</sup>(Γ)</sub></th></tr></thead><tbody><tr><td>Elliptic</td><td>279</td><td>0.01026</td><td>0.00041</td><td>0.00027</td></tr><tr><td>Predictor-Corrector</td><td>351</td><td>0.01026</td><td>0.00042</td><td>0.00027</td></tr></tbody></table>

## 5.2 Flat function

A second example is studied to understand the capability of the scheme to correct a flat function, i.e. large area with $|\nabla\phi| \approx 0$. These type of flat functions are a challenge for the minimization based reinitialization scheme which includes a $\frac{1}{|\nabla\phi|}$ term.

Let the initial level set function be a step function (see Figure (9)) on the domain $D = [0, 1]^2$ defined by:

$$ \phi_{\text{init}}(x, y) = \begin{cases} -1 & \text{if } x > 0.5, \\ 0 & \text{if } x = 0.5, \\ 1 & \text{if } x < 0.5. \end{cases} \qquad (83) $$

<figure><img src="image_10.png" /><figcaption>Figure 9: Step-shape warped for reinitialization step with P.C. method, $\gamma_D = 10$ and $h = L/40$: a Level set function before reinitialization; b level set function after prediction problem; c level set function after P.C. reinitialization.</figcaption></figure>

Without a predictor problem, the initial minimization problems right hand side becomes infinity everywhere in the Picard fixed point iteration scheme. The solution of the predictor problem is a level set with the same sign as the initial level set, while correcting the zero gradient norms so that the corrector problem is well-defined. As observed in Figures (9) and (10) for the level set solution of the predictor problem, no element $K$ of the mesh satisfies: $|\nabla\phi(K)| = 0$. From this predictor problem, only a few iterations are necessary to obtain a reinitialized level set function.

To complete our study, we initialized a variant of the previous level set to assess the ability of the scheme to correct very irregular functions, containing areas with zero gradient norms and areas with very high gradient norms. The initialized level set is defined as:

$$ \phi_{\text{init}}(x, y) = \arctan(10(x - 0.5)) + \xi(x) \qquad (84) $$
<figure><img src="image_11.png" /><figcaption>Figure 10: Step-shape level set contours for reinitialization step with P.C. method, $\gamma_D = 10$ and $h = L/40$: a level set contours before reinitialization; b level set contours after predictor problem; c level set contours after P.C. reinitialization.</figcaption></figure>

with $\xi$ a function of Gaussian small perturbation generated randomly.

This example highlights the relevance of the predictor problem for highly irregular initial level set functions. Indeed, as shown in Figures (13a) and (13b), the errors obtained using the prediction-correction scheme are significantly lower than those resulting from solving the minimization problem alone. Looking at the level set function at convergence displayed in Figures (11) and (12), we observe that the minimization problem alone is unable to smooth out significant irregularities in the vicinity of the domain boundaries while the predictor corrector problem yields a smooth approximation of the signed distance. Furthermore, it is observed that the elliptic scheme requires more iterations to converge to the solution. Finally, the error on $e$ decreases over iterations in the prediction-correction method, whereas this is not the case for the elliptic method as shown in Figure (13b).

<figure><img src="image_12.png" /><figcaption>Figure 11: Warped level set for arctan-problem with $\gamma_D = 10$ and $h = L/41$: a level set function before reinitialization; b level set function after minimization based reinitialization; c level set function after P.C. reinitialization.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Perturbed arctan level set contours with $\gamma_D = 10$ and $h = L/41$: a level set contours before reinitialization; b level set contours after minimization based reinitialization; c level set contours after P.C. reinitialization.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Error evolution for $\gamma_D = 10$ and $h = L/41$: a error $\|e_{eik}\|_{L^2(D)}$; b error $\|e\|_{L^2(D)}$.</figcaption></figure>

To illustrate the ability of our method to reinitialize level set with a complex interface, we'll look at an example next that features an irregular boundary.

## 5.3 Complex interface: star level set

In this section, we test our scheme on a star-shaped level set function (see Figures (14a) and (15a)) defined on the domain $D = [0, 1]^2$ given by

$$ \phi_{\text{init}}(x, y) = \frac{1}{2} \left( \sqrt{(x - x_0) + (y - y_0)} - \frac{1}{10} \cos \left( n \arctan \left( \frac{y - y_0}{x - x_0} \right) + 2 \right) + \frac{1}{4} \right) \quad (85) $$

with $(x_0, y_0)$ the center of the star and $n = 8$ the number of rays. The distance function associated with this star-shaped form is particularly complex to compute. The interface exhibits large areas with high curvature and a lot of equidistant points along lines between the star rays for which the norm of the gradient is zero. The preservation of the geometry as well as the large number of equidistant lines makes this a challenging test problem.
<figure><img src="image_15.png" /><figcaption>Figure 14: Warped star level set function with $\gamma_D = 1000$ and $h = L/320$:a before reinitialization; b after predictor; c after corrector reinitialization.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Star level set contours $\gamma_D = 1000$ and $h = L/320$:a before reinitialization; b predictor problem; c after corrector reinitialization.</figcaption></figure>

## 5.4 Simulation in three dimensions

As our reinitialization scheme is purely based on PDEs and a cut integration of the interface local to each cell, our method can be parallelized easily and can hence be used for complex 3D geometries. In this section, we consider a doughnut-shaped level set function defined on the domain $D = [-L/2, L/2]^3$, with $L = 2.5$ defined as follows:

$$ \phi_{\text{init}}(x, y, z) = \left(R - \sqrt{x^2 - y^2}\right) + z^2 - r^2 \qquad (86) $$

with, $R = 1.8$ and $r = 0.55$. We use a tetrahedral mesh with size $h = L/101$, resulting in over two million unknowns. Figure 16 shows the initial level-set which is not a signed distance function as can be seen by the densely packed surfaces close to the zero contour surface. After solving the predictor problem followed by 10 iterations of the corrector problem, Figure 17 illustrates that the method gives a very good approximation of the signed distance function associated with the initial interface. The simulation was completed in 7 minutes and 54 seconds using 6 processors based on FEniCSx [4], CutFEMx [10] as well as PETSc [3].
<figure><img src="image_17.png" /><figcaption>Figure 16: Doughnut level set contours before prediction-correction reinitialization: a view of a 2D cross-section of the level lines with the interface in 3D dimensions; b view of a cross section of the level set.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Doughnut level set contours after prediction-correction reinitialization: a view of a 2D cross-section of the level lines with the interface in 3D dimensions; b view of a cross section of the level set.</figcaption></figure>

# 6 Conclusions

In this article, we have presented a new finite element-based technique to compute the signed distance function starting from a level set function of general shape. Our method is based on a predictor-corrector formulation, where the predictor step is a linear diffusion problem with a discontinuous right hand side and the corrector step is a nonlinear diffusion problem based on the least squares minimization of the Eikonal equation. We presented both a finite element formulation for a mesh that fits to the zero contour line of the initial level set function as well as a cut finite element formulation that allows for a mesh that does not fit to the zero contour line.

We demonstrated that our predictor-corrector method is robust and converges to a signed distance function for a variety of initial level set functions, including smooth, discontinuous, and irregular functions. The method is also able to handle complex geometries in three dimensions and is easy to parallelize. We showed that the method is able to reinitialize the level set function with very few iterations, even for large meshes with millions of unknowns.
# Acknowledgements

This research was supported by the French Aerospace Lab ONERA under internal research funding. The authors gratefully acknowledge this support.

# References

[1] Thomas Adams, Stefano Giani, and William M. Coombs. A high-order elliptic pde based level set reinitialisation method using a discontinuous galerkin discretisation. *Journal of Computational Physics*, 379:373–391, 2019. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2018.12.003. URL https://www.sciencedirect.com/science/article/pii/S0021999118307915.

[2] Grégoire Allaire, François Jouve, and Anca-Maria Toader. A level-set method for shape optimization. *Comptes Rendus. Mathématique*, 334(12):1125–1130, 2002. doi: 10.1016/S1631-073X(02)02412-3. URL https://www.numdam.org/articles/10.1016/S1631-073X(02)02412-3/.

[3] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, Stefano Zampini, and Hong Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.6, Argonne National Laboratory, 2015. URL http://www.mcs.anl.gov/petsc.

[4] I. A. Baratta, J. P. Dean, J. S. Dokken, M. Habera, C. N. Richardson, J. S. Hale, M. E. Rognes, M. W. Scroggs, N. Sime, and G. N. Wells. Dolfinx: The next generation fensics problem solving environments. *Zenodo*, 2023. doi: https://doi.org/10.5281/zenodo.10447666.

[5] Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth N. Wells. Dolfinx: The next generation fensics problem solving environment, December 2023. URL https://doi.org/10.5281/zenodo.10447666.

[6] Christopher Basting and Dmitri Kuzmin. A minimization-based finite element formulation for interface-preserving level set reinitialization. *Computing*, 95:13–25, 2013. URL https://api.semanticscholar.org/CorpusID:17007867.

[7] Erik Burman, Susanne Claus, Peter Hansbo, Mats Larson, and Andre Massing. Cutfem: Discretizing geometry and partial differential equations. *International Journal for Numerical Methods in Engineering*, 104, 12 2014. doi: 10.1002/nme.4823.

[8] Elisabetta Carlini, Maurizio Falcone, and Philippe Hoch. A generalized fast marching method on unstructured triangular meshes. *SIAM Journal on Numerical Analysis*, 51, 01 2013. doi: 10.1137/110833610.

[9] David L. Chopp. Computing minimal surfaces via level set curvature flow. *Journal of Computational Physics*, 106(1):77–91, 1993. ISSN 0021-9991. doi: https://doi.org/10.1006/jcph.1993.1092. URL https://www.sciencedirect.com/science/article/pii/S0021999183710922.

[10] Susanne Claus. CutFEMx: A library for unfitted finite element methods on complex geometries. https://github.com/sclaus2/cutfemx, 2022. Accessed: May 15, 2025.
[11] Susanne Claus and Pierre Kerfriden. A cutfem method for two-phase flow problems. *Computer Methods in Applied Mechanics and Engineering*, 348:185-206, 2019.

[12] Susanne Claus, Pierre Kerfriden, Faezeh Moshfeghifar, Sune Darkner, Kenny Erleben, and Christian Wong. Contact modeling from images using cut finite element solvers. *Advanced Modeling and Simulation in Engineering Sciences*, 8(1), December 2021. doi: 10.1186/s40323-021-00197-2. URL https://hal.science/hal-03507472.

[13] Charles Dapogny and Florian Feppon. Shape optimization using a level set based mesh evolution method: an overview and tutorial. *Comptes Rendus. Mathématique*, 361:1267-1332, 2023. doi: 10.5802/crmath.498.

[14] Kenny Erleben and Henrik Dohlmann. *Signed Distance Fields Using Single-Pass GPU Scan Conversion of Tetrahedra*, pages 741–763. Addison-Wesley, 2008. ISBN 9780321515261.

[15] Zhisong Fu, Robert Kirby, and Ross Whitaker. A fast iterative method for solving the eikonal equation on tetrahedral domains. *SIAM journal on scientific computing : a publication of the Society for Industrial and Applied Mathematics*, 35:c473-c494, 01 2013. doi: 10.1137/120881956.

[16] Felix Henri. *Améliorations des méthodes Level Set pour l'impact de goutte de pluie*. Theses, Université de Bordeaux, December 2021. URL https://theses.hal.science/tel-03527520.

[17] R. Kimmel and J. A. Sethian. Computing geodesic paths on manifolds. *Proceedings of the National Academy of Sciences*, 95(15):8431-8435, 1998-07-21. doi: 10.1073/pnas.95.15.8431. URL https://www.pnas.org/doi/10.1073/pnas.95.15.8431. Publisher: Proceedings of the National Academy of Sciences.

[18] Chunming Li, Chenyang Xu, Changfeng Gui, and M.D. Fox. Distance regularized level set evolution and its application to image segmentation. *Image Processing, IEEE Transactions on*, 19:3243–3254, 01 2011. doi: 10.1109/TIP.2010.2069690.

[19] Neil Molino, Robert Bridson, Joseph Teran, and Ronald Fedkiw. A crystalline, red green strategy for meshing highly deformable objects with tetrahedra. *12th Int. Meshing Roundtable*, pages 103–114, 01 2003.

[20] Johannes C. C. Nitsche. Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. *Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg*, 36:9–15, 1971. URL https://api.semanticscholar.org/CorpusID:122609590.

[21] Stanley Osher and Ronald Fedkiw. *Level set methods and dynamic implicit surfaces*. Springer Science & Business Media, 2003.

[22] Jianliang Qian, Yong-Tao Zhang, and Hong-Kai Zhao. Fast sweeping methods for eikonal equations on triangular meshes. *SIAM Journal on Numerical Analysis*, 45(1):83–107, 2007-01. ISSN 0036-1429. doi: 10.1137/050627083. URL https://epubs.siam.org/doi/abs/10.1137/050627083. Publisher: Society for Industrial and Applied Mathematics.

[23] J A Sethian. A fast marching level set method for monotonically advancing fronts. *Proceedings of the National Academy of Sciences of the United States of America*, 93(4):
1591-1595, 1996-02-20. ISSN 0027-8424. URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC39986/.

[24] J.A. Sethian. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 1999. ISBN 9780521645577. URL https://books.google.fr/books?id=ErpOoynE4dIC.

[25] James A Sethian. A fast marching level set method for monotonically advancing fronts. *proceedings of the National Academy of Sciences*, 93(4):1591-1595, 1996.

[26] Modesar Shakoor, Benjamin Scholtes, Pierre-Olivier Bouchard, and Marc Bernacki. An efficient and parallel level set reinitialization method - Application to micromechanics and microstructural evolutions. *Applied Mathematical Modelling*, 39(23-24):7291-7302, December 2015. doi: 10.1016/j.apm.2015.03.014. URL https://minesparis-psl.hal.science/hal-01139858.

[27] M Sussman, P Smereka, and S Osher. A level set approach for computing solutions to incompressible two-phase flow. *Journal of Computational Physics*; (United States), 9 1994. ISSN 0021-9991. doi: 10.1006/jcph.1994.1155. URL https://www.osti.gov/biblio/7075774.

[28] Richard Tsai and Stanley Osher. Level set methods and their applications in image science. *Communications in mathematical sciences*, 1, 12 2003. doi: 10.4310/CMS.2003.v1.n4.a1.

[29] Yen-Hsi Richard Tsai, Li-Tien Cheng, Stanley Osher, and Hong-Kai Zhao. Fast sweeping algorithms for a class of hamilton-jacobi equations. *SIAM Journal on Numerical Analysis*, 41(2):673-694, 2003. doi: 10.1137/S0036142901396533. URL https://doi.org/10.1137/S0036142901396533.

[30] Tianju Xue, WaiChing Sun, Sigrid Adriaenssens, Yujie Wei, and Chuanqi Liu. A new finite element level set reinitialization method based on the shifted boundary method. *Journal of Computational Physics*, 438:110360, 2021. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2021.110360. URL https://www.sciencedirect.com/science/article/pii/S0021999121002552.

[31] Tianju Xue, WaiChing Sun, Sigrid Adriaenssens, Yujie Wei, and Chuanqi Liu. A new finite element level set reinitialization method based on the shifted boundary method. *Journal of Computational Physics*, 438:110360, 2021. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2021.110360. URL https://www.sciencedirect.com/science/article/pii/S0021999121002552.

[32] Hongkai Zhao. A fast sweeping method for eikonal equations. *Mathematics of Computation*, 74(250):603-627, 2004-05-21. ISSN 0025-5718, 1088-6842. doi: 10.1090/S0025-5718-04-01678-3. URL https://www.ams.org/mcom/2005-74-250/S0025-5718-04-01678-3/.