{'fraction_non_alphanumeric': 0.06351033219984306, 'fraction_numerical': 0.030133403086581217, 'mean_word_length': 3.6296924194720273, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Practical Phase Field Method for an Elliptic Surface PDE

We consider a diffuse interface approach for solving an elliptic PDE on a given closed hypersurface. The method is based on a (bulk) finite element scheme employing numerical quadrature for the phase field function and hence is very easy to implement compared to other approaches. We estimate the error in natural norms in terms of the spatial grid size, the interface width and the order of the underlying quadrature rule. Numerical test calculations are presented which confirm the form of the error bounds.


Introduction
Let Γ ⊂ R n+1 (n = 1, 2) be a closed hypersurface. In this paper we are concerned with a phase field approach for the numerical solution of the PDE −∆ Γ u + u = f on Γ (1.1) and more general elliptic PDEs on surfaces. Here, ∆ Γ denotes the Laplace-Beltrami operator and f is a given function on Γ. Apart from being of interest in their own right, elliptic surface PDEs may arise as subproblems in the time discretization of parabolic surface PDEs as well as in systems involving a coupling to a bulk PDE (see e.g. [18]). A major issue in the design and analysis of numerical methods for (1.1) lies in the fact that the simultaneous approximation of the PDE and of the surface Γ is required. Let us briefly review the various computational approaches that have been suggested in the literature. Further references can be found in the nice review articles [15] and [3].
In his seminal paper [14], Dziuk proposes and analyzes a method that employs continuous, piecewise linear finite elements on a regular simplicial partitioning of Γ h , a polyhedral approximation of Γ. This approach has been extended to higher order FEM spaces and higher order polynomial approximations of Γ by Demlow in [11], while an adaptive version of the method can be found in [12]. However, the construction of a regular polynomial approximation may be difficult in practice, in particular if the surface is given implicitly in terms of a level set function. The trace finite element method, proposed by Olshanskii, Reusken and Grande in [25], is based on a background mesh which induces an unfitted approximation Γ h of Γ and employs traces of bulk finite element functions. Even though Γ h is in general not regular, optimal error estimates for piecewise linear finite elements are obtained. Further developments and variants of this trace method (also called cut finite element method) can be found in [24], [9,Section 3], [28], [20], [13], [5] and [6]. In the case of a level set representation of Γ there is a class of methods that is based on extending the PDE (1.1) to an open neighborhood of Γ. Using earlier ideas of [2], Burger considers in [4,Section 2] an extension with the property that (1.1) is satisfied simultaneously on all neighboring level surfaces. This approach gives rise to a weakly elliptic bulk PDE, which is degenerate in the direction normal to the level surfaces and which can be solved numerically with the help of standard bulk finite elements. Error estimates have been derived in [4,Theorem 6], while [8] considers the problem in a narrow band of width h around Γ and provides an O(h) bound in H 1 (Γ). In both cases the corresponding error analysis is complicated by the degeneracy of the extended PDE; an extended PDE, which is uniformly elliptic, has been proposed in [7] and [26] and involves the mean curvature of Γ. A different method which leads to a uniformly elliptic bulk PDE, is obtained by considering the equation which is satisfied by a natural extension of the solution of the surface PDE. If Γ is given implicitly in terms of the signed distance function this extension is the function which is constant in normal direction, and one is led to the closest point method, see [23] for the parabolic case. In the case of a general level set function the corresponding PDE has been derived in [9], where unfitted sharp and narrow band finite element methods have been proposed and analyzed. Note that for schemes that are based on an implicit representation of the types described above the discrete surface or the boundary of a narrow band may cut arbitrarly through a bulk element. Locating these cuts and integrating over the discrete surface or partial elements is in general cumbersome. A way to circumvent these difficulties is offered by the use of a diffuse interface method. The starting point of this approach is again an extension of the surface PDE to a neighborhood of Γ, which is then localized to a thin layer of width proportional to ε with the help of a phase field function. The resulting problem can be solved using finite elements, where the geometry is now resolved by evaluating the phase field function. This approach was suggested and analyzed in [4,Section 3] in the elliptic case, and in [27] for a linear diffusion equation for a phase field function with nonlocal support. In [16], [17] and [10] a phase field function with compact support was used in the approximation of an advection diffusion equation on a moving surface. In practice, numerical integration needs to be used which now becomes an issue as estimates for the resulting error require derivatives of the phase field function, which scale with ε −1 . Our main contribution in this paper is a new, fully practical phase field method to solve (1.1) together with a corresponding error analysis in natural norms. Furthermore we shall present test calculations for hypersurfaces in two and three dimensions which confirm the form of our error bounds.

Notation and problem formulation
Let Γ ⊂ R n+1 (n = 1, 2) be a smooth, connected, compact and orientable hypersurface without boundary. In view of the Jordan-Brouwer separation theorem, Γ divides R n+1 into an interior and an exterior domain and we denote by d the signed distance function to Γ oriented in such a way that d < 0 in the interior, d > 0 in the exterior of Γ. It is well-known (see [19], Section 14.6) that there exists an open neighbourhood Ω of Γ such that d is smooth in Ω with |∇d(x)| = 1, x ∈ Ω as well as ∇d(x) = ν(x), x ∈ Γ, where ν(x) is the unit outer normal to Γ. Furthermore, the function p(x) := x − d(x) ∇d(x) assigns to every x ∈ Ω the closest point on Γ, so that whereη is an extension of η to an open neighborhood of Γ.
Let us consider the following elliptic PDE in divergence form We assume that a ij ∈ C 2 (Γ) and that A(x) = (a ij (x)) n+1 i,j=1 defines a symmetric, uniformly positive definite linear map from T x Γ into itself, so that there exists α > 0 with Furthermore, we suppose that a 0 and f belong to W 1,∞ (Γ) and that there exists α 0 > 0 such that It follows from the Lax-Milgram lemma that for every f ∈ L 2 (Γ) the PDE (2.3) has a unique weak solution u ∈ H 1 (Γ) in the sense that where dS p is the surface element of Γ. Furthermore, standard regularity theory implies that u ∈ H 2 (Γ) and In what follows we suppose that Γ is represented in implicit form, i.e. there exists a smooth function φ :Ω → R such that By choosing Ω smaller if necessary we may assume the existence of c 1 ≥ c 0 > 0 such that (2.10)

Extension
As already mentioned in the introduction our numerical approach is based on extending surface quantities and the surface PDE to a neighborhood of Γ. In what follows we abbreviate A common way to extend a given function g : Γ → R consists in setting g(x) := g( p(x)), often called the closest-point extension, and we shall use p in order to extend the data a ij , a 0 and f to a neighbourhood of Γ. However, in order to derive our scheme and in order to carry out the error analysis we require a further extension which is better adapted to the level set function φ and the diffusion matrix A, see in particular the relation (2.15) below. In what follows we generalize ideas from [9, Section 2.1]. Consider for p ∈ Γ the parameter-dependent system of ODEs γ p (s) = A(p) ∇φ(γ p (s)) A(p) ∇φ(γ p (s)) · ∇φ(γ p (s)) , γ p (0) = p. (2.11) It is not difficult to see that there is δ > 0 such that the solution γ p of (2.11) exists uniquely on (−δ, δ) for every p ∈ Γ, so that we may define the mapping F : Γ × (−δ, δ) → R n+1 by F (p, s) := γ p (s). Recalling that a ij ∈ C 2 (Γ) we infer with the help of well-known results on the differentiability of solutions of ODEs with respect to parameters and initial conditions that F ∈ C 2 (Γ × (−δ, δ); R n+1 ). Furthermore, (2.11) implies that d ds φ(γ p (s)) = ∇φ(γ p (s)) · γ p (s) = 1, which yields φ(γ p (s)) = s for |s| < δ, p ∈ Γ, since φ(γ p (0)) = φ(p) = 0. In particular, F is a bijection from Γ × (−δ, δ) onto U δ with inverse where p ∈ C 2 (U δ ; R n+1 ) satisfies It is not difficult to verify that p(x) = p(x) in the case A = I and φ = d.
Using p we may define an alternative extension for a given u : Γ → R to U δ by setting (2.14) It is easily seen that p(γ p (s)) = p, p ∈ Γ, so that s → u e (γ p (s)) is constant on (−δ, δ). Differentiation with respect to s, together with (2.11), then implies that Suppose in addition that u is a solution of the surface PDE (2.3). It is shown in Lemma A.4 of the Appendix that u e then satisfies the uniformly elliptic PDE with A e (x) := A(p(x)), a e 0 (x) := a 0 (p(x)), f e (x) := f (p(x)) and

Phase field approach and finite element approximation
Let us next derive a suitable localized weak formulation of (2.16), which we shall use later in order to formulate our numerical scheme. Let σ ∈ C 0 (R) be such that σ(r) > 0, |r| <r and σ(r) = 0, |r| ≥r. A concrete choice of σ will be made later. For ε ∈ (0, δ r ) we define the phase field function The restriction on ε ensures that supp( ) = U εr ⊂ U δ . For a function v ∈ L 1 (Ω) we obtain with the help of the coarea formula The latter expression explains the scaling factor ε −1 and the weight |∇φ|, which will frequently occur.
Let us now multiply (2.16) by v |∇φ| with v ∈ H 1 (U r ) for some 0 < r < δ and integrate over U r . For the leading term we obtain with the help of integration by parts where we have used (2.15) to see that A e ∇u e · ∇ = 0. For the same reason the boundary term vanishes as the unit outer normal to ∂U r is a multiple of ∇φ. Thus, we obtain that We now use this relation in order to introduce our numerical scheme. To do so, let us assume for simplicity that Ω is polyhedral and denote by T h a regular partitioning of Ω into simplices T , i.e.
We denote by I h : C(Ω) → V h the Lagrange interpolation operator. Note for q > n+1 m , m = 1 or 2 and = 0 or 1 that In particular we infer from (2.10) that there exists an Next, let T be the unit simplex in R n+1 and a quadrature rule which is exact for all polynomials of degree ≤ q. This gives rise to a quadrature rule on T via T and Φ T is the usual affine transformation from T onto T . Using a standard application of the Bramble-Hilbert lemma we obtain for the quadrature error E T (g) : The degree of exactness of the quadrature formula now enters our choice of profile function σ, which we define as (2.26) A straightforward calculation shows that the corresponding phase field function ( .

(2.27)
In order to set up our numerical scheme we define for h < ε giving rise to the computational domain Lemma. 2.1. Denote by r 0 ∈ (0, 1) the unique zero of the function r → arccos(r) − c 1 r and set .10)). Suppose that ε = γh for some γ > 1 r 0 and that c 2 ε < δ. Then we have ε > 0 and We may assume in addition that x ∈ U ε π 2 as otherwise (x) = 0. Then we have which yields the desired estimate.
Next, let us define the finite element space Motivated by (2.20), our fully practical scheme reads: where the forms a h and l h are given by Furthermore, we have abbreviated A(x) := A( p(x)), a 0 (x) := a 0 ( p(x)), f (x) := f ( p(x)) and remark that these are used in the scheme rather than A e , a e 0 and f e , since in practice the evaluation of p is easier compared to p.
Remark. 2.1. In contrast to other methods, which require the determination of and integration over an approximate surface Γ h or a suitable narrow band, the implementation of (2.32) is rather straightforward. The underlying geometry is incorporated through the level set function φ and the projection p. Note that p is only required at the grid points of T h .

Let us introduce
In view of (2.24), (2.5) and (2.6) there exists c 3 > 0, which is independent of h, such that In particular we have: Proof. It is sufficient to verify that the homogeneous problem only has the trivial solution.
Let us formulate the main result of this paper.
Theorem. 2.1. Let u ∈ H 2 (Γ) be the unique solution of (2.3) extended to u e via (2.14) and u h ∈ V h the unique solution of (2.32). Let ε = γh for γ > 1 r 0 . Then The proof of these results will be given in the next section.
Remark. 2.2. The three terms on the right hand side of (2.37) are related to the different approximations that are used in the discretization. The first term is due to the use of piecewise linear finite elements in order to discretize the solution and the level set function, while the second term arises from working with the extended PDE in a narrow band of width ε = γh.
Here, the factor γ > 1 roughly measures how many grid points are used across the narrow band, whereas γ −(q+1) reflects how well integrals involving the phase field function are approximated via the quadrature rule.

Error Analysis
Before we start with the actual error analysis, we first prove a useful auxiliary result.
Proof. Let us fix T ∈ T h . Using (2.27) and Young's inequality we have for every Taking the maximum with respect to x and recalling (2.36) we infer that Combining this bound with (3.2) and observing that L i=1 ω i = 1 we obtain which concludes the proof of the lemma after summation over T ∈ T h .
Let us now start the proof of the error bound. Define e h := (I h u e ) |D h − u h ∈ V h . We infer from (2.35) and (2.32) Recalling the definition of a h we may write [A e ∇u e · ∇e h + a e 0 u e e h ] |∇φ| dx =: Using (2.25) and (2.24) we obtain Using similar arguments we deduce that as well as (3.6) Since A ∈ C 1 (Γ), it follows from (A.12) and (2.29) and, similarly, | a 0 (x) − a e 0 (x)| ≤ C ε 2 . This implies together with (3.1) and (A.3) Combining (3.4)-(3.7) we infer that Next, it follows from (2.33b) that Arguing in a similar way as for S 1,i , i = 1, 3, 4, 5, we obtain (3.9) and hence where [A e ∇u e · ∇e h + a e 0 u e e h − f e e h ] |∇φ| dx = I + II.

Numerical Experiments
We investigate the experimental order of convergence (eoc) for the following errors: The corresponding calculations will be done for a circle (Example 1) and a sphere (Example 2) of radius 1, described as the zero level set of the function φ(x) := |x| 2 − 1. In this case one can verify without difficulty that the projection p constructed in Section 2.2 coincides with the closest point projection p, so that we have u e (x) = u( x |x| ) for x = 0. We use the finite element toolbox Alberta 2.0, [29], and implement a similar mesh refinement strategy to that in [1] with a fine mesh constructed in D h and a coarser mesh in Ω\D h . The resulting linear systems were solved using CG together with diagonal preconditioning. In all the examples we consider we set a ij = δ ij , i, j = 1, . . . , n + 1 and a 0 = 1 in (2.3).

Example 1
Let Ω = (−1.2, 1.2) 2 and take Γ = {x ∈ R 2 | |x| = 1} to be a circle of radius 1, described as the zero level set of the function φ(x) := x 2 1 + x 2 2 − 1. In addition to E 1 , E 2 we shall also investigate the errors appearing in (2.38). To do so, we approximate u − u h respectively, where we have chosen the quadrature points In our computations L = 200 turned out to be sufficient. We choose f so that u(x) := (x 2 1 − x 2 2 )/|x| 2 solves (2.3) and fixed ε = 5.333h. In Table 1 we display the values of E i , i = 1, . . . , 4, together with the eocs, for q = 2, while in Table 2 we display E i , i = 1, . . . , 4, together with the eocs, for q = 6. For the smaller value q = 2 we observe an eoc for E 2 which is lower than two indicating that in this case the term γ −(q+1) in (2.37) dominates. This effect disappears for the choice q = 6, where we see eocs close to two for E 2 and E 4 . Furthermore, we observe eocs close to four for E 1 and E 3 suggesting that the error analysis can be improved for the L 2 -errors.

Example 2
We set Ω = (−1.8, 1.8) 3 and take Γ = {x ∈ R 3 | |x| = 1} to be a sphere of radius 1, described as the zero level set of the function φ(x) := x 2 1 + x 2 2 + x 2 3 − 1. As in Example 1, in addition to E 1 , E 2 we shall also investigate the errors appearing in (2.38) which we approximate by the quadrature rules and Here,  Tables 3 and 4 we display the values of E i , i = 1, . . . , 4, together with the eocs, for q = 1 and q = 6 respectively and observe a similar behaviour as in the two-dimensional test example.
In Figure 1 we display the approximate solution u h plotted on the zero level surface of I h φ.

Results using piecewise quadratic finite elements
Even though we have restricted our error analysis to the case of piecewise linear finite elements it is not difficult to apply our approach to quadratic elements. In order to do so, we use 7.500e-02 0.4 1.134e-06 -8.248e-04 -1.439e-06 -2.079e-03 -3.750e-02 0.2 3.627e-08 4.97 1.440e-04 2.52 9.382e-08 3.94 5.034e-04 2.05 1.875e-02 0.1 1.721e-09 4.40 3.212e-05 2.16 6.245e-09 3.91 1.308e-04 1.94 9.375e-03 0.05 9.899e-11 4.12 7.789e-06 2.04 3.820e-10 4.03 3.197e-05 2.03 Table 4: Errors and experimental orders of convergence, q = 6 instead of (2.31) and define the forms a h and l h (for the case a ij = δ ij , a 0 = 1) by where I h denotes the Lagrange interpolation operator for piecewise quadratic finite elements. The results in Table 5 correspond to the setting outlined in Example 1. Using a quadrature rule of order q = 6 we see eocs close to order four for E 2 and E 4 in contrast to the eocs close to order two, that are displayed in Table 2, for the corresponding affine finite element approximation. The fact that the eocs for E 1 and E 3 are close to four (rather than six as expected for quadratic elements) is a consequence of the term ε 2 = γ 2 h 2 in (2.37) which now dominates.

A Appendix
The aim of this appendix is to derive certain properties of the projection p and the extension u e (x) = u(p(x)) which have been used in the analysis above. To begin, we infer from the definition of u e for 1 ≤ i, j ≤ n + 1 that Then Proof. Using the transformation F : Γ × (−r, r) → U r with Jacobian determinant µ, (A.1), (A.2) and the fact that p ∈ C 2 we obtain and the result follows.
Proof. Recalling that A(p)ν(p) = ν(p), p ∈ Γ as well as η x (0) = x, η x (1) = p(x) we obtain with the help of (A.4) Note that d A k ∈ C 2 , since this is true for x → η x and x → a kl (p(x)). The relation (A.6) immediately follows from (A.5). Next, observing that ∇ Γ f (p(x)) ∈ T p(x) Γ and ∇φ(p(x)) ⊥ T p(x) Γ we infer that which implies (A.7) in a similar way as above.
Taylor's theorem together with (A.9) and (A.10) implies the existence of d p,0 k ∈ C 2 with The relation (A.11) allows us to prove a bound between p(x) and the closest-point projection p(x), which is used in the error analysis.