A high-order elliptic PDE based level set reinitialisation method using a discontinuous Galerkin discretisation

In this paper, an eﬃcient, high-order accurate, level set reinitialisation method is proposed, based on the elliptic reinitialisation method (Basting and Kuzmin, 2013 [1]), which is discretised spatially using the discontinuous Galerkin (DG) symmetric interior penalty method (SIPG). In order to achieve this a number of improvements have been made to the elliptic reinitialisation method including; reformulation of the underlying minimisation problem driving the solution; adoption of a Lagrange multiplier approach for enforcing a Dirichlet boundary condition on the implicit level set interface; and adoption of a narrow band approach. Numerical examples conﬁrm the high-order accuracy of the resultant method by demonstrating experimental orders of convergence congruent with optimal convergence rates for the SIPG method, that is h p +1 and h p in the L 2 and DG norms respectively. Furthermore, the degree to which the level set function satisﬁes the Eikonal equation improves proportionally to h p , and the often ignored homogeneous Dirichlet boundary condition on the interface is shown to be satisﬁed accurately with a rate of convergence of at least h 2 for all polynomial orders.


Introduction
The level set method is a popular technique used for representing and tracking evolving interfaces in computer simulations which has found use across a wide range of areas interesting to computational physicists and engineers. This includes applications in fluid dynamics [2], shape optimisation [3], computer vision [4] and biomechanics [5] to name just a few; an extensive review into the methods surrounding and further applications of the level set method can be found in the textbooks written on the subject, in particular those by Sethian [4], and Osher and Fedkiw [6]. The aim of the work to be presented in this article is to extend the level set methodology, through the development of a level set reinitialisation method which employs a discontinuous Galerkin (DG) spatial discretisation. The remainder of the introduction is divided into the three following sections: 10 Section 1.1 provides an introduction to the level set method, Section 1.2 provides an introduction to DG methods, and finally Section 1.3 provides a review of the level set reinitialisation literature.

Level set method
The level set method was originally developed by Osher and Sethian [7], in 1988. The idea behind the level set method is to use a real scalar valued function, φ, called a level set function to 15 divide a problem domain, Ω, into a number of subdomains. In general, there will be two subdomains and an interface between the two which can be defined by the value of the level set function at any point in the domain. This can be stated as φ > 0 in Ω\D, φ = 0 on Γ and where Ω is the problem domain, Γ denotes the level set function's zero isocontour also known as the interface (i.e. the interface between the subdomains) and D is a subdomain implied by the level 20 set function. An example of a circular interface defined by a level set function in a square domain can be seen in Figure 1.
The level set function can be evolved through the solution of a scalar transport problem, sometimes called the level set equation, which can be stated as where t e is pseudotime as related to the evolution equation and b is the interface velocity. 25 When evolving a level set function through the solution of (2), the level set interface can only be transported along its normal, and as such it is a natural choice to initialise the level set function  positive for φ ∈ Ω\D and negative for φ ∈ D, using the notation from equation (1). One property 30 of a signed distance function is that it will satisfy the Eikonal equation, which can be stated as The example level set function shown in Figure 1(b), is a signed distance function to the circular interface.
For a given velocity field, b, it is unlikely that after any given time step in the solution of the evolution equation, that the level set function will maintain the properties of a signed distance 35 function. Whilst it is not required that the level set function must maintain these properties, it is often preferred, as it has been shown that large variations in the gradient of the level set function, can cause numerical instability during the solution of the level set evolution problem, [8]. The desire to maintain the level set function as a signed distance function led Chopp [9], to introduce the idea of reinitialisation, by which, between iterations of the level set evolution problem, the level 40 set function can be reinitialised as a signed distance function to the new interface. Reinitialisation makes it possible to ensure that for all time, t e , the level set function is, as defined by the user, a 'good' approximation of a signed distance function, and therefore allows one to generate numerically stable results.
3 In terms of a full level set methodology it can be observed that there is a further advantage to 45 always ensuring that the level set function satisfies the Eikonal equation. Given that the advection velocity can be written, b = bn φ where, n φ = ∇φ |∇φ| , is the normal of the level set function, and b, is the scalar magnitude of the advection velocity normal to the interface, then the evolution equation (2) can be simplified as follows ∂φ ∂t e = −b|∇φ| = −b.
In this way, much of the complexity in solving the evolution equation is now exchanged for the 50 complexity of solving the reinitialisation problem (with the potential added expense of having to reinitialise more often).

Discontinuous Galerkin methods
DG methods are a class of non-conforming finite element methods where the test and trial functions are not continuous across the faces and edges of the mesh, [10]. This decoupling be-55 tween adjacent elements allows the methods to be both highly parallelisable and well suited for hp−adaptivity, which in turn should allow one to develop methodologies which are both efficient and high-order accurate, [11]. This would be advantageous for use with some of the more expen- outlined in the conclusions of [13]. In particular, we highlight the well established optimal rates of convergence which are not possible for the nonsymmetric DG methods (which includes the nonsymmetric interior penalty method (NIPG), [14], and the method of Baumann-Oden (BO), [15] It should be noted that the aim of any of the reinitialisation methods to be discussed below is to ensure that at the beginning of each iteration of the solution of the level set equation, the level set function is a signed distance function to the current position of the level set interface; i.e. the interface returned at the end of the previous iteration of the evolution equation. In this sense, all of 90 the reinitialisation techniques discussed below are equivalent, however, they vary in terms of their computational efficiency, stability and accuracy, especially when applied to a DG level set methods.
The original reinitialisation method introduced by Chopp [9], reinitialises the level set function as a signed distance function using a direct geometric approach. This approach works by first explicitly discretising the interface, Γ h ∼ Γ, and then at each point in the problem domain setting 95 the value of the level set function equal to the minimum distance from that point, to the discretised interface multiplied by the sign of the original level set function at that point, which can be stated where φ 0 is the pre-reinitialisation level set function, and sign(·) denotes the signum function.
Whilst conceptually simple, there are a number of issues with the geometric reinitialisation 100 method. One of the key advantages of the level set method in terms of computational efficiency, is the implicit nature of the evolving interface. Not only is this advantage surrendered by discretising the interface, but the expense required to both discretise the interface and compute the minimum distance at each mesh node, to the discrete interface, increases with mesh density, with the number of points used to discretise the interface and with the length of the interface itself. Chopp notes 105 that the complexity of this reinitialisation method is O(n 6 ), [9]. Furthermore, adaptation of such an approach to a discontinuous Galerkin discretisation also poses some additional difficulties. As there is no longer a requirement of continuity across element edges, the zero isocontour can be discontinuous and thus the computation of the distance from a point to the interface can become problematic. Similarly, the sign of the level set function at each degree of freedom for a given 110 node will not necessarily be well-defined, particularly if a node is near to the interface. These problems will either cause strong discontinuities to develop in the level set function, or lead to a smoothing of the level set function which will cause movement in the position of the interface post-reinitialisation. Lastly, when using the geometric reinitialisation method, the approximation of the interface on each element will only be first-order, and any benefits arising from the high-115 order approximations possible through the use of DG methods will be surrendered each time the reinitialisation routine is called.
One of the most popular methods of reinitialisation is a PDE based, pure reinitialisation method, often referred to as the hyperbolic reinitialisation method, which was introduced by Sussman et al. [17]. This method involves solving a hyperbolic PDE which can be stated as where t r denotes pseudotime as related to the reinitialisation problem. The steady state solution to (6) will be achieved once the level set function provides a sufficient approximate solution to the Eikonal equation, (3). The multiplication by the sign of the pre-reinitialisation level set function, φ 0 , works as a weak Dirichlet boundary condition on the interface of the level set function, and thus the resulting function will be a signed distance function to the interface, Γ(φ 0 ).

125
The main issue with the hyperbolic reinitialisation method is that the (potentially poor) characteristics of the original level set function, φ 0 , can be propagated during the reinitialisation process, which most often presents as a 'smearing' of the interface [18]. Mousavi [19], presented a solution to (6), using a DG discretisation along with a third-order Runge-Kutta scheme in time. Mousavi outlines quite clearly the difficulties encountered in trying to produce a stable solution using these 130 methods of discretisation; the results presented demonstrate that at some point in pseudotime the solution will always gradually begin to diverge. Independent work done by the author of this paper, instead using an explicit Euler discretisation in time, found a similar issue when trying to solve the hyperbolic reinitialisation problem using a spatial DG discretisation. Mousavi [19] found that it was possible to create a method which was practically viable by utilising a severe time step restriction, a 135 sufficiently smoothed signum function and including an artificial viscosity term. Such a solution to the reinitialisation problem is less than ideal however, as a large number of iterations are required to return a signed distance function everywhere in the domain, which could be considered prohibitively expensive. Similar issues were found by Karakus et al. [20], in which the author takes advantage of the high level of parallelisation possible with DG methods to speed up the computation of the such that it was no longer a Hamilton-Jacobi equation, thus developing the first all-in-one type 145 method. Whilst theoretically such a formulation should force the level set function to maintain it's signed distance properties, it was found that once discretised there could still be a drift in the level sets leading to a loss of the signed distance property over time [22]. This idea however, prompted other all-in-one type methods whereby the evolution equation is modified to include a signed distance constraint, such that at each time step the resulting level set function is a solution 150 to both the evolution problem and the Eikonal equation. For example, Weber et al. [22], set up their evolution problem as an optimisation problem driven by an error functional which minimises deviations in the desired interface movement and also deviations from the signed distance property.
A similar solution was presented by Li et al. [23] whereby the level set evolution problem was reframed as an optimisation problem including an energy driving the evolution and a penalty term where α is a penalty parameter. Later, Li et al. [24], named this, distance regularised level set evolution (DRLSE).
Basting and Kuzmin [1], took the distance regularisation part of the DRLSE, and considered it 160 as a pure reinitialisation problem, which is a parabolic PDE and can be stated as By removing the time dependent part, so as to avoid pseudotime stepping, and also including an appropriate boundary condition, Basting and Kuzmin reformulated the problem as a quasilinear elliptic PDE to be solved iteratively which can be stated as where γ is a penalty parameter. The work presented in this paper provides a solution to the elliptic 165 reinitialisation problem using a DG method for the spatial discretisation.
Whilst the work presented in this paper was completed independently, Utz et al. [25] recently presented a similar DG solution to the elliptic reinitialisation problem. However, a number of issues were found with the work presented in [1] and [25], solutions to which are discussed here. Explicitly, in this paper issues are addressed concerning: boundary conditions on an implicit surface; experi-170 mental orders of convergence which align with the theoretically optimal rates of convergence in the relevant norms through the use of a narrow band approach; and the construction of a new potential function which removes the issues when reinitialising level set functions with small gradients, i.e. |∇φ| ≤ 0.5.

175
The rest of the paper is organised as follows. Section 2 presents the proposed elliptic reinitialisation method. Section 3 presents three numerical examples as demonstrations of the efficacy of the proposed method. The article is then concluded in Section 4.

Discontinuous Galerkin elliptic level set reinitialisation
Section 2 consists of the following subsections. Section 2.1 presents the mathematical prelimi-180 naries required for the discussion of DG methods. Section 2.2 presents the elliptic reinitialisation problem and the proposed DG discretisation. Section 2.3 presents a discussion on the use of a narrow band approach, which is required to allow one to demonstrate optimal rates of convergence. 8

Symmetric interior penalty discontinuous Galerkin method preliminaries
Let T h denote any partition of a domain, Ω, into nonoverlapping quadrilateral elements, τ , with element size, h, such that the computational domain can be defined, Ω = ∪ τ ∈T h τ , with boundary vertices on ∂Ω. The skeleton of the mesh, S, is defined as the set of all interior edges, that is The unit outward normal on the boundary, ∂τ , of a given element, τ , is denoted asn. For any mesh T h of Ω, with elements of maximum polynomial degree, p, the DG finite element space is defined as where Q p (τ ) denotes the space of polynomials of degree no more than, p, in each coordinate direction.
It should be noted that the work to be presented here is restricted to regular quadrilateral elements, on Cartesian grids. This is due to the Eulerian framework within which the level set method operates, which allows one to exploit the simplicity of such an approach. 195

Elliptic level set reinitialisation
The reinitialisation problem can be stated, for a given level set function, φ 0 , find a new level set function, φ, which is a signed distance function to the original position of the level set interface, Γ(φ 0 ). This can be stated mathematically as finding a solution to the Eikonal equation, stated in equation (3), relative to the following Dirichlet boundary condition As first presented by Basting and Kuzmin in [1], the elliptic reinitialisation method aims to solve the level set reinitialisation problem by minimising the least squares residual to the Eikonal Taking the derivative of the objective functional (13), leads to a strong formulation of the problem which can be stated as 9 The first equation forming (14) Applying a Picard linearisation to the terms which are nonlinear with respect to ∇φ, allows one to rewrite the above diffusion equation as where the superscript denotes the m th iteration. Discretising the problem spatially using the SIPG method then leads to a variational formulation which can be stated as; find φ m h ∈ V DG , as m → ∞, such that the following weak form statement of equilibrium is satisfied were µ is a penalty parameter, henceforth referred to as the discontinuity penalisation parameter, which for elliptic problems can be chosen as µ = Cp 2 e /h e , where C is a constant usually equal to 10, p e is the maximum polynomial order of the two elements sharing that edge, and h e is the length of the edge. For further information on the discontinuity penalisation parameter in quasilinear elliptic problems, please refer to [26]. The jump and average operators denoted by [[·]] and {{·}} respectively, are as defined in [10] and are reproduced here as follows: for an arbitrary scalar valued function, ψ, and vector valued function, Ψ, on adjacent elements, τ + and τ − , which share an edge The resulting linear system can then be solved using a fixed point iterative method as follows; where the matrix K = (k ij ), has elements given by and the column vector F = (f i ), has elements given by It should be noted that the above formulation is incomplete as it does not enforce the Dirichlet 220 boundary condition. This will be discussed separately in Section 2.2.2. The homogeneous Neumann boundary condition is, however, naturally satisfied in the above formulation.
Modifications to the above formulation will be discussed in the following subsections. Section where Choosing the objective functional to be equal to t 1 , leads to a diffusion term in the weak formulation which could be stated as where It can be seen that the diffusion functional, d 1 , becomes singular as |∇φ| → 0. To combat this 240 problem, authors such as Li [24] and Basting [1], have modified the objective functional such that it minimises the least squares residual to the Eikonal equation everywhere except in the region where |∇φ| is small. For example, [1] presents the following functional which leads to a diffusion term shows that for the corresponding diffusion functional, d 2 , that where the gradient is small, i.e.
|∇φ| < 0.5, the diffusion is positive, which corresponds to forcing the level set function towards the solution at |∇φ| = 0.

250
In order to overcome these two issues, here we propose a new objective functional which both avoids the singularity at |∇φ| = 0 and always has negative diffusion for |∇φ| < 1. One such  functional could be stated as which leads to a diffusion term It should be stated that conceptually any function which satisfies these conditions would suffice.
255 Figure 2 demonstrates that the objective functional, t 3 , does indeed satisfy both of these conditions.
To include any of the above defined diffusion functionals using the formulation presented in Section 2.2, the only modification required to the linear system stated in (20) is the entries to the F vector, which can now be written  When using the objective functional, t 1 , it can be observed that the solution immediately begins to oscillate and does not converge. Figure 3(b) shows a snapshot of the level set function after 50 iterations when using t 1 . It can be seen that the attempt to correct the almost zero gradients in the centre element, leads to an overcorrection causing the level set function to twist as it tries to force the gradient back to unity, after which the solution breaks down and continues to get worse the analytical solution as far as possible given the coarseness of the mesh. Therefore the objective functional adopted in this work is that defined as t 3 .

Boundary conditions on implicit surfaces
In both [1] and [25], the Dirichlet boundary condition on the level set interface is enforced using a penalty method. As such the weak formulation would be stated as, find φ m h ∈ V DG , as m → ∞ such that where γ is a penalty parameter, henceforth referred to as the interface penalisation parameter.     Figure 4(b), demonstrates that if the value of the interface penalisation parameter is too large, there will be boundary locking, [29], in elements intersected by the interface.
Evidence is provided in [25], which supports the idea that an appropriate choice for the value of the interface penalisation parameter for a given mesh, is equal to the discontinuity penalisation parameter, µ such that, γ = µ. Whilst it can be observed that the interface penalisation parameter 305 is problem dependent, it is not necessarily apparent that it is related to the mesh size in the same way as the discontinuity penalisation parameter. Repeating the example problem from the previous paragraph, with a mesh of linear elements, the interface penalisation parameter would therefore be computed, γ = 10p 2 h = 50. As evidenced at a glance by the solution in Figure 5 this is an appropriate value for this penalty parameter in this case. Increasing the order of the 310 elements to p = 5 causes an increase in this value to γ = 1250; Figure 5(c) shows that this value is too large and causes locking/spurious oscillations in the elements at the boundary and therefore is not appropriate. However, once again using quintic elements, but choosing γ = 50 allows one to return a solution which no longer displays locking at the boundary as shown in Figure 5(d).
The same is true when changing the number of elements used to discretise the problem. This  The literature highlights four main approaches for the imposition of Dirichlet boundary conditions on implicit surfaces; the aforementioned penalty method, Nitsche's method [30], the method of Lagrange multipliers [31], and methods involving enrichment or modification of shape functions, for example [32]. Nitsche's method is akin to the penalty method in that there is a penalty term which imposes the prescribed value on the boundary. Without re-presenting the evidence, the same 325 arguments against using the penalty method described above were found to also be true of Nitsche's method when applied to the implicit surface. The methods involving the modification of the shape functions require a priori knowledge of the position of the interface, whereas the methodology here deals with evolving and implied interfaces only, and therefore methods such as these are not appropriate in this context.

330
The method of Lagrange multipliers involves the reformulation of the weak form of the problem such that a new unknown, the Lagrange multiplier, λ, is to be solved for in addition to the original unknown, in this case the level set function, φ, such that the solution on the Dirichlet boundary is  constrained by a prescribed value. The weak form of the elliptic reinitialisation problem can thus be reformulated: find φ m h ∈ V DG and λ ∈ L, as m → ∞ such that and One of the difficulties of using such a method, is choosing the correct interpolation space for the Lagrange multipliers, L. One natural choice is choosing the space, L, as follows where that is, T Γ h denotes the subset of elements in T h which are intersected by the level set interface, Γ.

335
This means that the Lagrange multiplier space will consist of the same basis functions as the finite element space, and therefore one can solve for one Lagrange multiplier per degree of freedom on any element intersected by the interface.
When choosing the Lagrange multiplier interpolation space, it is necessary that the space is rich enough such that it contains the approximate solution, but not so large as to overconstrain 340 the problem. It is a known phenomena, [33], that boundary locking or spurious oscillations can occur when the approximation spaces V DG and L are chosen to be of equal order. Repeating the previous experiment, using a Lagrange multiplier approach to enforce the boundary condition with the Lagrange multiplier space defined as in (35) gives the results shown in in Figure 6, which confirms that such a choice will in fact lead to boundary locking.

345
In order to rectify this problem, the order of the Lagrange multiplier space has been reduced to the space of piecewise constant functions with one degree of freedom per element intersected by the interface. This can be stated as  where 1 τ is the indicator function defined as follows This choice of space means that for each element, τ ∈ T Γ h , the integral of the level set function over 350 the portion of the interface contained within that element, averages to be zero over the element.
In other words, this reduction in the order of the constraint space allows some movement to occur at the interface (limited by the size of the element), which is a sufficient relaxation to remove the boundary locking observed above and allows the boundary condition to be satisfied without affecting the signed distance property. It should be noted that even for higher order elements, i.e.

355
p ≥ 2, choosing the Lagrange multiplier space, as the space of piecewise constants, is required to ensure that there is no boundary locking.
As such, the preferred method of the author therefore for enforcing a homogeneous Dirichlet boundary condition on an implicit surface, is to use a Lagrange multiplier approach, where the Lagrange multiplier space is the space of piecewise constant functions. Using this formulation, to 360 enforce the Dirichlet boundary condition, (12), the linear system, (20), will be modified as follows where A = (a ij ) is a matrix, where the number of rows is equal to the number of elements in T Γ h and the number of columns is the total number of degrees of freedom in the problem, with elements given by and K and F are defined in equations (21) and (31) respectively. 365

Integration on immersed implicit surfaces
Regardless of the method chosen to impose a boundary condition on an immersed implicit surface, it will require a method for integrating a function on that surface. There are three general approaches found in the literature: explicit reconstruction of the interface through mesh refinement [34]; implicit reconstruction of the interface using an approximate Dirac delta function such as in 370 the original immersed boundary method, [35]; and methods which generate a new quadrature rule over the volume of an element τ ∈ T Γ h , which is equivalent to integrating an arbitrary function over the implicit surface [36,27].
As the Eulerian nature of the level set method allows one to take advantage of the use of Cartesian meshes, methods involving r-adaptivity to explicitly reconstruct the interface are not 375 appropriate in the context of this work. Such methods also suffer from extreme computational expense, especially when the desired level of accuracy is high. Methods involving the use of an approximate Dirac delta function, allow one to replace the surface integral over the interface with an equivalent volume integral weighted by the Dirac delta function. Whilst this method is simple to implement, and has found use in other works, even prompting research into high order approxi-380 mations of the delta function [37], the method depends on the global cancellation of errors over the domain. Thus such an approach has limited accuracy when working with piecewise discontinuous level set functions.
The final group of methods are able to provide arbitrarily high-order elementwise approximations of integrals on implicit interfaces. One such method presented by Müller et al. [27], involves the 385 construction of a new quadrature rule based on the solution to the moment-fitting equations [38], 21 and can be stated as follows For a given set of divergence free vector valued functions, g , an integral over the unknown interface, Γ(φ 0 ), can be transformed to an integral over the known surface, ∂τ , using the divergence theorem, which forms the RHS of (41) which can then be approximated using a standard Gauss quadrature 390 rule. Then the weights, w, for a new quadrature rule over the element, τ , using the standard 2D Gauss quadrature abscissae, p, equivalent to integrating these functions, g , over the interface can be solved for, which can then be used to compute the integral of any function over the interface.
Müller et al. [27], chose the functions, g, to be the monomial basis functions, where the derivatives g are orthonormalised using a Gram-Schmidt procedure. The maximum order of these func-395 tions, determines the number of equations, M , to be solved and it is noted that care should be taken to ensure that the number of quadrature points, N , is chosen such that the resulting linear system is underdetermined, i.e. N > M . The full details of the method can be found in [27].
The integration method presented in [27] is the preferred method of the author, with two caveats.
Firstly, it was found that the accuracy of this integration method depends heavily on the accuracy 400 with which one is able to compute the terms on the RHS of equation (41). The Heaviside function, H(−φ), in each of the integrals is present such that the integral is computed only along the part of the edge where, φ < 0. When using standard 1D Gauss quadrature along element edges, the discontinuity present in the Heaviside function is smoothed to such an extent that it becomes difficult to predict whether a given quadrature rule will be sufficient to ensure that the method is sufficiently 405 accurate, without using a (potentially prohibitively) high-order quadrature rule. As such for edges intersected by the interface, a Newton/bisection method is used to find the intersection point(s) and a standard quadrature rule is used to integrate over this newly defined interval. Secondly, as the number of quadrature points is chosen to ensure that the system is underdetermined, the linear system will likely be rank deficient and ill-conditioned. Thus the numerically stable singular 410 value decomposition approach is used to solve for the least squares solution. Any singular values, s, deemed too small, that is s < max(s)/10 12 , are removed to further improve stability.
These two choices have proven imperative in ensuring the quadrature rule produced is then able to accurately integrate a function on the interface. As a final note, whilst generally robust, this integration method is problem dependent and small perturbations in the relative position between 415 the mesh and the immersed surface will have an influence on the accuracy for a given problem.

Narrow band level set methods
When using the level set method for problems involving evolving interfaces it can be noted that the maximum amount of movement of the interface at each time step will be a known value limited by the Courant-Friedrichs-Lewy (CFL) condition, which, if the level set function is always 420 a signed distance function, will be a function of the smallest element size, h min . In other words, the evolution can only occur within a small banded region around the interface, and therefore the information about the level set function outside of this band can effectively be ignored. Narrow band strategies such as that presented in [4], can therefore be a useful tool in reducing the computational expense when using level set methods, as the computation of both the evolution problem and the 425 reinitialisation problem can be restricted to a set of elements, defined by some measure as being close to the interface.
Computational efficiency isn't the only benefit of using a narrow band approach. One of the issues with choosing the level set function to be a signed distance function, is that if the zero isocontour of the level set function, Γ(φ), has at least one loop surrounding a simply-connected 430 subdomain, D, there will always be a singularity which occurs in the level set function, this can be observed in Figure 1 for example. An added benefit of narrow band strategies is that, for a 'sufficiently refined' mesh, almost all of these areas would be far enough away from the level set interface so as to fall outside of the narrow band. When using SIPG, it is known that optimal convergence rates are a function of the smoothness of the problem [39]. Since these singularities will 435 always occur, the use of a narrow band approach is therefore necessary to allow one to demonstrate optimal orders of convergence when using an SIPG discretisation.
'Sufficient refinement' is, as of right now, a poorly defined term. In order to capture a given interface to a prescribed level of accuracy, there is some requirement on the number and order of the elements present along the interface. In the literature, where adaptive meshes are used, the 440 general refinement strategy can be stated as "split any cell whose edge length exceeds its minimum distance to the interface", [40]. Whilst the simplicity of such a strategy is attractive and will result in high levels of h-refinement close to the interface, it is unlikely that such a refinement strategy is optimal. One area of future work therefore could be to develop appropriate error estimators and 23 refinement strategies for the level set reinitialisation problem. For the purposes of this article it 445 will suffice to demonstrate that the combination of sufficient mesh refinement and a narrow band approach, are required to return optimal convergence; this will be demonstrated in Section 3.
For the purposes of this article, it is noted that when deciding on an appropriate width for the narrow band one needs to consider that in order to satisfy the CFL condition, the furthest that the interface should be able to move to maintain stability is from the element within which it currently 450 resides to one of its neighbours. As such the best case scenario for a narrow band is the union of the be considered in this article we always start with a level set function which is not a signed distance 460 function however, and as such a slightly more conservative value is used, equal to four element widths, to account for the variation in the gradient either side of the interface. Furthermore, it is noted that the boundary conditions stated in Equation (14) naturally extend to the narrow band region. Since the narrow band will contain the set of elements cut by the interface T Γ h , no change has to be made to the imposition of the Dirichlet condition presented in Equation (39). Finally, 465 whilst the set of element edges constituting the Neumann boundary, ∂Ω, will change, it is the same homogeneous Neumann condition, (15), which is to be applied on this new set of edges.

Error measures
Where the analytical solution, φ, is known, the error for the example problems in this section 470 is given in the L 2 norm which can be stated as 24 the L ∞ norm which can be stated as and the DG norm which can be stated as For elliptic problems discretised using SIPG the optimal convergence rates in the L 2 norm are known to be h p+1 , and in the DG norm, h p , [10], assuming the problem is sufficiently smooth.

475
Similarly, it has been shown that optimal convergence when using the L ∞ norm is proportional to ln(h −1 )sh p+1 , wheres = 1 for p = 1, ands = 0 otherwise, [41]. It is shown in [39], that for a problem which lacks sufficient smoothness, the convergence rates fall back equal to the linear case for all p.
When the analytical solution is not known, there are two additional error measures which can 480 demonstrate the efficacy of the reinitialisation method. The first is an error measure which measures globally, the degree to which the computed solution satisfies the Eikonal equation, that is This signed distance error measure acts similarly to the H 1 seminorm, computing the difference between measures of the gradient of the solution. As such it would be reasonable to expect optimal convergence rates to be equivalent to optimal convergence in the H 1 seminorm, which is known to 485 be h p , once again assuming sufficient smoothness.
The second of these, is a measure of the movement of the interface in the L 2 norm, which is evaluated by integrating the difference between the computed and desired value of the solution along the original position of the interface, that is which will be referred to as the interface error measure.

490
For all of the numerical experiments presented in this section the following statements are true.
The objective functional defining the problem to be solved is that defined in equation (29), i.e. t 3 .
The fixed point iterative method is considered to have converged when E SD (φ m ) − E SD (φ m−1 ) < 25 10 −8 , or the number of iterations required exceeds 1000. The Dirichlet boundary condition is 495 enforced using the Lagrange multiplier approach, with an interpolation space consisting of piecewise constant functions. The method of Müller et al. [27] is used to compute the integral along the interface, with the maximum order of the divergence free basis functions, g , equal to 10. This is much higher than that required in practice, for the problems to be presented, however, it allows as much as possible one to remove the error associated with the mesh/problem dependency of the 500 integration method and thus better evaluate the reinitialisation method.

Circular interface
The first test case presented is that of a circular interface defined initially by a level set function, φ 0 , which can be described analytically as in the domain Ω = (−2, 2) 2 . The corresponding signed distance function, and therefore the analyt-505 ical solution to the problem can thus be stated For this problem, the zero isocontour of the level set function can also be described analytically as follows, for 0 ≤ θ ≤ 2π, x = cos(θ), y = sin(θ).
As such the interface error measure will be computed using the trapezium rule, to remove any error associated with the methods for integrating over an implicit surface.      experimental orders of convergence, using the four aforementioned error measures, are congruent with those expected for a non-smooth problem.

520
The convergence rate using the interface error measure does show an increase between p = 1 and p = 2, but remains constant beyond that point. For the purposes of our discussion it is useful to observe that the presence of a singularity in the mesh, constrains the rate at which the L 2 error at the interface decreases when using high-order elements.

Circular interface with narrow band 525
For the previous example problem, the analytical solution is known to be singular, and thus the computed experimental orders of convergence are limited. In order to demonstrate optimal rates of convergence one needs to change the domain such that everywhere within the domain the solution is smooth, which, as discussed in Section 2.3, can be achieved through the use of a narrow band approach. For this somewhat trivial example, the position of the singularity is known to be at the 530 origin and thus a naive implementation of a narrow band approach, is to simply repeat the previous experiment in the domain, Ω = (−2, 2) 2 \(−0.4, 0.4) 2 , such that the singularity at the origin is removed.
The same h-convergence study is computed on the new domain leading to the results shown in Figure 8. As expected, removing the origin from the problem domain, allows the solution to 535 be smooth enough everywhere to display optimal convergence rates in all of the relevant norms.
This includes a convergence rate using the interface error measure of h p+1 , which suggests that one might expect this to be the optimal rate of convergence for this error measure. It can be noted that the quoted orders of convergence for all measures and polynomial orders are computed using the difference between the results for h = 0.4 and h = 0.05.

540
It should also be noted here that for this example the number of iterations taken to reach the convergence criterion is often few, for this simple example; for the mesh with h = 0.05 and p = 5, just 6 iterations are required.

Smooth star shaped interface
The remaining examples will be of a more arbitrary nature than the simple circle example, 545 thus the rule determining the width of the narrow band will be defined as follows: remove from the mesh any element which has a minimum absolute nodal value greater than four times the size 28 (e) Convergence in the L ∞ norm.  of the smallest element, h min . This will also mean that analytical solutions to the problems will be unknown and as such the convergence data presented will be using the signed distance and interface error measures only. The interface error will be computed using the method of Müller [27] 550 instead of the trapezium rule, and as such the error computed will be a measure of the movement of the interface from it's initial projection as opposed to the distance from the analytical solution (although in practice, calculating the error in these two ways gives similar results except for the coarsest meshes tested).
The first of the arbitrary interfaces will be a smooth six pointed star interface, shown in Figure   555 9(a), which has an initial level set function which can be defined everywhere by on a domain of maximum size Ω = (−2, 2) 2 , however for a given element size, h, the narrow band within which the reinitialisation problem is solved will be a subset of the full domain.
In this case, an h-convergence study will be computed on a sequence of Cartesian meshes    is because the criterion defining the narrow band, is yet to be sufficient to remove the part of the mesh which is singular, see Figure 9(b). As h becomes smaller, the narrow band becomes 565 narrower and the singular part of the solution is no longer part of the mesh, allowing for optimal rates of convergence. The quoted orders of convergence for all measures and polynomial orders are computed using the difference between the results for h = 0.2 and h = 0.025.
The rate of convergence for the interface error increases slightly between the meshes of linear and quadratic elements, however increasing the polynomial order of the elements beyond that, no 570 longer results in an increase in the accuracy of the solution at the interface, despite the improving gradient solution.

Multiple arbitrary interfaces
The final example to be presented consists of multiple nested interfaces, which more closely resembles an arbitrary level set function which one might encounter in practice. In particular,   where q 1 = 1.5 x 2 + y 2 − 1 + 0.8 sin arctan y The original configuration of this mesh can be seen in Figure 11(a). These curves have been chosen somewhat arbitrarily, however, considerations were made such that across the domain, the problem   The interface error is displayed in Figure 12(b). It shows almost equivalent errors for a given element size, h, regardless of polynomial order, p, with a small increase in accuracy between p = 1 and p = 2 which was also the case for the previous example. As has been the case for all of the 595 presented examples, it is difficult to explain the behaviour of this error measure for this problem.
As such we restrict our comments to the following; for all examples the demonstrated movement of the level set function at the interface is small (in comparison to other reinitialisation methods), and furthermore can be decreased predictably by controlling the element size with order ∼ h 2 .
For this example, it should also be noted that for the denser higher-order meshes, the number 600 of iterations required to satisfy the convergence criterion grows large, for this problem when p = 3 it takes an average of 920 iterations. However, it can also be noted that, for the most dense, high-order mesh tested, it takes just 5 iterations to improve the gradient solution by 3 orders of magnitude, and 34 iterations for an improvement of 4 orders of magnitude. This suggests that in practical applications of the method, it would be up to the user to decide where to strike the 605 balance between expense and accuracy.

Conclusions
A practical method for level set reinitialisation using an SIPG discretisation has been presented, based on the elliptic reinitialisation method originally presented by Basting and Kuzmin, [1]. The proposed method is able to demonstrate optimal convergence in the relevant norms and overcomes 610 a number of issues found with other similar reinitialisation techniques. This is achieved through the adoption of a Lagrange multiplier technique, with an appropriate interpolation space, for imposing a Dirichlet boundary condition on an immersed implicit boundary; through a reformulation of the problem by the introduction of a new objective functional driving the problem; and through the adoption of a narrow band approach. This reinitialisation method can be combined with a 615 much simplified level set transport problem, to create a full DG level set methodology. It was demonstrated that a combination of sufficient refinement and a narrow band approach allow one to return optimal convergence rates, as such future work will focus on the development of error estimators and strategies for driving mesh adaptivity, based on the reinitialisation problem.