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

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 ﬁnite difference methods. However, ﬁnite difference methods are less suited for irregular domains. Moreover, it seems slightly awkward to use ﬁnite differences for the capturing of a discontinuity, while in a subsequent stress analysis ﬁnite elements are normally used. For this reason, we here present a ﬁnite 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 ﬁnite domain, and the stabilization of the resulting ﬁnite 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 ﬁnite 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.


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 . 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 1 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][9][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.

LEVEL SET METHODS
In level set methods a hypersurface = (x, t) is defined on the domain of interest . The basic idea is that the intersection of this hypersurface with the zero level, i.e. = 0, defines the internal discontinuity. From a Lagrangian point of view, stationarity of the level set field requires that or, equivalently, 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 , Figure 1: for all points x at the zero-isolevel contour, and n the normal vector at x: With these definitions, and adding the initial condition (x, t = 0) = f (x), with f (x) known, the initial value problem that describes the propagation of the zero-isolevel contour becomes [1] ∀x ∈ , Mourad et al. [10] assume that ∇ = 1 is maintained everywhere in for t>0, so that the formulation of the initial value problem (5) simplifies to ∀x ∈ , A difficulty resides in the fact that the velocity V n is defined only at the internal discontinuity . A standard methodology to circumvent this difficulty is to extend the velocity field V n at the discontinuity such that 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 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].

Propagation equation
To obtain a finite element discretization, we multiply Equation (6) 1 by a test function , integrate over , and use Equation (7) 2 to obtain With a standard interpolation N being the array that contains the interpolation polynomials and / the array that contains the nodal values of ; using a Bubnov-Galerkin method in the sense that the test function is also interpolated in the sense of Equation (9), Equation (8) becomes with the mass matrix M defined as Using a first-order Euler time discretization scheme we obtain the following discrete evolution equation:

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: 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 ∇ t one obtains where Considering that 1, a Newton-Raphson procedure can now be used. At iteration k +1 a field k+1 t is sought such that Subsequently, the problem is cast into a weak format. With U the space of sufficiently regular fields, e.g. those that belong to H 1 ( ), the weak formulation of Equation (17) becomes A major difficulty regarding the solution to Equation (19) arises from the essential condition that is imposed at the internal boundary , 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][16][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: where N i (x) are the traditional polynomial shape functions, N j (x) the enrichment functions at node i, i and 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 , and we have m = 1 and which represents the value of the level set function 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 0 of (x) that is imposed at . With the usual interpolation, we can rewrite this enrichment as and substitution into Equation (22) subsequently yields We now remove the standard part from the interpolation, i.e. from which a value at can stem that is unequal to 0 t (x)+ 0 , and obtain The first term enforces a vanishing value of (x) at the zero-isolevel contour, while the second enrichment sets the value of (x) at the zero-isolevel contour equal to 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 /, we can express Equation (26) in a more compact manner, e.g. [17] = N e / where, for notational simplicity, the explicit dependence of 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 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 , and using the discretization (27), Equation (19) can be expressed in a standard matrix-vector format as follows:

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 A k [12], so that it transforms into where and the stabilization parameter 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 , and using the discretization (27), Equation (32) can be expressed in a standard matrix-vector format as follows: with A k and b k defined in Equations (21), and The set (36) can now be solved in a standard Newton-Raphson manner by repeatedly solving the set of equations The criterion used for convergence is | ∇ −1|< , with a sufficiently small number. The method is applied only over a small area around the front, i.e. 0 ∪ 1 , Figure 1.

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, n t , onto the classical Galerkin space: To further improve the accuracy of the approximation of ∇ = 1, the two values carried by P are related to the median element size, i.e. = 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 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.

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 , 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 where x andx 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 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(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(x) has different values, are needed so that mesh refinement results in an error that is monotonically decreasing.   Figure 4 presents the evolution of the error with respect to the signed-distance function at initialization, measured by e ∇ , for progressive iterations. Graphs are presented for different reference domains around the initial position of the discontinuity, namely 0 (left) and 0 ∪ 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 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 obtained for this case on a 10×10 mesh are already better than those computed for the finest mesh (40×40) with linear interpolants.

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