# A 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 Finite Element Method for Level Sets. ECCOMAS Multidisciplinary Jubilee Symposium, 2008, Vienne, Austria. pp.95-106, 10.1007/978-1-4020-9231-2_7. hal-00400233

HAL Id: hal-00400233

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

Submitted on 26 Mar 2021

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 CC BY 4.0 - Attribution - International License
# A Finite Element Method for Level Sets

S. Valance, R. de Borst, J. Réthoré, and M. Coret

**Abstract** Level set methods have recently gained much popularity to capture discontinuities, including their possible propagation. In this contribution we present a finite element approach for solving the governing equations of level set methods. After a review of the governing equations, the initialisation of the level sets, the discretisation on a finite domain and the stabilisation 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.

**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 propagation, especially in

S. Valance, J. Réthoré, and M. Coret
LaMCoS, UMR CNRS 5514, INSA de Lyon, 69621 Villeurbanne, France

R. de Borst
Eindhoven University of Technology, Department of Mechanical Engineering, P.O. Box 513, 5600
MB Eindhoven, The Netherlands; E-mail: r.d.borst@tue.nl
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 somewhat 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 utilised 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 initialisation 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 initialisation of the level sets, and the discretisation 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 nonsymmetric character, stability and uniqueness are not ensured. For this reason a stabilisation term is added using a Galerkin least-square formalism [12]. The ensuing algorithm is presented and the section is concluded by a description how to initialise the level set method from measured (discrete) data.

## 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) $$
**Fig. 1** Domains around the internal discontinuity Γ with:
$\Omega^\alpha = \Omega_+^\alpha \cup \Omega_-^\alpha$

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

or, equivalently,

$$ \frac{\partial \phi}{\partial t} + \nabla \phi \cdot \frac{\partial \mathbf{x}}{\partial t} = 0. \qquad (2) $$

Equation (2) defines the propagation of the zero-isolevel contour, or equivalently, of the evolving discontinuity Γ. We next define the velocity $V_n$ that is normal to the discontinuity Γ, Fig. 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, \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, \left\{ \begin{array}{l} \frac{\partial \phi}{\partial t} + V_n = 0 \\ \phi(\mathbf{x}, t) = 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

$$ \left\{ \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} \right. \qquad (7) $$

with $V_e$ the extended velocity field [10, 13]. The orthogonality requirement implied in Eq. (7a) 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 discretisation, we multiply Eq. (6) by a test function $\psi$, integrate over $\Omega$, and use Eq. (7) 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$, and using a Bubnov-Galerkin method in the sense that the test function $\psi$ is also interpolated in the sense of Eq. (9), Eq. (8) becomes:

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

with the mass matrix $\mathbf{M}$ defined as

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

Using a first-order Euler time discretisation 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 \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{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 $||∇φ|| = 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:

$$ \left\{ \begin{array}{ll} ||∇φ|| = 1 & ⇐⇒ ∇φ · ∇φ_t = 1, \quad \forall \mathbf{x} \in Ω/Γ \\ φ = 0, & \forall \mathbf{x} \in Γ \end{array} \right. \qquad (14) $$

Since the first equation of the set (14) is nonlinear, an iterative procedure at each time step is necessary for its solution. Adopting a first-order Taylor series at $∇φ_t$ one obtains:

$$ g(∇φ_t + ε_φ) ≈ g(∇φ_t) + \frac{∂g}{∂∇φ}(∇φ_t)ε_φ \qquad (15) $$

with

$$ g(∇φ_t) = ∇φ_t · ∇φ_t - 1 \qquad (16) $$

Considering that $ε_φ ≪ 1$, a Newton-Raphson procedure can now be used. At iteration $k+1$ a field $φ_t^{k+1}$ is sought, such that:

$$ \left\{ \begin{array}{ll} 2∇φ_t^k · ∇φ_t^{k+1} = 1 + ∇φ_t^k · ∇φ_t^k, & \forall \mathbf{x} \in Ω/Γ \\ φ_t^{k+1} = 0, & \forall \mathbf{x} \in Γ \end{array} \right. \qquad (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(Ω)$, the weak formulation of Eq. (17) becomes:

$$ \forall ψ \in U = \{φ \in \mathcal{W} | φ = 0 \text{ on } Γ\} \qquad (18) $$

find $φ_t^{k+1}$ such that

$$ \mathcal{A}^k(ψ, φ_t^{k+1}) = ℓ^k(ψ) \qquad (19) $$

with

$$ \mathcal{A}^k(ψ, φ_t^{k+1}) = 2 \int_{Ω} ψ ∇φ_t^k · ∇φ_t^{k+1} dV \qquad (20) $$

$$ ℓ^k(ψ) = \int_{Ω} ψ(1 + ∇φ_t^k · ∇φ_t^k) dV \qquad (21) $$

A major difficulty regarding the solution of Eq. (19) arises from the essential condition that is imposed at the internal boundary $Γ$, cf. Eq. (14). General approaches to impose boundary conditions, or linear relations between degrees-of-freedom, in finite element methods are either to reduce the projection space 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}^{\mathcal{N}} N_i(\mathbf{x}) \left( \bar{\phi}_i + \sum_{j=1}^{m} \tilde{N}_j(\mathbf{x}) \tilde{\phi}_{ij} \right) \quad (22) $$

with $N_i(\mathbf{x})$ the traditional polynomial shape functions, $\tilde{N}_j(\mathbf{x})$ the enrichment functions at node $i$, $\tilde{\phi}_i$ and $\tilde{\phi}_{ij}$ the degrees-of-freedom that relate to the standard and the enhanced interpolations, respectively, $\mathcal{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 \quad (23) $$

which represents the value of the level set function $\phi$ at the initialisation 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 \quad (24) $$

and substitution into Eq. (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 \quad (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 \quad (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 write Eq. (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 $x$ has been dropped. The enrichment is only used in the area of interest, i.e. at the nodes that belong to $\Omega^0$, Fig. 1. A well-known issue of partition-of-unity based finite element methods is that the enrichment should be regularised outside the area of interest in order to avoid perturbations [16, 18]. For the present purpose, we only use 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 discretisation (27), Eq. (19) can be written in a standard matrix-vector format, as follows:

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

with

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

$$ b^k = \int_{\Omega} \left( 1 + (\nabla N_e \hat{\phi}_t^k)^T (\nabla N_e \hat{\phi}_t^k) \right) N_e^T dV \qquad (30) $$

### 3.3 Stabilisation

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

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

find $\phi_t^{k+1}$ such that:

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

with

$$ \mathcal{A}_{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_{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 stabilisation parameter

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

In Eqs. (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 discretisation (27), Eq. (32) can be written in a standard matrix-vector format, as follows:

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

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

$$ \mathbf{A}_{\text{GLS}}^k = \sum_E 4 \tau^e \int_{\Omega_e} \nabla \mathbf{N}_e^T (\nabla \mathbf{N}_e \hat{\boldsymbol{\phi}}_t^k)^T (\nabla \mathbf{N}_e \hat{\boldsymbol{\phi}}_t^k) \nabla \mathbf{N}_e \, dV \qquad (37) $$

$$ \mathbf{b}_{\text{GLS}}^k = \sum_E 2 \tau^e \int_{\Omega_e} \left( 1 + (\nabla \mathbf{N}_e \hat{\boldsymbol{\phi}}_t^k)^T (\nabla \mathbf{N}_e \hat{\boldsymbol{\phi}}_t^k) \right) \nabla \mathbf{N}_e^T (\nabla \mathbf{N}_e \hat{\boldsymbol{\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{\boldsymbol{\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$, Fig. 1.

## 3.4 The Complete Algorithm

As the preceding re-initialisation 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-initialisation step, $\hat{\phi}_t^n$, onto the classical Galerkin space:

$$ \forall \mathbf{x} \in \Omega, \min(\hat{\phi}_t - \phi_t) \qquad (40) $$

In a discrete format we then have:

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

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

$$ \boldsymbol{\phi}_{t+\Delta t} = \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^T (\hat{\boldsymbol{\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 initialisation step. This property is interesting since it prevents the appearance of artificial shocks without need for auxiliary measures [9].

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}$, defined by:

$$ \forall \mathbf{x} \in \Omega, 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 minimised on the domain $\Omega$:

$$ \forall \mathbf{x} \in \Omega, \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}) \mathrm{d}V \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 initialisation stage – results in the initial field $\phi$ and includes a propagation

**Algorithm 1** Algorithm for finite element level set computations

**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}) \mathrm{d}V $$

**repeat**

$$ \hat{\phi}_{t}^{k+1} \leftarrow (\mathbf{A}^{k} + \mathbf{A}_{\mathrm{GLS}}^{k})^{-1} (\mathbf{b}^{k} + \mathbf{b}_{\mathrm{GLS}}^{k}) $$

$$ \mathbf{until} \quad \left\| \nabla \hat{\phi}_{t}^{k+1} \right\| - \left\| \nabla \hat{\phi}_{t}^{k} \right\| < \varepsilon $$

$$ \phi_{t=0} \leftarrow \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} \hat{\phi}_{t}^{n} \mathrm{d}V $$

**repeat**

$$ \hat{\phi}_{t}^{k+1} \leftarrow (\mathbf{A}^{k} + \mathbf{A}_{\mathrm{GLS}}^{k})^{-1} (\mathbf{b}^{k} + \mathbf{b}_{\mathrm{GLS}}^{k}) $$

$$ \mathbf{until} \quad \left\| \nabla \hat{\phi}_{t}^{k+1} \right\| - \left\| \nabla \hat{\phi}_{t}^{k} \right\| < \varepsilon $$

$$ \phi_{t+\Delta t} \leftarrow \mathbf{M}^{-1} \int_{\Omega} \mathbf{N}^{\mathrm{T}} (\hat{\phi}_{t}^{n} - \Delta t V_{n}) \mathrm{d}V $$

**until** $t = T$
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Initial position of the two circular contours</figcaption></figure>

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

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

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Expanding circles at a constant velocity for different discretisations with linear shape functions. From left to right: 10 × 10, 20 × 20, and 40 × 40 elements</figcaption></figure>

step with a zero velocity. This will initiate the propagation with an already projected field. The second step is the propagation and the initialisation to the signed distance. The initialisation to the signed distance is done every step *a priori*, notwithstanding whether the error criterion on the gradient norm is satisfied or not.

# 4 Example

As example we show the computation of two expanding circular contours. The initial positions of the circles, which are used for the initialisation of the level set field are shown in Fig. 2. After initialisation, the velocity field is chosen constant and the time step depends on the discretisation according to the CFL criterion.

When linear polynomials are used as the finite element shape functions, the projection space for signed distance field is rather poor. Indeed, the only signed distance function that can be represented properly is that which represents straight lines. Figure 3 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, Fig. 4. Indeed, the results obtained for this case on a 10 × 10 mesh are already better than those computed for the finest mesh (40 × 40) with linear interpolants.
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Expanding circles at a constant velocity for different discretisations with quadratic shape functions. From left to right: 10 × 10, 20 × 20, and 40 × 40 elements</figcaption></figure>

# 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 that is normally known only at the zero-isolevel contour to the entire domain of interest. This 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-initialisation 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-initialisation. 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 fashion. The resulting equations can be solved in a manner that is now well established for the partition-of-unity method. To improve stability the discretised equations are augmented by Galerkin-Least Squares terms. 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 (1988) Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation. *Journal of Computational Physics* **79**: 12–49.

2. Sethian JA (1998) *Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Sciences*. Cambridge University Press, Cambridge.
3. Gravouil A, Moës N, Belytschko T (2002) Non-planar 3D crack growth by extended finite elements and level sets – Part I: Mechanical model. *International Journal for Numerical Methods in Engineering* **53**: 2549–2568.

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

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

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

7. Duddu R, Bordas S, Chopp D, Moran B (2008) A combined extended finite element and level set method for biofilm growth. *International Journal for Numerical Methods in Engineering*. **74**: 848–870.

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

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

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

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

12. Hughes TJR, Franca LP, Hulbert GM (1989) 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* **73**: 173–189.

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

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

15. Babuska I, Melenk JM (1997) The partition of unity method. *International Journal for Numerical Methods in Engineering*, **40**: 727–758.

16. Belytschko T, Moës N, Usui S, Parimi C (2001) Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* **50**: 993–1013.

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

18. Fries, TP (2008) A corrected XFEM approximation without problems in blending elements. *International Journal for Numerical Methods in Engineering* **75**: 503–532.