# A partition-of-unity-based finite element method for level sets

Stéphane Valance, René de Borst, Julien Réthoré, Michel Coret

► To cite this version:

Stéphane Valance, René de Borst, Julien Réthoré, Michel Coret. A partition-of-unity-based finite element method for level sets. International Journal for Numerical Methods in Engineering, 2008, 76 (10), pp.1513-1527. 10.1002/nme.2371. hal-01004925

HAL Id: hal-01004925

https://hal.science/hal-01004925

Submitted on 11 Mar 2018

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 partition-of-unity-based finite element method for level sets

Stéphane Valance¹, René de Borst²,*,†, Julien Réthoré¹ and Michel Coret¹

¹LaMCoS, CNRS UMR 5259, INSA-Lyon, Villeurbanne, France

²Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands

Level set methods have recently gained much popularity to capture discontinuities, including their possible propagation. Typically, the partial differential equations that arise in level set methods, in particular the Hamilton-Jacobi equation, are solved by finite difference methods. However, finite difference methods are less suited for irregular domains. Moreover, it seems slightly awkward to use finite differences for the capturing of a discontinuity, while in a subsequent stress analysis finite elements are normally used. For this reason, we here present a finite element approach to solving the governing equations of level set methods. After a review of the governing equations, the initialization of the level sets, the discretization on a finite domain, and the stabilization of the resulting finite element method will be discussed. Special attention will be given to the proper treatment of the internal boundary condition, which is achieved by exploiting the partition-of-unity property of finite element shape functions. Finally, a quantitative analysis including accuracy analysis is given for a one-dimensional example and a qualitative example is given for a two-dimensional case with a curved discontinuity.

KEYWORDS: level sets; finite elements; partition of unity; evolving discontinuities

## 1. INTRODUCTION

In the late 1980s, Osher and Sethian [1] have suggested an elegant method to numerically model hypersurfaces. The starting point is the definition of a scalar level set function $\phi$. The zero-isolevel contour of this function describes the hypersurface, while the signed distance provided by the level set function enables the simulation of the evolution of the hypersurface.

Initially, level set methods were applied to the computation of phase changes in flows as driven by a diffusion equation. Subsequent applications have also included weather predictions and image analysis [2]. More recently, they have also been used in conjunction with finite element methods that exploit the partition-of-unity property of finite element shape functions to capture crack

<footnote>*Correspondence to: René de Borst, Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, NL-5600 MB Eindhoven, The Netherlands.</footnote>

<footnote>†E-mail: R.d.Borst@tue.nl</footnote>
propagation, especially in three-dimensional cases [3, 4], to model holes and inclusions [5, 6], or to model biofilm growth [7].

Originally, the evolution equations that arise in a level set method have been solved using finite difference methods. However, difference methods are less suited for irregular domains, and it seems less elegant and even slightly awkward to use finite differences to capture a discontinuity, while in a subsequently stress analysis finite elements are employed, e.g. those that exploit the partition-of-unity property of finite element shape functions. Moreover, for the stress analysis of bodies with propagating cracks or other evolving discontinuities the construction of the enriched functions that are utilized in partition-of-unity-based finite element functions is intimately related to the geometry of the propagating discontinuity. For these reasons, the integration of the level set method that describes the discontinuity in a finite element method, which also analyses the effects of the evolving discontinuity, may have advantages.

Finite element schemes for solving the equations that describe the level set evolution are encountered less frequently in the literature, see [8–10] for exceptions. These contributions give a solid framework, but can be improved further with respect to the initialization of the zero isolevel and the treatment of the internal boundary condition that arises from the very condition that the level set function must vanish at the propagating discontinuity. The issue how to impose essential boundary conditions in enriched finite element methods that exploit the partition-of-unity property of the standard polynomial shape functions was addressed in detail in [11], but the procedures described therein do not seem to be directly applicable to an evolving internal boundary.

This paper begins with a concise review of the governing equations for level set methods, including the initialization of the level sets, and the discretization on a finite domain using a finite element method. To properly capture the internal essential boundary condition, the finite element method is enhanced by exploiting the partition-of-unity property of finite element shape functions. Since the resulting equations have a non-symmetric character, stability and uniqueness are not ensured. For this reason, a stabilization term is added using a Galerkin least-squares formalism [12]. The ensuing algorithm is presented and the section is concluded by a description how to initialize the level set method from measured (discrete) data. Finally, an accuracy analysis is given for a one-dimensional example and a qualitative example is provided for a two-dimensional case with a curved discontinuity.

## 2. LEVEL SET METHODS

In level set methods a hypersurface $\phi = \phi(\mathbf{x}, t)$ is defined on the domain of interest $\Omega$. The basic idea is that the intersection of this hypersurface with the zero level, i.e. $\phi=0$, defines the internal discontinuity. From a Lagrangian point of view, stationarity of the level set field requires that

$$ \frac{D\phi}{Dt} = 0 \qquad (1) $$

or, equivalently,

$$ \frac{\partial \phi}{\partial t} + \nabla \phi \cdot \frac{\partial \mathbf{x}}{\partial t} = 0 \qquad (2) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Domains around the internal discontinuity $\Gamma$ with $\Omega^x = \Omega_+^x \cup \Omega_-^x$.</figcaption></figure>

Equation (2) defines the propagation of the zero-isolevel contour or, equivalently, of the evolving discontinuity $\Gamma$. We next define the velocity $V_n$ that is normal to the discontinuity $\Gamma$, Figure 1:

$$ \frac{\partial \mathbf{x}}{\partial t} \cdot \mathbf{n} = V_n \qquad (3) $$

for all points $\mathbf{x}$ at the zero-isolevel contour, and $\mathbf{n}$ the normal vector at $\mathbf{x}$:

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

With these definitions, and adding the initial condition $\phi(\mathbf{x}, t=0) = f(\mathbf{x})$, with $f(\mathbf{x})$ known, the initial value problem that describes the propagation of the zero-isolevel contour becomes [1]

$$ \forall \mathbf{x} \in \Omega, \quad \left\{ \begin{array}{l} \frac{\partial \phi}{\partial t} + V_n \|\nabla \phi\| = 0 \\ \phi(\mathbf{x}, t=0) = f(\mathbf{x}) \end{array} \right. \qquad (5) $$

Mourad et al. [10] assume that $\|\nabla\phi\|=1$ is maintained everywhere in $\Omega$ for $t>0$, so that the formulation of the initial value problem (5) simplifies to

$$ \forall \mathbf{x} \in \Omega, \quad \left\{ \begin{array}{l} \frac{\partial \phi}{\partial t} + V_n = 0 \\ \phi(\mathbf{x}(t), t=0) = f(\mathbf{x}) \end{array} \right. \qquad (6) $$

A difficulty resides in the fact that the velocity $V_n$ is defined only at the internal discontinuity $\Gamma$. A standard methodology to circumvent this difficulty is to extend the velocity field $V_n$ at the discontinuity $\Gamma$ such that

$$ \begin{array}{ll} \nabla V_e \cdot \nabla \phi = 0, & \mathbf{x} \in \Omega/\Gamma \\ V_e = V_n, & \mathbf{x} \in \Gamma \end{array} \qquad (7) $$
where $V_e$ is the extended velocity field [10, 13]. The orthogonality requirement implied in Equation (7)$_1$ has some interesting consequences. As shown by Zhao et al. [13], cf. [10], it preserves the initial signed-distance function for sufficiently smooth $\phi$ and $V_e$. Keeping the signed-distance function smoothens the zero-isolevel contour around the discontinuity and enables a more regular motion of the front [14].

### 3. WEAK FORMULATION AND FINITE ELEMENT IMPLEMENTATION

#### 3.1. Propagation equation

To obtain a finite element discretization, we multiply Equation (6)$_1$ by a test function $\psi$, integrate over $\Omega$, and use Equation (7)$_2$ to obtain

$$ \int_{\Omega} \psi \left( \frac{\partial \phi}{\partial t} + V_e \right) dV = 0 \qquad (8) $$

With a standard interpolation

$$ \phi = N\phi \qquad (9) $$

$N$ being the array that contains the interpolation polynomials and $\phi$ the array that contains the nodal values of $\phi$; using a Bubnov-Galerkin method in the sense that the test function $\psi$ is also interpolated in the sense of Equation (9), Equation (8) becomes

$$ M \frac{\partial \phi}{\partial t} + \int_{\Omega} N^T V_e dV = 0 \qquad (10) $$

with the mass matrix $M$ defined as

$$ M = \int_{\Omega} N^T N dV \qquad (11) $$

Using a first-order Euler time discretization scheme

$$ \frac{\partial \phi}{\partial t} \approx \frac{\phi_{t+\Delta t} - \phi_t}{\Delta t} \qquad (12) $$

we obtain the following discrete evolution equation:

$$ \phi_{t+\Delta t} = \phi_t - \Delta t M^{-1} \int_{\Omega} N^T V_e dV \qquad (13) $$

#### 3.2. Enforcement of constant gradient norm and zero-isolevel contour

The preceding time-marching scheme holds subject to the requirement that $||\nabla\phi||=1$, which ensures that the signed-distance function remains valid, and to the constraint that at the evolving discontinuity the level set function vanishes. Accordingly, the following supplementary conditions must be imposed locally:

$$ \begin{aligned} \|\nabla\phi\| = 1 &\iff \nabla\phi \cdot \nabla\phi = 1 & \forall \mathbf{x} \in \Omega/\Gamma \\ \phi = 0 &\iff \mathbf{x} \in \Gamma \end{aligned} \qquad (14) $$
Since the first equation of set (14) is non-linear, an iterative procedure at each time step is necessary for its solution. Adopting a first-order Taylor series at $\nabla\phi_t$ one obtains

$$g(\nabla\phi_t + \varepsilon_\phi) \approx g(\nabla\phi_t) + \frac{\partial g}{\partial\nabla\phi}(\nabla\phi_t)\varepsilon_\phi \quad (15)$$

where

$$g(\nabla\phi_t) = \nabla\phi_t \cdot \nabla\phi_t - 1 \quad (16)$$

Considering that $\varepsilon_\phi \ll 1$, a Newton–Raphson procedure can now be used. At iteration $k+1$ a field $\phi_t^{k+1}$ is sought such that

$$\begin{aligned} 2\nabla\phi_t^k \cdot \nabla\phi_t^{k+1} &= 1 + \nabla\phi_t^k \cdot \nabla\phi_t^k & \forall \mathbf{x} \in \Omega/\Gamma \\ \phi_t^{k+1} &= 0 & \forall \mathbf{x} \in \Gamma \end{aligned} \quad (17)$$

Subsequently, the problem is cast into a weak format. With $\mathcal{W}$ the space of sufficiently regular fields, e.g. those that belong to $H^1(\Omega)$, the weak formulation of Equation (17) becomes

$$\forall \psi \in \mathcal{W} = \{\phi \in \mathcal{W} | \phi = 0 \text{ on } \Gamma\} \quad (18)$$

Find $\phi_t^{k+1}$ such that

$$\mathcal{W}^k(\psi, \phi_t^{k+1}) = \ell^k(\psi) \quad (19)$$

where

$$\mathcal{W}^k(\psi, \phi_t^{k+1}) = 2 \int_{\Omega} \psi \nabla \phi_t^k \cdot \nabla \phi_t^{k+1} \, dV \quad (20)$$

$$\ell^k(\psi) = \int_{\Omega} \psi(1 + \nabla\phi_t^k \cdot \nabla\phi_t^k) \, dV \quad (21)$$

A major difficulty regarding the solution to Equation (19) arises from the essential condition that is imposed at the internal boundary $\Gamma$, cf. Equation (14)$_2$. General approaches to impose boundary conditions, or linear relationships between degrees-of-freedom, in finite element methods are to reduce the projection space either by suppressing degrees of freedom or by using Lagrange multipliers, which effectively enlarges the projection space. Recently, Moës et al. [11] have discussed solutions within the context of partition-of-unity-based finite element methods, but these solutions do not seem readily applicable to the case of an evolving internal boundary.

Herein, we impose the internal boundary condition by exploiting the partition-of-unity property of finite element shape functions [15–17]. In this method, the projection space is enriched by new functions, which are multiplied by the classical polynomial functions on the support of a node:

$$\phi(\mathbf{x}) = \sum_{i=1}^{N} N_i(\mathbf{x}) \left( \bar{\phi}_i + \sum_{j=1}^{m} \tilde{N}_j(\mathbf{x}) \tilde{\phi}_{ij} \right) \quad (22)$$

where $N_i(\mathbf{x})$ are the traditional polynomial shape functions, $\tilde{N}_j(\mathbf{x})$ the enrichment functions at node $i$, $\bar{\phi}_i$ and $\tilde{\phi}_i$ the degrees of freedom that relate to the standard and the enhanced interpolations, respectively, $N$ the number of nodes per element, and $m$ the number of enhanced functions per
node. In the present case the essential boundary condition is imposed at the evolving internal boundary $\Gamma$, and we have $m=1$ and

$$ \tilde{N}_1(\mathbf{x}) = \phi_t^0(\mathbf{x}) + \delta^0 \qquad (23) $$

which represents the value of the level set function $\phi$ at the initialization of the iterative procedure, i.e. the solution obtained at the end of the previous propagation step, plus a possibly non-zero value $\delta^0$ of $\phi(\mathbf{x})$ that is imposed at $\Gamma$. With the usual interpolation, we can rewrite this enrichment as

$$ \tilde{N}_1(\mathbf{x}) = \sum_{k=1}^{\mathcal{N}} N_k(\mathbf{x}) \phi_t^0(\mathbf{x}_k) + \delta^0 \qquad (24) $$

and substitution into Equation (22) subsequently yields

$$ \phi(\mathbf{x}) = \sum_{i=1}^{\mathcal{N}} N_i(\mathbf{x}) \left( \bar{\phi}_i + \sum_{k=1}^{\mathcal{N}} N_k(\mathbf{x}) \phi_t^0(\mathbf{x}_k) + \delta^0 \right) \tilde{\phi}_i \qquad (25) $$

We now remove the standard part from the interpolation, i.e. from which a value at $\Gamma$ can stem that is unequal to $\phi_t^0(\mathbf{x})+\delta^0$, and obtain

$$ \phi(\mathbf{x}) = \sum_{i=1}^{\mathcal{N}} N_i(\mathbf{x}) \left( \sum_{k=1}^{\mathcal{N}} N_k(\mathbf{x}) \phi_t^0(\mathbf{x}_k) + \delta^0 \right) \tilde{\phi}_i \qquad (26) $$

The first term enforces a vanishing value of $\phi(\mathbf{x})$ at the zero-isolevel contour, while the second enrichment sets the value of $\phi(\mathbf{x})$ at the zero-isolevel contour equal to $\delta^0$. Accordingly, this technique enables the imposition of Dirichlet conditions at an internal, possibly evolving, boundary. For the present purpose of enforcing the zero-isolevel contour, it suffices to include the first term.

Assembling the standard and enriched interpolation functions in a matrix $N_e$, and assembling the degrees of freedom needed for the standard interpolation and those that relate to the enrichment in a single array $\hat{\phi}$, we can express Equation (26) in a more compact manner, e.g. [17]

$$ \phi = N_e \hat{\phi} \qquad (27) $$

where, for notational simplicity, the explicit dependence of $\phi$ and $N_e$ on $\mathbf{x}$ has been dropped. The enrichment is only used in the area of interest, i.e. at the nodes that belong to $\Omega^0$, Figure 1. A well-known issue of partition-of-unity-based finite element methods is that the enrichment should be regularized outside the area of interest in order to avoid perturbations [16, 18]. For the present purpose, we use only the enriched interpolation functions, and this issue does not arise.

Requiring that the set (19)–(21) holds for any variationally admissible $\psi$, and using the discretization (27), Equation (19) can be expressed in a standard matrix-vector format as follows:

$$ A^k \hat{\phi}_i^{k+1} = b^k \qquad (28) $$

where

$$ A^k = 2 \int_{\Omega} N_e^T (\nabla N_e \hat{\phi}_i^k)^T \nabla N_e dV \qquad (29) $$

$$ b^k = \int_{\Omega} (1 + (\nabla N_e \hat{\phi}_i^k)^T (\nabla N_e \hat{\phi}_i^k)) N_e^T dV \qquad (30) $$
### 3.3. Stabilization

The weak form (19) is not completely satisfactory owing to its non-symmetric nature. As a consequence, uniqueness and stability of the solution are not ensured. To overcome this, a Galerkin least-squares stabilization term is added to the bilinear form $\mathcal{A}^k$ [12], so that it transforms into

$$ \forall \psi \in \mathcal{U} = \{\phi \in \mathcal{U} | \phi = 0 \text{ on } \Gamma\} \qquad (31) $$

Find $\phi_t^{k+1}$ such that

$$ \mathcal{A}^k(\psi, \phi_t^{k+1}) + \mathcal{A}_{\text{GLS}}^k(\psi, \phi_t^{k+1}) = \ell^k(\psi) + \ell_{\text{GLS}}^k(\psi) \qquad (32) $$

where

$$ \mathcal{A}_{\text{GLS}}^{k}(\psi, \phi_{t}^{k+1}) = \sum_{E} 4\tau^{e} \int_{\Omega_{e}} (\nabla\psi \cdot \nabla\phi_{t}^{k})(\nabla\phi_{t}^{k} \cdot \nabla\phi_{t}^{k+1}) dV \qquad (33) $$

$$ \ell_{\text{GLS}}(\psi) = \sum_{E} 2\tau^{e} \int_{\Omega_{e}} (\nabla\psi \cdot \nabla\phi_{t}^{k})(1 + \nabla\phi_{t}^{k} \cdot \nabla\phi_{t}^{k}) dV \qquad (34) $$

and the stabilization parameter

$$ \tau^e = \frac{h_e}{2 \| \nabla \phi_t^k \|} \qquad (35) $$

In Equations (32)–(35) E denotes the set of elements of the mesh, and $h_e$ is the element size. Requiring that the set (32)–(35) holds for any variationally admissible $\psi$, and using the discretization (27), Equation (32) can be expressed in a standard matrix-vector format as follows:

$$ (\mathbf{A}^k + \mathbf{A}_{\text{GLS}}^k) \hat{\mathbf{\Phi}}_t^{k+1} = \mathbf{b}^k + \mathbf{b}_{\text{GLS}}^k \qquad (36) $$

with $\mathbf{A}^k$ and $\mathbf{b}^k$ defined in Equations (21), and

$$ \mathbf{A}_{\text{GLS}}^{k} = \sum_{E} 4\tau^{e} \int_{\Omega_{e}} \nabla \mathbf{N}_{e}^{\text{T}} (\nabla \mathbf{N}_{e} \hat{\mathbf{\Phi}}_{t}^{k})^{\text{T}} (\nabla \mathbf{N}_{e} \hat{\mathbf{\Phi}}_{t}^{k}) \nabla \mathbf{N}_{e} dV \qquad (37) $$

$$ \mathbf{b}_{\text{GLS}}^{k} = \sum_{E} 2\tau^{e} \int_{\Omega_{e}} (1 + (\nabla \mathbf{N}_{e} \hat{\mathbf{\Phi}}_{t}^{k})^{\text{T}} (\nabla \mathbf{N}_{e} \hat{\mathbf{\Phi}}_{t}^{k})) \nabla \mathbf{N}_{e}^{\text{T}} (\nabla \mathbf{N}_{e} \hat{\mathbf{\Phi}}_{t}^{k}) dV \qquad (38) $$

The set (36) can now be solved in a standard Newton-Raphson manner by repeatedly solving the set of equations

$$ \hat{\mathbf{\Phi}}_{t}^{k+1} = (\mathbf{A}^{k} + \mathbf{A}_{\text{GLS}}^{k})^{-1} (\mathbf{b}^{k} + \mathbf{b}_{\text{GLS}}^{k}) \qquad (39) $$

The criterion used for convergence is $||\nabla\phi|| - 1 < \varepsilon$, with $\varepsilon$ a sufficiently small number. The method is applied only over a small area around the front, i.e. $\Omega^0 \cup \Omega^1$, Figure 1.

### 3.4. The complete algorithm

As the preceding re-initialization step is achieved using an enrichment, a slight modification must be applied to the propagation scheme. First, we project the converged solution at the re-initialization step, $\hat{\phi}_t^n$, onto the classical Galerkin space:

$$ \forall \mathbf{x} \in \Omega, \quad \min(\hat{\phi}_t - \phi_t) \qquad (40) $$
In a discrete format we then have

$$ \phi_t = \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} \widehat{\phi}_t^n \, dV \qquad (41) $$

We subsequently combine this equation with the propagation equation (13) to obtain

$$ \phi_{t+\Delta t} = \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} (\widehat{\phi}_t^n - \Delta t V_n) \, dV \qquad (42) $$

The zero-isolevel contour is not completely maintained when the last equation is used with a vanishing velocity. However, we have found that it generally removes high-frequency oscillations that can arise from the initialization step. This property is interesting since it prevents the appearance of artificial shocks without need for auxiliary measures [9].

<caption>**Algorithm 1** Algorithm for finite element level set computations</caption>

**Require:** $\forall \mathbf{x} \in \Omega$, $P(\mathbf{x})$ and $V_n(\mathbf{x})$ known
$$ \phi_{t=0} \leftarrow \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} P(\mathbf{x}) \, dV $$
**repeat**
$$ \widehat{\phi}_t^{k+1} \leftarrow (\mathbf{A}^k + \mathbf{A}_{\mathrm{GLS}}^k)^{-1} (\mathbf{b}^k + \mathbf{b}_{\mathrm{GLS}}^k) $$
**until**
$$ \|\| \nabla \widehat{\phi}_t^{k+1} \|\| - \|\| \nabla \widehat{\phi}_t^k \|\| < \varepsilon $$
$$ \phi_{t=0} \leftarrow \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} \widehat{\phi}_t^n \, dV $$
**repeat**
$$ \widehat{\phi}_t^{k+1} \leftarrow (\mathbf{A}^k + \mathbf{A}_{\mathrm{GLS}}^k)^{-1} (\mathbf{b}^k + \mathbf{b}_{\mathrm{GLS}}^k) $$
**until**
$$ \|\| \nabla \widehat{\phi}_t^{k+1} \|\| - \|\| \nabla \widehat{\phi}_t^k \|\| < \varepsilon $$
$$ \phi_{t+\Delta t} \leftarrow \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} (\widehat{\phi}_t^n - \Delta t V_n) \, dV $$
**until** $t=T$

The method is completed by an approach to construct an initial field, $\phi(\mathbf{x}, t=0)$, from a measured discrete field, e.g. a bitmap picture, where the interface is defined as a difference in contrast. For convenience, we reduce the picture to a binary scalar function $P: \Omega \to \mathbb{R}$, which is defined by

$$ \forall \mathbf{x} \in \Omega, \quad P(\mathbf{x}) = \begin{cases} \beta & \text{if the pixel value at } \mathbf{x} \text{ is superior to } s \\ -\beta & \text{else} \end{cases} \qquad (43) $$

where $s$ is a threshold that states the limit of the contrast difference. Next, the distance between $P$ and $\phi$ is minimized on the domain $\Omega$:

$$ \forall \mathbf{x} \in \Omega, \quad \min_{\phi: \Omega \to \mathbb{R}} (\phi(\mathbf{x}) - P(\mathbf{x})) \qquad (44) $$

This expression can be cast in a variational formulation along the same lines as before and leads to a discrete (finite element) formulation given by

$$ \phi = \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} P(\mathbf{x}) \, dV \qquad (45) $$
To further improve the accuracy of the approximation of $\|\nabla\phi\| = 1$, the two values carried by $P$ are related to the median element size, i.e. $\beta = H_e/2$, where $H_e$ stands for the mean element size of the mesh.

The resulting algorithm (Algorithm 1) thus exhibits two main steps. The first step—the initialization stage—results in the initial field $\phi$ and includes a propagation step with a zero velocity. This will initiate the propagation with an already projected field. The second step is the propagation and the initialization to the signed distance. The initialization to the signed distance is done every step *a priori*, notwithstanding whether the error criterion on the gradient norm is satisfied or not.

## 4. EXAMPLES

First, a quantitative study is carried out for the example of Figure 2, where a straight discontinuity is initialized at $x = 1.5$. Quadrilateral elements with bilinear shape functions are used and a single element is employed over the height of the bar. The length of the bar is chosen as $\pi$, so as to prevent the occurrence of trivial solutions. Since there is no discretization error due to a curved front, it is possible to quantify the performance of the numerical scheme, including its dependence on the time step, the fineness of the discretization and the spatial numerical integration of the resulting vectors and matrices.

Two criteria are used to assess the performance of the numerical scheme. The computed error with respect to the position of the front is measured through

$$e_{\phi} = \frac{|x - \tilde{x}|}{L} \qquad (46)$$

where $x$ and $\tilde{x}$ stand for the exact and numerically computed positions of the front, with $L$ the length of the bar. For the error in the norm of the gradient, we use

$$e_{\nabla} = \frac{\int_{\Omega} \|\nabla\phi\| - 1 |dV|}{\int_{\Omega} dV} \qquad (47)$$

Factors that influence the error with respect to the position of the front at initialization are the discretization and the number of Gauss points in the elements that are crossed by the zero-isolevel contour. The latter is important since this number determines the accuracy of the sampling of the binary scalar function $P(\mathbf{x})$, Equation (43), and therefore provides a spatial filter to the image. Figure 3 shows the influence of these two factors. Clearly, a minimum number of 64 Gauss points in the elements, where $P(\mathbf{x})$ has different values, are needed so that mesh refinement results in an error that is monotonically decreasing.

<figure><img src="image_3.png" /><figcaption>Figure 2. Initial configuration for the first test case.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. Initialization error in terms of $e_{\phi}$ as a function of the discretization for different numbers of Gauss points.</figcaption></figure>

Figure 4 presents the evolution of the error with respect to the signed-distance function at initialization, measured by $e_{\nabla}$, for progressive iterations. Graphs are presented for different reference domains around the initial position of the discontinuity, namely $\Omega^0$ (left) and $\Omega^0 \cup \Omega^1$, see Figure 1 for the definitions of the domains, in order to bring out the quality of convergence on elements that contain the front as well as on adjacent elements.

Next, the properties of the finite element schemes are studied for a front that is propagating with a constant unit velocity. Figure 5 shows that the error is not amplified irrespective of the size of the elements. Furthermore, the classical Courant-Friedrichs-Lewy (CFL) criterion holds as is shown in Figure 6. Indeed, this figure shows that for time steps at or below the critical time step defined as

$$ \Delta t_{CFL} = \min_{x \in \Omega} \frac{h_e}{V_n} \qquad (48) $$

with $h_e$ being the size of the element that contains the point with spatial coordinates $x$, the error is conserved, but the error is amplified when the CFL criterion as formulated above is violated.

The second example concerns a two-dimensional computation of two expanding circular contours. Owing to the increased complexity of the problem, a quantitative assessment is no longer possible. The initial positions of the circles, which are used for the initialization of the level set field, are shown in Figure 7. After initialization, the velocity field is chosen constant and the time step depends on the discretization according to the CFL criterion, Equation (48).

When linear polynomials are used as the finite element shape functions, the projection space for the signed-distance field is rather poor. Indeed, the only signed-distance function that can be represented properly is that which represents straight lines. Figure 8 underscores this, since the performance in terms of a neat representation of the curved contour is not impressive, although the results of course improve with mesh refinement. The results show a marked improvement when quadratic polynomials are used as the finite element base functions, Figure 9. Indeed, the results
<figure><img src="image_5.png" /><figcaption>Figure 4. Error in the signed-distance function: $e_∇$ as a function of the number of iterations for different discretizations. Top: for initialization on $\Omega^0$ and bottom: for initialization on $\Omega^0 \cup \Omega^1$.</figcaption></figure>

obtained for this case on a 10 × 10 mesh are already better than those computed for the finest mesh (40 × 40) with linear interpolants.

## 5. CONCLUDING REMARKS

In this contribution we have presented a finite element approach to solve the evolution equation for level set functions. The first element of the approach is the extension of the velocity field
<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of the error $e_{\phi}$ during propagation of the front as a function of the discretization for a constant velocity.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Evolution of the error $e_{\phi}$ during propagation of the front as a function of the time step for a fixed discretization.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7. Initial position of the two circular contours.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Expanding circles at a constant velocity for different discretizations with linear shape functions. From left to right: 10×10, 20×20, and 40×40 elements.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Expanding circles at a constant velocity for different discretizations with quadratic shape functions. From left to right: 10×10, 20×20, and 40×40 elements.</figcaption></figure>
that is normally known only at the zero-isolevel contour to the entire domain of interest. This
is done via a standard orthogonality requirement between the gradients of the level set function
and the extended velocity field, which is tantamount to the requirement that the norm of the
gradient of the level set function maintains a unit value. The latter equation is solved throughout
the domain using a finite element method. Secondly, a unit constant gradient norm is maintained
in the propagation equation. Prior to the propagation step a re-initialization to the signed distance
is performed for this. Here, a major difficulty is to enforce the vanishing of the level set function at
the zero-isolevel contour during the re-initialization. Indeed, this requirement induces a Dirichlet
condition at an internal boundary, for which a straightforward solution was not readily available.
Herein, this internal boundary condition has been imposed by exploiting the partition-of-unity
property of finite element shape functions, which enables the imposition of a possibly non-zero
internal boundary condition of the Dirichlet type in a rigorous and elegant manner. The resulting
equations can be solved in a manner that is now well established for the partition-of-unity method.
To improve stability the discretized equations are augmented by Galerkin least-squares terms. The
accuracy of the method is assessed for a one-dimensional example, while a qualitative example
that involves expanding circular contours is given to demonstrate the versatility of the method in
two-dimensional applications with curved discontinuities.

REFERENCES

1. Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulation. *Journal of Computational Physics* 1988; **79**:12–49.

2. Sethian JA. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Sciences. Cambridge University Press: Cambridge, 1998.

3. Gravouil A, Moës N, Belytschko T. Non-planar 3D crack growth by extended finite elements and level sets—Part I: mechanical model. *International Journal for Numerical Methods in Engineering* 2002; **53**:2549–2568.

4. Gravouil A, Moës N, Belytschko T. Non-planar 3D crack growth by extended finite elements and level sets—Part II: level set update. *International Journal for Numerical Methods in Engineering* 2002; **53**:2569–2586.

5. Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:6183–6200.

6. Duflot M. A study of the representation of cracks with level sets. *International Journal for Numerical Methods in Engineering* 2007; **70**:1261–1302.

7. Duddu R, Bordas S, Chopp D, Moran B. A combined extended finite element and level set method for biofilm growth. *International Journal for Numerical Methods in Engineering* 2007. DOI: 10.1002/nme.2200.

8. Barth T, Sethian JA. Numerical schemes for the Hamilton–Jacobi and level set equations on triangulated domains. *Journal of Computational Physics* 1998; **145**:1–40.

9. Chessa J, Smolinski P, Belytschko T. The extended finite element method for solidification problems. *International Journal for Numerical Methods in Engineering* 2002; **53**:1959–1977.

10. Mourad HM, Dolbow J, Garikipati K. An assumed-gradient finite element method for the level set equation. *International Journal for Numerical Methods in Engineering* 2005; **64**:1009–1032.

11. Moës N, Béchet E, Tourbier M. Imposing essential boundary condition in the extended finite element method. *International Journal for Numerical Methods in Engineering* 2006; **67**:1641–1669.

12. Hughes TJR, Franca LP, Hulbert GM. A new finite element formulation for computational fluid dynamics. Part VIII: the Galerkin least square method for advective–diffusive equations. *Computer Methods in Applied Mechanics and Engineering* 1989; **73**:173–189.

13. Zhao H-K, Chan T, Merriman B, Osher S. A variational level set approach to multiphase motion. *Journal of Computational Physics* 1996; **127**:179–195.

14. Gomes J, Faugeras O. Reconciling distance functions and level sets. *Journal of Visual Communication and Image Representation* 2000; **11**:209–223.

15. Babuska I, Melenk JM. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**:727–758.
16. Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 2001; **50**:993–1013.

17. Réthoré J, de Borst R, Abellan MA. A two-scale approach for fluid flow in fractured porous media. *International Journal for Numerical Methods in Engineering* 2007; **71**:780–800.

18. Fries TP. A corrected XFEM approximation without problems in blending elements. *International Journal for Numerical Methods in Engineering* 2007. DOI: 10.1002/nme.2259.