Higherorder meshing of implicit geometries—part II:
Approximations
on manifolds
Abstract
A new concept for the higherorder accurate approximation of partial differential equations on manifolds is proposed where a surface mesh composed by higherorder elements is automatically generated based on levelset data. Thereby, it enables a completely automatic workflow from the geometric description to the numerical analysis without any userintervention. A master levelset function defines the shape of the manifold through its zeroisosurface which is then restricted to a finite domain by additional levelset functions. It is ensured that the surface elements are sufficiently continuous and shape regular which is achieved by manipulating the background mesh. The numerical results show that optimal convergence rates are obtained with a moderate increase in the condition number compared to handcrafted surface meshes.
Keywords: higherorder FEM, manifold, surface PDEs, levelset method
0] ] ] ] ] ] ]
Institute of Structural Analysis
Graz University of Technology
Lessingstr. 25/II, 8010 Graz, Austria
www.ifb.tugraz.at
fries@tugraz.at
Contents
1 Introduction
Many challenging applications in engineering and natural sciences are characterized by physical phenomena taking place on curved surfaces in the threedimensional space. There are numerous examples for transport and flow phenomena on biomembranes or bubble surfaces [26, 39]. Examples in structures are membranes and shells [7, 3]. Phenomena on surfaces may also be coupled to processes in the surrounding volume such as in surfactant transport, hydraulic fracturing, reinforced structures etc. As an additional challenge, the surfaces may be moving [18, 40, 22], i.e. the domain of interest changes. The modeling of such phenomena naturally leads to boundary value problems where partial differential equations are formulated on manifolds. For the solution of such models, customized numerical methods are needed.
The first application of the finite element method for the solution of the LaplaceBeltrami operator on manifolds is reported in 1988 by Dziuk [17]. Since then, the topic has attracted a tremendous research interest leading to a variety of numerical methods for PDEs on surfaces existing today, see [20] for an overview. The most straightforward approach is to generate surface meshes on the manifold and extend the finite element method in a natural way using tangential differential calculus. That is, standard gradients of the planar twodimensional case are replaced by surface gradients on the manifold. It is interesting to note that this approach has been chosen from the beginning in the simulation of transport phenomena, e.g. [2, 1, 16, 20]. However, for the modeling of membranes and shells, a less intuitive path using local coordinate systems and Christoffel symbols is standard since a long time [7, 3]. It is rather recent that these models have been recasted in the frame of global tangential operators [10, 11, 28, 29].
Another approach is to only employ an implicit description of the manifold and solve the model equations on all isosurfaces at once [19, 4]. Then, the problem is naturally set up in the threedimensional space embedding the manifolds, i.e. volumetric elements and shape functions are employed. However, typically only the solution on one isosurface, say the zeroisosurface, is of interest. One may then restrict the surrounding domain to a narrow band around the manifold [2, 8, 21]. There are interesting similarities to phase field and diffuse interface approaches [37]. A recent approach is to collaps the narrow band to the manifold itself. Then, shape functions of the volumetric background elements are used, however, the integration takes place on the trace of the manifold only [32, 31, 28]. The resulting approaches are labeled TraceFEM [9, 25, 32, 35, 38] or CutFEM [27, 28]. Higherorder approximations of PDEs on manifolds have been reported in different contexts before: For explicit handcrafted surface meshes in [12] and in the context of the TraceFEM in [38]. Adaptivity is considered e.g. in [13, 14].
Herein, we propose a higherorder accurate approach for the approximation of PDEs on manifolds. The manifold is described implicitly based on the levelset method. A surface mesh composed by mixed higherorder quadrilateral and triangular finite elements is automatically generated from a background mesh and given levelset functions. A master levelset function defines the shape of the manifold. However, as the implied zeroisosurface may be infinite, it is restricted by additional (slave) levelset functions. That is, several levelset functions imply the bounded manifold being the domain of interest in the BVP. As a model problem, we consider the LaplaceBeltrami operator and an instationary advectiondiffusion problem. Based on this, the extension of the approach to more advanced transport problems on surfaces and in the simulation of membranes and shells will be reported in forthcoming publications.
The automatic detection of higherorder surface elements has been reported by the authors in [23, 24] in the context of integration and interpolation. There, only a set of surface elements is needed featuring double nodes and not necessarily fulfilling continuity. In order to be suited for the approximation of PDEs on surfaces as discussed herein, (1) continuity requirements have to be fulfilled, (2) the elements must be sufficiently shape regular, and (3) connectivity information in the usual FEM sense has to be provided, enabling the concept of nodal degrees of freedom. These issues are addressed herein with emphasis on higherorder accurate approximations. Also, the concept of using several levelset functions for the definition of the bounded manifold is new and an extension of [24].
The paper is organized as follows: In Section 2 we outline the geometric description of the bounded manifold based on several levelset functions defined on a background mesh composed by higherorder elements. The automatic generation of suitable higherorder surface meshes is described in Section 3: The reconstruction of surface elements approximating the zeroisosurface of the master levelset function, the restriction by means of additional (slave) levelset functions, the extraction of a continuous surface mesh from the element set, and the manipulation of the background mesh to achieve shape regular elements. Section 4 shortly recalls the standard finite element approach for approximations on meshed surfaces. Numerical results are presented in Section 5 for curved lines in two dimensions and curved surfaces in three dimensions. The LaplaceBeltrami operator is considered as well as instationary advectiondiffusion on manifolds. Finally, a summary and outlook is given in Section 6.
2 Preliminaries
The task is to solve a boundary value problem (BVP) on a surface in three dimensions. Let the surface be possibly curved, sufficiently smooth, orientable, connected (so there is only one surface), and feature a finite fixed area. The surfaces discussed herein are defined following the concept of the levelset method. That is, there is a continuous scalar function with and its zerolevel set
(2.1) 
defines a (zero)isosurface. We note that a general isosurface of may be unconnected, see Fig. 1(a). Then, it is necessary to select the one zerolevel set of interest for the solution of a BVP, so this is not a problem. Furthermore, isosurfaces of levelset functions are orientable by default, i.e. there is a consistent normal vector pointing from the negative to the positive side, see Fig. 1(b).
There is no concept of a boundary when only one levelset function is considered with . As a consequence, isosurfaces with finite area must be closed (compact) as e.g. shown in Fig. 1(c), hence, there is no boundary. Otherwise, in the more general case, they are open and feature an infinite area without a boundary. It is thus seen that the manifold of interest, i.e. the surface where a BVP is to be solved, is typically only a subset of , hence, .
There are two alternatives to define bounded isosurfaces with finite area, see Fig. 2(a), which fully represent the manifold of interest: One is to evaluate the levelset function only in a subregion , hence
(2.2) 
See Fig. 2(b) for an example, where is only evaluated in the gray subregion instead of . The other is to employ additional levelset functions to restrict , hence
(2.3) 
See Fig. 2(c) for an example, where two additonal isosurfaces of and are shown in blue. The same bounded isosurface from Fig. 2(a) is thereby defined.
Further examples are shown in Fig. 3. In Fig. 3(a) to (d), zeroisosurfaces are shown, in yellow of and in blue of . It is seen that complex, bounded manifolds are thereby implied, see Fig. 3(e) to (h). It is clear, that it would be very difficult to define such complex manifolds by using only and defining subregions . Therefore, we prefer the second alternative of using multiple levelset functions which are evaluated in . Of course, a combination of these two alternatives is straightforward: This is in particular useful when the outer boundary of the manifold is straight and follows from a boxlike subdomain, e.g. , but features internal boundaries which are rather described by .
An important consequence of using several levelset functions to restrict is that corners in the boundaries of the manifold may be defined straightforwardly (although each is smooth). This seen in Fig. 3(f) to (h) where more than one is present. It is noted that the employed levelset functions and do not have to be signed distance functions.
3 Mesh generation of the manifold
The procedure is outlined as follows: First, a background mesh is introduced such that the manifold of interest is completely immersed. The levelset data is evaluated at the nodes of the higherorder elements of this background mesh. Second, the isosurface is identified and meshed which is called reconstruction. As mentioned before, this isosurface is either closed or is bounded by the boundary of the background mesh. Third, the meshed isosurface is restricted by the additonal levelset functions .
3.1 Background mesh
A background mesh is introduced such that the manifold of interest is completely immersed. The mesh is composed by higherorder background elements of the Lagrange class and we restrict ourselves to tetrahedral elements for simplicity. Hexahedral elements could be decomposed into tetrahedra and treated similarly, see [23]. There are no restrictions on the background mesh, i.e. it may be unstructured and the element faces may be curved.
We shall later discuss the issue of manipulating the background mesh by moving nodes in order to improve the resulting surface elements of the manifold mesh. We find that the concept of universal meshes [36], extended to three dimensions as shown in Fig. 5(b), allows for a maximum of flexibility in the node movements. But natural choices are also the blocktype meshes shown in Fig. 5(d). This is in particular useful when the manifold of interest features straight edges and follows from Eq. (2.2) with a blocktype subdomain, e.g. .
3.2 Reconstruction of
In each element of the background mesh which is cut by the zeroisosurface of , the task is to identify element nodes of a higherorder surface element on the isosurface. This has been described in detail in the first part of this series of publications, see [24] and also [23]. The procedure is shortly summarized as: In each tetrahedral background element, it is confirmed that valid levelset data is present which implies only two possible cut scenarios of the isosurface: The isosurface is either of triangular shape, see Fig. 6(a), then three edges are cut. Otherwise, four edges are cut and it is of quadrilateral shape, see Fig. 6(b). Depending on the topology, customized search algorithms are evaluated in the reference background element.
First, element nodes on the faces of the tetrahedron are identified leading to the edge nodes of the sought surface element approximating the isosurface. This is, in fact, realized by reconstructing a higherorder line element in a twodimensional, triangular reference element with the levelset data from the respective face nodes of the tetrahedron. This line is mapped to the corresponding face of the tetrahedron. The inner nodes of the sought surface element are then determined in the reference tetrahedron wherefore the search algorithm also makes use of the position of the edge nodes determined beforehand. Finally, the determined surface element may be mapped to the physical background element using the isoparametric concept. An example is seen in Fig. 6(c) and (d) for the two different topological cases of a triangular or quadrilateral surface element, respectively.
Obviously, applying the described procedure in all cut background elements leads to a set of higherorder surface elements representing an accurate approximation of . Lateron, a mesh suitable for a finite element analysis of the BVP on the manifold will be automatically generated from this element set, see Section 3.4. An important property of such a mesh is that it is at least continuous so that there are no gaps in the representation of the manifold. It is interesting to note that the reconstruction described above does, in fact, not assure this property without further considerations.
The reason for this is found in the orientation of the face coordinate systems, see Fig. 7: Only when the same coordinate system is used for both tetrahedra sharing a face, this will lead to exactly the same edge nodes of the surface elements. Otherwise, only the corner nodes of the surface elements perfectly match. Fortunately, there is a simple solution to the problem: The face nodes of a tetrahedron which are extracted for the reconstruction of the edge nodes of the sought surface element should be permuted such that they imply the same coordinate systems on the shared face. The resulting set of surface elements is then, as desired, a globally continuous, higherorder approximation of . This is an important difference to the reconstructions discussed in [23, 24].
3.3 Restriction of by additional levelset functions
The focus is on one background element cut by the zerolevel set of . A triangular or quadrilateral surface element, approximating the zero isosurface, has been succesfully reconstructed. Assume that the background element and the surface element are also cut by another level set function , see Fig. 8(a) to (c) for different examples. A decomposition of the surface element is then sought into surface elements on the two sides, see Fig. 8(d) to (f) for the same examples, respectively.
For this decomposition, it is necessary to first interpolate the nodal values of in the background element (crosses) at the nodes of the reconstructed surface element (dots). Then, the decomposition is carried out in twodimensional reference elements with these nodal values, see Fig. 8(g) to (i). As discussed in detail in [23, 24], different cut scenarios may be faced: In triangles, always one of the corners is on the other side than the other two and one subtriangle and one subquadrilateral are obtained, see Fig. 8(g). In quadrilaterals, either two neighboring nodes are on the other side leading to two subquadrilaterals, see Fig. 8(h), or only one node, leading to four triangles, see Fig. 8(i). The decomposed twodimensional elements are mapped to the cut tetrahedral element using an isoparametric map implied by the original surface element (from the reconstruction of ). Typically, considering Eq. (2.3), only the subelements on the negative side of are needed and the others are neglected, this may be called restriction.
The procedure is very similar when several levelset functions cut the same element, possibly implying one or more corners in the boundary of the manifold. The starting point is then a set of surface elements resulting from a previous restriction with respect to other levelset functions. Then, at each surface element, the procedure from above is repeated: The nodal values of the current levelset function at the current surface element are interpolated from the tetrahedral nodes. Then, the decomposition is carried out in the twodimensional reference elements and mapped to the current surface element. Finally, the resulting subelements are restricted to those with . See Fig. 9 for some examples. It is emphasized that the reconstructions are always realized in reference elements (before mapping them to the situation in global coordinates): For the levelset function , these are 3D reference tetrahedra and for the levelset functions , these are reference triangular or quadrilateral elements. This enables a very efficient and robust implementation as the reference elements always feature straight edges and planar faces, respectively. On the contrary, it seems impossible to reconstruct directly in general, possibly curved, physical elements.
3.4 Mesh generation
Starting point is a set of higherorder surface elements which has been reconstructed from and restricted to the manifold of interest using the additional levelset functions . In order to approximate a BVP on the manifold, a finite element mesh has to be generated from the element set. The aim is to extract a set of unique nodes, introduce a global node numbering and set up a connectivity matrix defining the elements. It is noted that the generated element set from above contains triangular and quadrilateral elements, which leads to a mixed finite element mesh. Of course, the decompositon of quadrilateral elements to triangular elements (or vice versa) is possible if meshes of one element type are preferred.
There are two obvious alternatives for the automatic mesh generation from the element set: In alternative 1, the nodes are first sorted with respect to their coordinates and then repeated nodes, such as they appear on the element edges and corners, are unified and associated with a node number. This may be implemented very efficiently, however, it requires some userdefined threshold characterizing a distance within which nodes are considered equal. Because the overall set of elements may contain extremely small elements (e.g. when the zeroisosurface cuts very close to a corner node of a background element) this may interfere with the threshold leading to invalid data. For example, all element nodes may be unified to one node in an extreme situation.
Therefore, we prefer alternative 2 where the connectivity information of the background mesh is employed. One may easily determine the neigboring elements of each tetrahedron. For the surface elements in the reconstructed set of elements, the information of the corresponding background elements is stored, respectively. One may then directly compute which nodes on a shared face of two tetrahedra must match. The same holds for the nodes on a shared edge of several tetrahedra. No distance computations are needed in this approach. It is clear that an advanced implementation is possible where the computation of nodes with the same coordinates (at the edges or corners of the surface elements belonging to different background elements) is completely avoided.
3.5 Manipulation of the background mesh
Some example meshes for complex manifolds are seen in Fig. 4. It is obvious that awkward element shapes result from the automatic meshing. Very small inner angles may occur and the size of neighboring elements can vary strongly so that the shape regularity is an issue. This may lead to the assumption that the elements are not suitable for the analysis, even less in a higherorder context. However, as shown in [31, 33] in a loworder context, only a few manipulations ensure that the “maximal angle condition” is fulfilled and the spectral condition number is bounded uniformly. More generally, it has recently been shown in several publications by the groups around Hansbo, Burman and Reusken, e.g. [5, 33, 38], that optimal results for BVPs on implicitly defined manifolds are possible with such meshes. They use finite element shape functions based on the background mesh which are only evaluated on the manifold; these approaches are labelled TraceFEM or CutFEM. They introduce tailored stabilization terms which avoid negative effects resulting from the illshaped elements. We shall report on the use of such methods in this context in a future publication. Herein, we wish to avoid the detailed (mathematical) discussion of these approaches.
A simple and intuitive attempt is to manipulate the background mesh such that the resulting reconstructed elements remain shape regular. In particular, the area of the surface elements shall be bounded from below. When the condition numbers in a finite element analysis obtained by a tailored, manufactured surface mesh, ideally featuring the same element shapes and areas in the whole mesh, are compared to those of the automatically generated meshes from above, it is seen that there is typically an increase by the factor of the largest to the smallest element area . Because typically scales with the size of the (very regular) background elements, this clearly motivates the need to bound from below.
For simplicity, we focus on the set of elements resulting from a reconstruction with respect to , i.e. the set of elements approximating on a finite background mesh. An algorithm is proposed which moves the nodes of the background mesh away from the zerolevel set of . Basically, this ensures that the zeroisosurface is not too close to the corner nodes of the background mesh and, consequently, that the surface elements are not too small. Only the nodes in a close band around the zero isosurface are moved.
Although the elements of the background mesh may be arbitrarily curved, it is quite useful to enforce straight edges for simplicity of the node manipulation. The algorithm for the node movement is described for one concrete node at position in the background mesh; the same holds for all other nodes. It may be seen as a fixpoint iteration and the following steps are repeated until all nodes have been moved sufficiently away from the zerolevel set.

Approximate the distance of the node at to the zerolevel set. When the levelset function is given in analytical form, this is realized by a NewtonRaphson procedure,
using as the starting point. This algorithm converges to some position on the zerolevel set which is not necessarily the closest point. However, it is found that is approaching the closest point quickly the closer is to the zerolevel set. Because the node manipulation is only realized for nodes which are close to the zerolevel set, it turned out to be sufficient to continue with the computed position . The estimated signed distance is
and the direction
is computed as well; there holds . In the case that is only given at the nodes of a higherorder background mesh (and no analytic information is available), one may estimate the distance through a linear reconstruction of the zero isosurface which is much cheaper than the desired higherorder reconstruction performed after the node manipulation.

Move the node provided that it is sufficiently close to the zerolevel set, i.e. when . scales with the resolution of the background mesh and we typically set with being some characteristic element length. The node is moved as
with
That is, the maximum distance a node is moved in one step is which is often set as . Obviously, scales linearly with the distance of a node from the zerolevel set.

Set and repeat the two steps from above until where, typically, .
This algorithm works very fast and scales linearly with the number of nodes in the background mesh. It is controlled by the three parameters which are summarized as: and define local regions around the zerolevel sets (“bands”): Nodes within the band controlled by are considered for movement and defines the band within which no nodes shall be present. defines the maximum moving length in one iteration which is scaled linearly in the band controlled by . Obviously, we need .
It is mentioned that in our first approximations of BVPs on manifolds with the automatically generated surface meshes, we were surprised by the accuracy of the results even without any mesh manipulations. Actually, it seems that as long as direct solvers are employed (in our case Matlab’s backslash solver), there is a large range of condition numbers still leading to accurate results even in the frame of systematic convergence studies. Nevertheless, the node manipulations are used by default in the numerical studies presented in Section 5. It is an important advantage that they also render recursive refinements of the background elements unnecessary as long as the resolution of the background elements is reasonably adjusted to the complexity of the zerolevel sets. The effect of the proposed node manipulation is now demonstrated in one, two an three dimensions.
3.5.1 Node manipulation in 1D
Consider the onedimensional example shown in Fig. 10. In Fig. 10(a), a levelset function is plotted in black and the corresponding signeddistance function , which stores the (signed) distance to the closest root of , is shown in blue. On the axis, there are regularly distributed nodes marked as blacked dots, the distance between two nodes is . The aim is to move these nodes away from the roots (blue circles). In Fig. 10(b), we set and . It is then seen that after about 10 steps of the fixpoint iteration, the nodes in the region are squeezed to the boundary between the band scaled by and the rest of the domain. The resulting nodal distribution leads to undesired, locally concentrated background elements. Therefore, we set and next. Results with and are shown in Figs. 10(c) and (d), respectively. Obviously, the larger the value for , the faster the nodes are moved out of the band controlled by . Nevertheless, in more than one dimension, setting considerably smaller than allows for curved node paths throughout the nodal movements leading to smoother element distributions in the background meshes.
3.5.2 Node manipulation in 2D and 3D
Examples for manipulations of background meshes in 2D are shown in Fig. 11. Two different levelset functions on a coarse and a fine mesh are considered, respectively. The red lines indicate the zerolevel sets which are meshed by higherorder line elements. It is seen that due to the nodal movements, the lengths of the line elements in the background mesh are quite regular, hence, the ratio is bounded. The results are achieved with , , and .
A similar study is repeated in 3D. In Fig. 12, the resulting surface meshes are shown for two different levelset functions on different background meshes, respectively. The manipulated background meshes are omitted for clarity and the focus is on the shape and size of the surface elements. The element areas of the largest and smallest elements are computed and given in Tables 1 and 2. As can be seen, the ratio is well bounded; it is typically less than . The increase in the condition number when using the automatically generated meshes as proposed herein (compared to manually constructed meshes) is in the same range. This is not a problem as long as direct solvers are employed, however, will increase the effort in the context of iterative solvers.
coarse mesh  

medium mesh  
fine mesh 
coarse mesh  

medium mesh  
fine mesh 
4 PDEs on manifolds
Herein, we consider the stationary Poisson equation and the instationary advectiondiffusion equation on curved manifolds. Both involve the LaplaceBeltrami operator. The PDEs are conveniently defined based on the tangential differential calculus, i.e. on surface gradients (sometimes refered to as intrinsic gradients). The manifold is implicitly defined by the zerolevel set of a sufficiently smooth levelset function. The tangential gradient operator of a scalar function on the manifold is defined as
is the discrete tangential projector, is the unit normal vektor on the manifold and is a smooth extension in a neighborhood of the manifold . With the tangential divergence of a vector , , the LaplaceBeltrami operator is [31, 20]
(4.1) 
Assume a sucessfull reconstruction of the manifold , the discrete surface is the union of the elements , where is the set of automatically generated Lagrange elements as described above. The standard surface FEM as described in [12, 20] is used to solve the PDEs on the discrete manifold. In the context of the FEM the function is often a shape function living in a reference element with coordinates . Using the standard isoparametric mapping from the reference element to the physical element, , and being nodal coordinates, the finite element space on with the polynomial degree is then defined by
(4.2) 
is the polynomial basis in a triangular element and in a quadrilateral element, both being complete of order . This space is spanned by the nodal basis , with being the total number of nodes, which is given by
A base function is defined by . Note that the function space is a subspace of the continuous space . For more information we refer to [12, 20].
In general, it is desirable to compute the surface gradient without explicitly computing an extension . Assume that the surface element in the physical space follows from the standard isoparametric mapping , then the surface gradient is directly obtained from
(4.3) 
with being the ()Jacobi matrix and being the first fundamental form.
4.1 Stationary Poisson equation on manifolds
The complete boundary value problem in strong form for the Poisson equation on manifolds states that we seek such that:
(4.4a)  
(4.4b)  
(4.4c) 
where is a source function and is the conormal vector being normal to the boundary yet in the tangent plane of the manifold at each point on the boundary.
For a closed (compact) manifold, where no boundary exists, one needs an additional condition for the problem to be wellposed. Therefore, typically the zero mean constraint is imposed,
(4.5) 
The boundary value problem from above is converted to a weak form and discretized using the surface FEM. Therefore, we introduce the following trial and test function spaces
(4.6)  
(4.7)  
(4.8) 
Then, the discrete weak form for the Poisson equation on a manifold is stated as: Find such that
(4.9) 
The existence and uniqueness of this BVP is shown in [20, 6]. The zero mean constraint for compact manifolds is enforced by a Lagrange multiplier.
4.2 Instationary advectiondiffusion equation on manifolds
The strong form of the instationary advectiondiffusion equation in the time interval from to reads
on  (4.10a)  
on  (4.10b)  
on  (4.10c)  
on  (4.10d) 
with , being the advection velocity in the tangent space of the manifold and the diffusion coefficient. Using the same function spaces from above, the semidiscrete weak form becomes: Find with such that
(4.11) 
The zero mean constraint from (4.8) for compact manifolds is not needed here. The existence and uniqueness of the solution is guaranteed by the LaxMilgramm lemma and is proven for the stationary problem in [34].
After the FEM discretization in space, a firstorder system of ordinary differential equations in time is obtained as usual. We use a thorder accurate, implicit, step RungeKutta method based on GaussLegendre points to discretize in time. The corresponding Butcher tableau for this RungeKutta scheme is given as
(4.12) 
It is useful to employ this highly accurate scheme in time to virtually eliminate the temporal error in the numerical studies and focus on the spatial approximation error of the surface finite elements (with different orders).
5 Numerical results
As defined in the previous section, the stationary Poisson equation and the instationary advectiondiffusion equation are considered here. The focus above is on curved surfaces in three domensions, however, results are also shown for curved lines in two dimensions. Convergence studies for different element orders up to are conducted where the analytic solutions are generated by the method of manufactured solutions. In each test case, we compare the performance of handcrafted meshes composed by very regular elements with the automatically generated meshes as described in Section 3. The manipulation of the background mesh from Section 3.5 is enabled in all examples.
According to the method of manufactured solutions in the context of the Poisson equation, the source function on the right hand side is computed such that Eq. (4.4a) is fulfilled for a chosen solution . That is, the LaplaceBeltrami operator (in local coordinates)
(5.1) 
is applied to to achieve . There, with entries is the first fundamental form and with entries its inverse.
For the convergence studies, the error is shown over the element lengths of the corresponding meshes. For the handcrafted meshes, is the element length of a typical element of the surface mesh (where elements are quite uniformly shaped). For an automatically reconstructed surface mesh, the element lengths vary largely so that is the characteristic element length of the background mesh (from which the surface mesh was reconstructed).
5.1 Poisson equation on curved lines in 2D
5.1.1 Circular manifold
The Poisson equation is first solved on a circle with radius . The circle is defined by the implicit levelset function . Using polar coordinates , the solution which satisfies the zeromean condition from Eq. (4.5) is chosen as [30]
(5.2) 
with the corresponding right hand side resulting as
(5.3) 
The element sizes of the background meshes are systematically varied as for the convergence studies. One example of a coarse background mesh and a reconstruction by quadratic elements is shown in Fig. 13(a). The corresponding numerical solution is shown in Fig. 13(b) together with the analytical solution. Fig. 13(c) shows the ratio of the largest element length of the line mesh over the smallest for the different background meshes. These numbers are very similar for the different element orders, so that we only show one curve being representative for all orders. It is seen that the worst ratio is less than in this example which shows the success of the manipulation of the background mesh. Of course, for handcrafted meshes on the circle, it is not a problem to achieve the optimal ratio of in this test case.
Fig. 14 shows the convergence rates and condition numbers. Fig. 14(a) and (c) are achieved on handcrafted meshes with equidistant line elements. As expected, optimal convergence rates are seen in the norm and the condition number scales with . Fig. 14(b) and (d) show the corresponding results on the automatically generated meshes. The convergence rates are again optimal (the dotted lines show the optimal slope) and the condition numbers also scale with . Obviously, compared to the handcrafted meshes, the error level is modestly larger and the condition number is larger by about two order of magnitude.
5.1.2 Flowershaped manifold
Next, the LaplaceBeltrami operator is solved on a closed, flowershaped manifold implied by the levelset function
with and . The function is defined as and the right hand side is computed accordingly. See Fig. 15(a) for a coarse background mesh and a reconstructed line mesh of order . The exact solution is seen in Fig. 15(b) and the largest ratio of the reconstructed line elements resulting from the different background meshes in Fig. 15(c).
In the convergence studies, we use background meshes with which is fine enough to ensure that the reconstruction is successfully achieved without any (recursive/adaptive) refinements of the background elements. For the handcrafted meshes, significantly larger element lengths may be used. The results presented in Fig. 16 follow the same style than in the previous test case. It is confirmed that optimal convergence rates are achieved on handcrafted and automatically reconstructed meshes and the condition numbers behave as expected.
5.1.3 Sshaped manifold
In this example, we consider an open manifold where the boundary is defined by additional levelset functions . The manifold in parametric form is given as
(5.4) 
The corresponding levelset function may be defined as
(5.5) 
The boundaries are defined by the two additional levelset functions and . An example for a reconstructed line mesh is seen in Fig. 17(a) where also the zeroisolines of all three levelset functions are shown. The exact solution is plotted in Fig. 17(b). Again, the ratio of the largest element length of the automatically reconstructed line mesh to the smallest is shown in Fig. 17(c). This ratio may be worse than in the previous examples due to the presence of several levelset functions which renders the manipulation of the background mesh more complicated.
In Fig. 18, the results of the convergence analysis is presented in the same style than above with the same conclusions to be drawn. Hence, we conclude that approximations of the Poisson equation with optimal convergence rates are possible for the automatically generated meshes discretizing closed and open curved lines in 2D.
5.2 Poisson equation on curved surfaces in 3D
5.2.1 Quarter Cylinder
Here, the LaplaceBeltrami operator is solved on the surface of a quarter cylinder with the radius and height . The analytical solution in cylindrical coordinates is given in [30] based on the following parameters and functions:
(5.6)  
(5.7)  
(5.8)  
(5.9) 
The analytic solution is then defined as . Applying the LaplaceBeltrami operator yields the source function
The Dirichlet boundary condition is . In this special case, it is particularly simple to use a tensorproduct background mesh in the volume to reconstruct the higherorder surface mesh with the desired boundaries. Hence, the use of addional levelset functions would unnecessarily complicate the situation here.
An example of a reconstructed, higherorder surface mesh is shown in Fig. 19(a). The exact solution is given in Fig. 19(b). The ratio of the largest surface element to the smallest is given in Fig. 19(c) and is, again, clearly bounded (in the range of for this example).
For the convergence studies we use background meshes with and results are shown in Fig. 18. As for the previous convergence studies for curved lines, Fig. 18(a) and (c) refer to convergence rates and condition numbers of handcrafted surface meshes, respectively. Fig. 18(b) and (d) give the results of the automatically generated surface meshes. The resulting convergence plots show optimal convergence rates which are surprisingly smooth given the irregularity of the reconstructed surface elements. The condition numbers are about two orders of magnitude larger than in the handcrafted meshes.
5.2.2 Sphere
Next, we consider a sphere with radius . Using spherical coordinates , an exact solution fulfilling the zeromean constraint from Eq. (4.5) is chosen as
(5.10) 
Applying the LaplaceBeltrami operator yields the source function
(5.11) 
The test case is sketched in Fig. 21 showing (a) a reconstructed example mesh and (b) the exact solution on the sphere. The area ratios are given in Fig. 21(c). In Fig. 22, convergence results and condition numbers are only shown for the automatically reconstructed surface meshes because the findings coincide with the previous test case. It is wellknown that the inner angles in the elements play an important role for the quality of a mesh. In particular, large inner angles may hinder optimal convergence [15]. Therefore, for all meshes the maximum inner angles have also been computed. The generated meshes are typically mixed, i.e., composed by quadrilateral and triangular elements and the maximum angles have been computed for these two element types individually. Results are shown in Fig. 23 for first order and thorder meshes only because all other orders yield very similar results. As can be seen the maximum inner angles in quadrilaterals are below and in triangles below . This was also confirmed in the other test cases where one levelset function defines the manifold.
5.2.3 Hyperbolic paraboloid with bumps
Now we consider the LaplaceBeltrami operator on a hyperbolic paraboloid with bumps. The surface is defined by the function
with . The analytic solution is chosen as
(5.12) 
and the source function determined accordingly. The levelset function implying the manifold is
(5.13) 
If the reconstruction is achieved based on a cubelike background mesh with suitable dimensions, the manipulation of the background mesh is restricted in that nodes on the outer contour of the background mesh cannot be moved freely without changing the boundary of the implied manifold. Therefore, we use four additional levelset functions which restrict the manifold properly. The mesh manipulation of the (universal) background mesh, see Section 2, is then applied with respect to all levelset functions and the element ratios are sufficently bounded.
The additional levelset functions are
where is a scalar product and
The resulting zerolevel sets are planes restricting the zerolevel set implied by Eq. (5.13). See Fig. 24(a) for a visualization of the zerolevel sets and (b) the resulting bounded manifold. Fig. 24(c) and (d) show the exact solution and the area ratios of the resulting surface elements, respectively. Convergence results and condition numbers for the automatically reconstructed meshes are seen in Fig. 25.
5.3 Instationary transport equation in 2D and 3D
The instationary advectiondiffusion equation on manifolds is considered here. The semidiscrete weak form is given in Eq. (4.11). The time discretization is carried out using a thorder accurate, implicit, step RungeKutta method as defined in Eq. (4.12). We use a constant number of 4096 time steps in all upcoming examples which virtually eliminates the discretization error in time. First, we consider pure advection problems on a circle and a sphere, i.e. without any diffusion. Then, two test cases are considered with diffusion where no analytic solution is available. In the error studies, the error is measured in space at the final time step. Of course, this error is effected by the discretizations in space and time, wherefore we employ the highly accurate time stepping scheme mentioned above.
5.3.1 Pure advection on a circular manifold
The instationary transport problem is solved on a circle with radius in the time interval . The advection velocity is in tangential direction of the circle and the diffusion coefficient . The initial condition is
(5.14) 
This initial distribution is simply advected around the circle without changing its shape. That is, the exact solution is (5.14) with a shift in the angle of . See Fig. 26 for a graphical representation of the solution at different instances in time at .
For this study, the same meshes from Section 5.1.1 are used. Results are presented in Fig. 27 in the same style than above. It is seen that higherorder convergence rates are achieved. However, for this pure advection problem, one order of the optimal convergence rate is lost for even orders of the elements, i.e. This is the same for handcrafted meshes as well as for the automatically reconstructed meshes. We have confirmed that the same behavior also occurs for planar advection problems so this has nothing to do with the fact that the problem is solved on a curved manifold nor that no boundaries are present.
5.3.2 Pure advection on a sphere
A similar problem without diffusion is now solved on a sphere with radius in the time interval . The advection velocity is
taking place only in the plane tangential to the sphere. The initial state is set to
(5.15) 
Analogously to the example above, the exact solution is obtained by the rotation angle in the plane. In Fig. 28, the analytical solution at is illustrated.
The same meshes from Section 5.2.2 are used and results are seen in Fig. 29. The same conclusions from above may be drawn, in particular the loss of one order in the convergence rate for elements with even orders. It is recalled that this behaviour is typical for pure advection and optimal convergence rates are recovered in the presence of diffusion. This is equivalent to results of planar problems with handcrafted meshes.
5.3.3 Advectiondiffusion on an Sshaped manifold
Next, the instationary transport problem is solved on the Sshaped manifold of Section 5.1.3 in the time interval . The advection velocity is and the diffusion coefficient . A Dirichlet boundary condition of is prescribed at the inflow. The initial condition is everywhere on . As there is no analytical solution available, no convergence study is performed and only a representative approximation is shown in Fig. 30.
5.3.4 Advectiondiffusion on a hyperbolic paraboloid with bumps
Finally, an advectiondiffusion problem is considered on the hyperbolic paraboloid with bumps introduced in Section 5.2.3. The time interval is again , the advection velocity tangential to the manifold in the direction of , and the diffusion coefficient . The initial condition is given by . A representative result is seen in Fig. 31.
6 Conclusions
A higherorder accurate method for PDEs on surfaces is proposed. It enables a completely automatic workflow from the geometric description based on levelsets to the numerical analysis without any userintervention. This is an important advantage over other methods which are based on handcrafted surface meshes. Furthermore, these meshes are often composed by flat triangles leading to a loworder representation of the geometry and the resulting approximation of the BVP. Compared to methods which solve BVPs on all zeroisosurfaces at once by using volumetric elements, the proposed approach is more efficient as the effort scales with standard planar, twodimensional BVPs. Compared to methods which employ shape functions of the background mesh as in TraceFEM and CutFEM, it is an important advantage that boundary conditions are enforced in the standard way without additional technqiues for general constraints (Lagrange multipliers, Nitsche’s method etc.).
The proposed method is characterized by the following key ingredients:

A geometry description of the bounded manifold based on several levelset functions.

The automatic generation of surface elements which enable a continuous, higherorder accurate representation of the levelset geometry. A conforming surface mesh used for the finite element approximation is automatically extracted from this set of surface elements.

The manipulation of the background mesh by node movements in order to ensure the shape regularity of the resulting surface elements and a bounded condition number of the resulting system of equations. This step may later be replaced by a stabilization similar to what is done in the TraceFEM and CutFEM.
The numerical results confirm that higherorder accurate solutions of BVPs on curved surfaces in three dimensions are achieved. The next steps will be to investigate the proposed technique for the solution of more advanced BVPs on manifolds such as flow problems on curved surfaces and the numerical analysis of membranes and shells.
References
 [1] Adalsteinsson, D.; Sethian, J.A.: Transport and diffusion of material quantities on propagating interfaces via level set methods. J. Comput. Phys., 185, 271 – 288, 2003.
 [2] Bertalmio, M.; Cheng, L.T.; Osher, S.; Sapiro, G.: Variational problems and partial differential equations on implicit surfaces: The framework and examples in image processing and pattern formation. J. Comput. Phys., 174, 759 – 780, 2001.
 [3] Blaauwendraad, J.; Hoefakker, J.H.: Structural Shell Analysis, Vol. 200, Solid Mechanics and Its Applications. Springer, Berlin, 2014.
 [4] Burger, M.: Finite element approximation of elliptic partial differential equations on implicit surfaces. Computing and Visualization in Science, 12, 87 – 100, 2009.
 [5] Burman, E.; Hansbo, P.; Larson, M.G.: A stabilized cut finite element method for partial differential equations on surfaces: The LaplaceBeltrami operator. Comp. Methods Appl. Mech. Engrg., 285, 188 – 207, 2015.
 [6] Burman, E.; Hansbo, P.; Larson, M.G.; Larsson, K.; Massing, A.: Finite Element Approximation of the LaplaceBeltrami Operator on a Surface with Boundary. ArXiv eprints, Sep 2015.
 [7] Chapelle, D.; Bathe, K.J.: The Finite Element Analysis of Shells – Fundamentals. Computational Fluid and Solid Mechanics. Springer, Berlin, 2011.
 [8] Deckelnick, K.; Dziuk, G.; Elliott, C.M.; Heine, C.J.: An narrow band finiteelement method for elliptic equations on implicit surfaces. IMA Journal of Numerical Analysis, 30, 351 – 376, 2010.
 [9] Deckelnick, K.; Elliott, C.M.; Ranner, T.: Unfitted finite element methods using bulk meshesfor surface partial differential equations. SIAM J. Numer. Anal., 52, 2137 – 2162, 2014.
 [10] Delfour, M.C.; Zolésio, J.P.: A Boundary Differential Equation for Thin Shells. Journal of Differential Equations, 119, 426 – 449, 1995.
 [11] Delfour, M.C.; Zolésio, J.P.: Tangential Differential Equations for Dynamical ThinShallow Shells. Journal of Differential Equations, 128, 125 – 167, 1996.
 [12] Demlow, A.: Higherorder finite element methods and pointwise error estimates for elliptic problems on surfaces. SIAM J. Numer. Anal., 47, 805 – 827, 2009.
 [13] Demlow, A.; Dziuk, G.: An adaptive finite element method for the LaplaceBeltrami operator on implicitly defined surfaces. SIAM J. Numer. Anal., 45, 421 – 442, 2007.
 [14] Demlow, A.; Olshanskii, M.A.: An adaptive surface finite element method based on volume meshes. SIAM J. Numer. Anal., 50, 1624 – 1647, 2012.
 [15] Deuflhard, P.; Weiser, M.: Adaptive numerical solutions of PDEs. De Gruyter, Berlin, 2012.
 [16] Du, Q.; Ju, L.; Tian, L.: Finite element approximation of the CahnHilliard equation on surfaces. Comp. Methods Appl. Mech. Engrg., 200, 2458 – 2470, 2011.
 [17] Dziuk, G.: Finite Elements for the Beltrami operator on arbitrary surfaces, Chapter 6, 142 – 155. Springer Berlin Heidelberg, Berlin, Heidelberg, 1988.
 [18] Dziuk, G.; Elliott, C.M.: Finite elements on evolving surfaces. IMA Journal of Numerical Analysis, 27, 262 – 292, 2007.
 [19] Dziuk, G.; Elliott, C.M.: Eulerian finite element method for parabolic PDEs and on implicit surfaces. Interfaces and Free Boundaries, 10, 2008.
 [20] Dziuk, G.; Elliott, C.M.: Finite element methods for surface PDEs. Acta Numerica, 22, 289 – 396, 2013.
 [21] Elliott, C.M.; Stinner, B.: Analysis of a diffuse interface approach to an advection diffusion equation on a moving surface. Mathematical Models and Methods in Applied Sciences, 19, 787 – 802, 2009.
 [22] Elliott, C.M.; Styles, V.: An ALE ESFEM for Solving PDEs on Evolving Surfaces. Milan Journal of Mathematics, 80, 469 – 501, 2012.
 [23] Fries, T.P.; Omerović, S.: Higherorder accurate integration of implicit geometries. Internat. J. Numer. Methods Engrg., 106, 323 – 371, 2016.
 [24] Fries, T.P.; Omerović, S.; Schöllhammer, D.; Steidl, J.: Higherorder meshing of implicit geometries—part I: Integration and interpolation in cut elements. Comp. Methods Appl. Mech. Engrg., 313, 759 – 784, 2017.
 [25] Grande, J.; Reusken, A.: A Higher Order Finite Element Method for Partial Differential Equations on Surfaces. SIAM J. Numer. Anal., 54, 388 – 414, 2016.
 [26] Gross, S.; Reusken, A.: Numerical Methods for Twophase Incompressible Flows, Vol. 40, Springer Series in Computational Mathematics. Springer, Berlin, 2011.
 [27] Hansbo, A.; Hansbo, P.: A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comp. Methods Appl. Mech. Engrg., 193, 3523 – 3540, 2004.
 [28] Hansbo, P.; Larson, M.G.: Finite element modeling of a linear membrane shell problem using tangential differential calculus. Comp. Methods Appl. Mech. Engrg., 270, 1 – 14, 2014.
 [29] Hansbo, P.; Larson, M.G; Larsson, F.: Tangential differential calculus and the finite element modeling of a large deformation elastic membrane problem. Comput. Mech., 56, 87 – 95, 2015.
 [30] Kamilis, D.: Numerical methods for the PDEs on curves and surfaces. Masters thesis, Umeå university, Umeå, Sweden, 2013.
 [31] Olshanskii, M.A.; Reusken, A.: A finite element method for surface PDEs: matrix properties. Numer. Math., 114, 491 – 520, 2010.
 [32] Olshanskii, M.A.; Reusken, A.; Grande, J.: A finite element method for elliptic equations on surfaces. SIAM J. Numer. Anal., 47, 3339 – 3358, 2009.
 [33] Olshanskii, M.A.; Reusken, A.; Xu, X.: On surface meshes induced by level set functions. Comput. Visual Sci., 15, 53 – 60, 2012.
 [34] Olshanskii, M.A.; Reusken, A.; Xu, X.: A stabilized finite element method for advectiondiffusion equations on surfaces. IMA Journal of Numerical Analysis, 34(2), 732 – 758, 2014.
 [35] Olshanskii, M.A.; Safin, D.: Numerical integration over implicitly defined domains for higher order unfitted finite element methods. Lobachevskii Journal of Mathematics, 37, 582 – 596, 2016.
 [36] Rangarajan, R.; Lew, A.J.: Universal Meshes: A new paradigm for computing with nonconforming triangulations. arxiv:1201.4903, 2012.
 [37] Rätz, A.; Voigt, A.: PDEs on surfaces – A diffuse interface approach. Communications in Mathematical Sciences, 4, 575 – 590, 2006.
 [38] Reusken, A.: Analysis of trace finite element methods for surface partial differential equations. IMA J Numer Anal, 35, 1568 – 1590, 2014.
 [39] Xu, J.J.; Yang, Y.; Lowengrub, J.: A levelset continuum method for twophase flows with insoluble surfactant. J. Comput. Phys., 231, 5897 – 5909, 2012.
 [40] Xu, J.J.; Zhao, H.K.: An Eulerian Formulation for Solving Partial Differential Equations along a Moving Interface. Journal of Scientific Computing, 19, 573 – 594, 2003.