Higherorder meshing of implicit geometries—part I:
Integration
and interpolation in cut elements
Abstract
An accurate implicit description of geometries is enabled by the levelset method. Levelset data is given at the nodes of a higherorder background mesh and the interpolated zerolevel sets imply boundaries of the domain or interfaces within. The higherorder accurate integration of elements cut by the zerolevel sets is described. The proposed strategy relies on an automatic meshing of the cut elements. Firstly, the zerolevel sets are identified and meshed by higherorder interface elements. Secondly, the cut elements are decomposed into conforming subelements on the two sides of the zerolevel sets. Any quadrature rule may then be employed within the subelements. The approach is described in two and three dimensions without any requirements on the background meshes. Special attention is given to the consideration of corners and edges of the implicit geometries.
Keywords: Numerical integration, levelset method, fictitious domain method, XFEM, GFEM, interface capturing
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
The approximation of boundary value problems (BVPs) based on the finite element method (FEM) has achieved an enormous importance in engineering, physics, and related fields. The goal to achieve higherorder accurate approximations of BVPs can be traced back to the early days of the FEM, see e.g. [58] or the text books [52, 5, 7, 64]. This is typically labeled FEM and there are numerous references found. In fact, the isoparametric concept leads to a conceptually simple approach to achieve optimal results. Therefore, an accurate geometry representation based on higherorder elements which consider boundaries and interfaces is necessary, see Fig. 1(a) for an example. The generation of such meshes is, however, not a simple task. Moreover, elements may have to be refined during the analysis, for instance in the context of adaptivity and convergence studies. The interplay of the FEM software, the meshing tool, and the CAD program is far from trivial and hardly automated, in particular using higherorder elements. Things become even worse in case of moving boundaries and interfaces.
Hence, it is not surprising that a number of methods have been developed which do not require conforming meshes. In the case where boundaries are not meshed explicitly, we use the term “fictitious domain methods” (FDMs) and there is a large variety of variants. These methods are characterized by the fact that the physical domain in which a BVP is formulated, is completely embedded in a fictitious domain, also called background mesh, see Fig. 1(b). The background mesh is often structured and no geometry dependent meshing is needed. The shape functions are provided by the background mesh. Methods falling into this class include the unfitted or cut finite element method [9, 10, 8, 23], finite cell method [2, 15, 43, 47, 48], Cartesian grid method [59, 62], immersed interface method [30], virtual boundary method [45], embedded domain method [32, 40] etc.
Another class of methods has been developed for problems with interfaces rather than boundaries. Often, simple meshes are chosen to describe the overall geometry without taking care that internal interfaces such as material interfaces or crack surfaces are meshed conformingly, see Fig. 1(c) and (d). Nonstandard approximation spaces are often constructed near the interfaces. In this class, we mention extended or generalized finite element methods (XFEM/GFEM), see e.g. [6, 34, 17] for the XFEM and [55, 56] for the GFEM. They consider for innerelement jumps and kinks by adding enrichment functions based on the partition of unity concept [3, 4, 33]. The XFEM has also been used in the realm of FDMs e.g. in [57]. Other mesh independent approximations for nonsmooth solutions involve the weak element method [20, 44], the ultra weak variational method [11], the leastsquares method as proposed in [35], and the globallocal FEM as described in [36, 41]. Finally, we mention the manifold method [12, 51] which may be used for innerelement discontinuities.
A shared property of FDMs and XFEMrelated approaches is the numerical integration in elements cut by boundaries or interfaces. Approaches for the integration of elements with internal boundaries and interfaces may be distinguished based on the fact whether they rely on a decomposition of the cut elements into subelements or not. For the first class, the standard approach is to recursively decompose a cut element into polygonal subcells until the desired accuracy is obtained [1, 37, 14]. This typically leads to a very large number of integration points in the context of higherorder approximations [54, 63, 28, 26, 14]. It was already noted in [27, 13, 18, 46] that a decomposition into subelements with curved, higherorder edges or faces is a strategy which consistently enables the generation of higherorder accurate integration rules with only a modest number of integration points. Of course, the effort for generating these integration points is larger than for the polygonal approaches with reduced accuracy. The other class of approaches is built by methods that do not decompose the cut elements [60, 61, 39, 38]. Typically, the number of integration points is small, however, the generation of proper integration weights is often involved and requires the solution of small systems of equations e.g. in the context of moment fitting and Lasserre’s technique. The extendability to three dimensions, higherorder accuracy, general integrands, and to the presence of corners and edges often poses problems in this class of approaches.
The higherorder accurate numerical integration in threedimensional, cut elements is subject of some very recent contributions and interested readers are referred to the references given therein: For those approaches relying on element decompositions, [18] presents a very general approach in the context of the levelset method which applies for all element types and has no restriction on the background meshes. The adaption to the finite cell method and spline based boundary representations is reported in [53]. “Smart octrees” are proposed in [25] and are employed for structured quadrilateral and hexahedral background meshes. Among the approaches which avoid decompositions, [24] presents an approach based on moment fitting. Finally, the approach in [29] manipulates the background mesh such that the interfaces are meshed conformingly. Comparisons of these advanced approaches with the classical use of (polygonal) recursive spacetrees are reported in [25, 53] and clearly demonstrate the advantages of the recent approaches.
Herein, we propose a new strategy for the numerical integration of elements cut by boundaries and interfaces which relies on the decomposition of cut elements into isoparametric, higherorder subelements. Let us assume the abstract domain represented by Fig. 2(a) featuring an internal interface implicitly defined by the levelset method. It does not matter here whether is an interface between two materials in and or a boundary if e.g. is considered to be a void region. A higherorder background mesh is introduced as shown in Fig. 2(b) and the levelset data is only given at the nodes. The implied zerolevel set is, in general, curved within the elements, see Fig. 2(c). The task is to integrate on one or both sides of the interface individually with higherorder accuracy. Therefore, the zerolevel sets are approximated by higherorder interface elements (reconstruction) and higherorder subelements result on the two sides of the interface (decomposition). See Fig. 3(a) for the resulting subelements in the overall mesh and Fig. 3(b) for a detail. Standard integration points, e.g. Gauss points may now be mapped into the subelements and used in FDMs or XFEMrelated methods for the integration in the corresponding background element, see Fig. 3(c). It is also shown how corners and edges may be considered properly by the proposed method using several levelset functions and performing the decomposition successively for each. Note that previous works in [18, 16] by (partly) the same authors do not offer this simple concept for corners and edges. Furthermore, we find that the approach herein unifies the treatment of the reconstruction and decomposition step as they now both rely on element generation, i.e. meshing.
This is the first paper in a sequence of contributions where it is shown that the resulting, automatically generated meshes, herein only used for the proper placement of integration points in cut elements and interpolation, may also be used for directly approximating BVPs. In the second part, this will be shown for BVPs on zerolevel sets, i.e. on manifolds. Therefore, only the mesh of reconstructed interface elements will be needed. In the third part, BVPs with innerelement boundaries and interfaces are considered using meshes consisting of regular and cut background elements, the latter decomposed into conforming subelements as described herein. Further studies will then outline how fully automatic, higherorder analyses are enabled for domains in the context of constructive solid geometries.
The paper is organized as follows: Section 2 describes some preliminaries and defines the higherorder background meshes used herein and the levelset method. In Section 3, the identification and meshing of the zerolevel set is described in detail (reconstruction). The decomposition into higherorder, conforming subelements is detailed in Section 4. The proper consideration of corners and edges with several levelset functions is outlined in Section 5. Numerical results are presented at the end of each section. Finally, the paper concludes in Section 6 with a summary and outlook. An appendix is provided where mappings are detailed which are frequently used throughout this work.
2 Background mesh and levelset data
A two or threedimensional domain of interest, , is defined implicitly based on the levelset method [42, 50]. That is, a scalar function , is given, called levelset function, which is positive on one side of the boundary/interface and negative on the other. The function is continuous and its zerolevel set defines the position of a smooth boundary or interface. Note that does not have to be a signed distance function herein. In Section 5, we shall also consider the situation where boundaries and interfaces are, in fact, not smooth but feature edges and corners which is of high practical relevance. Then, we suggest to use several levelset functions for the description.
The domain is fully immersed in a background mesh composed by higherorder Lagrange elements. The background mesh may be unstructured and feature curved elements. The focus herein is on triangular elements in 2D and tetrahedral elements in 3D. The levelset data is only given at the nodes of the background mesh, , and is interpolated by standard finite element shape functions inbetween, hence, . That is, once the nodal levelset data is given, no (further) use of any analytic knowledge of the levelset functions is made. Obviously, the background mesh enables a higherorder accurate description of the domain as long as the boundaries and interfaces are sufficiently smooth.
The task is to obtain integration points which allow a higherorder accurate integration in the implicit geometry. Therefore, the background elements are decomposed into subelements which conform with the interfaces and automatically consider also for corners and edges. It is then simple to use any desired integration rule (e.g. Gauss quadrature) within each subelement. The decomposition is successively applied in each element for the involved levelset functions. This is a twostep procedure for each zerolevel set that cuts an element: (1) Reconstruction, which is the meshing of the zerolevel set with higherorder interface elements. (2) Decomposition of cut elements into higherorder subelements on the two sides of the interface elements. These two steps are depicted in the left half of Fig. 4(a) and (b) for the two and threedimensional case, respectively. It is important to note that both steps are performed first in each cut reference element. This is, in fact, crucial for the ability to handle arbitrarily curved background elements and perform the decomposition successively for several levelset functions. Thereafter, (3) the subelements are mapped to the corresponding physical background element and (4) any desired quadrature points are mapped to these subelements, see the right half of Fig. 4.
Herein, the focus is on triangular and tetrahedral elements, but the extension to quadrilateral and hexahedral elements is truly straightforward and described in [18] in a similar context. For example, one may simply decompose a reference hexahedral element into tetrahedra which is not a problem because the faces are flat in the reference domain. Then, the decomposition is realized in each reference tetrahedron as described below.
3 Reconstruction
Let us focus on a reference element which is cut by the zerolevel set of a levelset function. The nodal values of the levelset function at the element nodes are taken from the physical element, i.e. . The task is to approximate the zerolevel set of by means of a higherorder interface element. That is, roots have to be identified on the zerolevel set which are then the element nodes of the interface element. Inbetween, the zerolevel set is only approximated. It is important to note that the location of the inner element nodes is not unique as shall be seen below.
The description here is partly along the lines of [18, 16], however, the following issues in the context of reconstructions are new: (i) The generation of start values for the rootfinding based on (cubic) Hermite interpolation in triangular elements, (ii) the use of tailored mappings for predicting suitable start values in tetrahedral elements, leading to optimal accuracy in contrast to results reported in [18, 16], and (iii) the study of integration and interpolation properties of the resulting reconstructed surface elements.
3.1 Valid levelset data and recursion
It is obvious that higherorder elements may feature very complicated topologies of zerolevel sets. For example, element edges may be cut several times, or the zerolevel set is completely inside an element, see Fig. 6 in 2D and Fig. 7 in 3D. Therefore, a recursive procedure may be required (typically only in very few elements) until “valid levelset data” inside the element is obtained. In two dimensions, by “valid” we refer to the situation where (i) each element edge is only cut once, (ii) the overall number of cut edges must be two, and (iii) if no edge is cut then the element is completely uncut. In three dimensions, the conditions (i) to (iii) are checked on each element face plus the additional condition that (iv) if no face is cut then the element is completely uncut. Whether these conditions are fulfilled is checked based on a sample grid in the reference element, see Fig. 5. Only the sign of the levelset data at the sample points is considered for these checks.
If the levelset data is not valid the reference element is recursively refined and the described procedures apply to the refined elements, see Figs. 6 and 7. The distinction between valid and invalid levelset data is done before the reconstruction and decomposition steps. However, the validity check does not necessarily ensure that the subsequent reconstruction and decomposition are successful. For example, when a resulting subelement after the decomposition features a negative Jacobian, a recursive refinement is still triggered although the initial levelset data was possibly valid in the above sense. This may, for example, be the case for strongly curved zerolevel sets, see Fig. 6(d). A levelset function must not be zero right at the corner nodes of the background element; if that is the case, the levelset value is simply perturbed in the range of without any noticeable change in the results.
3.2 Topological cases
From now on, it is assumed that valid levelset data is present in a (refined) reference element. Then, the zerolevel set cuts the reference element in a limited number of topologically different cases, see Fig. 8. Only the signs at the corner nodes of the reference element determine the situation. A cut triangle is decomposed into one subtriangle and one subquadrilateral; the latter may be further decomposed into subtriangles. A cut tetrahedron is decomposed either into one subtetrahedron and one subprism (topology 1), see Fig. 8(b), or into two subprisms (topology 2), see Fig. 8(c). Obviously, the subprisms may be further decomposed into tetrahedra. There is no need to avoid hanging nodes in the decomposition into subelements for integration purposes.
3.3 Root search and interface elements
Depending on the topological situation, different algorithms for the identification of the element nodes of the interface elements are employed. Each node is a root of the levelset function which is obtained by solving a nonlinear problem. The root finding is characterized by (i) the start values of the iteration, (ii) the search directions in the reference element, and (iii) the iteration method. The discussion is along the lines of [18].
3.3.1 Start values and search directions in 2D
The first step is to identify the intersection of the zerolevel set with the element edges. In triangular reference elements, these are two points. One may then define a linear interpolation of these two points or a Hermite interpolation which takes into account also the direction of the zerolevel set at the intersection points. In order to obtain the start guesses, an interface element of the desired order is mapped onto the linear or Hermite intermediate reconstruction. These are the blue lines and crosses in Fig. 9 (a) to (d) for the linear and (e) to (h) for the Hermite case.
For the search directions, four different variants are studied:

Towards the node on the other side than the other two, see Fig. 9(a) and (e).

In direction of the interpolated edge directions at the intersection points, see Fig. 9(b) and (f).

In normal direction to the linear or Hermite reconstruction, see Fig. 9(c) and (g).

In direction of the gradient of the levelset function, i.e. , see Fig. 9(d) and (h).
For later reference in the numerical studies, the 2D variants are characterized by twodigit numbers: The first digit is either for a linear or for a Hermite reconstruction. The second digit refers to the four variants of the search directions listed above.
3.3.2 Start values and search directions in 3D
In three dimensions, the element nodes of the higherorder interface element (triangle for topology 1, quadrilateral for topology 2) are separated into outer and inner nodes. The outer (i.e., edge) nodes are enforced to remain on the faces of the reference tetrahedron. In fact, this is achieved by treating one face after the other and performing a reconstruction in reference triangles as described above. Thereafter, these nodes are mapped to the corresponding face and define the edge nodes of the interface element.
For the inner nodes, the start values and search directions depend on an intermediate reconstruction of the zeroisosurface based on the mapping defined in [52]. This mapping defines a surface implied by the three or four curved edges (higherorder line elements) of the triangular or quadrilateral zeroisosurface within the cut tetrahedron, see Fig. 10. The definition of the mapping is outlined in the appendix 7.1 of this work. Next, a higherorder reference interface element is mapped onto the intermediate reconstruction to obtain the start values, see the blue crosses in Fig. 10. The search directions are then, as above, either the corresponding normal vectors or the gradients of the levelset function.
It is thus seen that also in 3D, the 2Dvariants play a major role as they fully define the outer nodes and are also important to define start values for the inner nodes through the intermediate reconstructions. It is noted that, in contrast to the definition of the start values in 3D as proposed in [18], optimal convergence rates are achieved with the procedure proposed herein.
3.3.3 Iterative procedure
An iterative procedure is required to identify positions on the zerolevel set from the starting points. One approach is to move strictly along the straight search paths as defined above, see Fig. 11(a) and (b). The other approach rather uses the gradient of at each intermediate position during the iteration, see Fig. 11(c). All approaches yield quadratic convergence rates of the iterative procedure but, of course, to (slightly) different positions of the element nodes on the zerolevel set.
Mathematically, the algorithm is described as follows. The task is to find the root of in the reference element, i.e. some position on the zerolevel set. The starting point of the iterative procedure is labeled . The NewtonRaphsontype algorithm for all approaches considered here is based on the following iteration:
(3.1) 
is the search direction. In the case of a straight search path, this is e.g. the normal vector of the intermediate reconstruction or the gradient at each starting point, hence . Otherwise, it may also be the current gradient .
3.4 Numerical results for reconstructions in 2D and 3D
In order to distinguish the quality of the different variants to recover reconstructions, some numerical results are presented resulting from an integration and interpolation on the zerolevel sets. This is also justified by the fact that many properties of the second step, i.e. the decomposition, as outlined in the next section, are immediately inherited from the reconstruction. If the reconstruction were not optimal, the resulting integration points in cut elements would also be suboptimal.
It is emphasized that in all studies presented herein, the respective orders of the background elements and the reconstructed surface elements coincide. Later, in Section 4, the same holds for the order of the resulting subelements after the decomposition. As confirmed in [18], using different orders would not be useful as the overall convergence rates would be determined by the lowest order.
3.4.1 Numerical results for zeroisolines in 2D
We start with the twodimensional situation. Following [16], two different levelset functions are considered:
(3.2)  
(3.3) 
with , and . See Fig. 12(a) for a visualization of the zerolevel sets of and . The zerolevel set of is a circle with radius and is frequently used in the literature, e.g. [18], and is similar to [31]. The function is defined as
(3.4) 
and is later integrated on the zerolevel sets.
In the convergence studies, elements are employed per dimension. We use standard Gauss rules for the integration on the zerolevel set with a rather high order of which is kept constant independently of the order of background and interface elements. Results in 2D are studied in different error norms which are related to integrating and interpolating a function on [16]:
(3.5)  
(3.6)  
(3.7)  
(3.8)  
(3.9) 
where are the integration points in the interface elements with integration weights . It is noted that, in contrast to the nodes of the interface elements, the integration points at are only approximately on the zerolevel set of .
The first error, in Eq. (3.5), is evaluated by summing up integration weights giving the length of the reconstructed interface which is compared with the exact length (for , this is the circumference of a circle). The second error, , integrates the levelset function on the zerolevel set, which, in the ideal case, would be zero. Nevertheless, because the nodes of the interface elements are on the zerolevel set of rather than and, furthermore, the integration points are only approximately on the zerolevel set of , is not exactly zero. The other error norms integrate either the exact or the interpolated function as given in Eq. (3.4). For , the function is evaluated exactly at the integration points . Instead, evaluates the function at the element nodes of the interface elements and interpolates at the integration points using the shape functions of the interface element. Equivalently, uses the element nodes of the background elements to evaluate the function and the corresponding shape functions to interpolate.




Results are shown in Fig. 13. In Fig. 13(a), it is seen that most of the variants for the search algorithm (start values and search directions) perform similar and with optimal convergence. Variants and , where the search directions are based on the gradient of the levelset function are split into two subversions 14a/14b and 24a/24b depending on whether the search path is constant, , or changes during the iteration, . A few variants actually underperform, in particular variant and . In these variants, the search directions depend on the topological situation only, i.e. the search directions are towards the one node on the other side or the interpolated edge directions but are, otherwise, independent of the levelset data. It is seen that there is only a negligible difference between the variants which are based on the linear reconstruction (first digit in the variant number is ) compared to the Hermite reconstruction (first digit is ). Because it is found that this also applies to the threedimensional situation, the Hermite reconstruction is neglected in the following. Of course, Fig. 13(a) is just one sample of a large parameter space that was actually studied where, as an example, is chosen and an order of for the background mesh and the reconstruction.
Fig. 13(b) shows that the error converges independently of the error norms in Eqs. (3.5) to (3.9). This was again confirmed for the whole parameter space that has been studied for the two levelset functions. Now that it is known which variants perform potentially optimal independently on the error norm, the full convergence results for different orders are shown in Fig. 13(c) and (d) for and , respectively. The figures look very similar for all optimal variants. It is confirmed that optimal convergence rates are achieved. It is noted that only data points are displayed where no recursive refinement was necessary (because this influences the element lengths and complicates the proper placement of the result in the convergence plots).
3.4.2 Numerical results for zeroisosurfaces in 3D
Numerical results are now presented for the integration on zeroisosurfaces in 3D. The procedure is similar to the studies in two dimensions. It is recalled that in 3D, inner and outer nodes of each interface element are treated differently. For the outer nodes, any of the 2D variants may be used, see Section 3.3.1. Here we restrict ourselves to variant but all other variants have been systematically studied as well. For the inner nodes, three variants are studied, see 3.3.2: Variant A and B use straight search paths, the first based on the normal vector of the intermediate reconstruction at the start value and the second based on the gradient of the levelset function at . Variant C is based on the gradient as well, but depending on the points during the iterative procedure.


The following two levelset functions are considered in 3D:
(3.10)  
(3.11) 
with . The zeroisosurfaces of these two levelset functions are visualized in Fig. 14. There, also a coarse tetrahedral background mesh is shown. For the integrand, the function
(3.12) 
is chosen. In the convergence studies, elements are used per dimension.
Again, the error is studied in different norms. The three norms in Eqs. (3.5) to (3.7) are directly applicable in this 3D study as well, the other two are slightly modified as
Therein, is the integrand interpolated by the shape functions of the reconstructed interface elements and is interpolated by the shape functions of the background elements.
Results are shown in Fig. 15. Fig. 15(a) shows that all three variants for the inner nodes lead to virtually identical results. Therefore, we concentrate on the variant A using the normal vector for the root search. This is consistent to variant 13 for the 2D case, here applied on each of the faces of the tetrahedron. Fig. 15(b) confirms that the same convergence rates are expected in all of the five error norms. Finally, Fig. 15(c) and (d) vary the order and resolution of the background mesh to show that optimal convergence rates are obtained for integrations on the zeroisosurfaces of and .




Hence, optimal convergence rates are found for the reconstruction regardless of the fact that a number of errors are involved: The reconstruction is only an approximation of the zerolevel set of the interpolated levelset which, in fact, is only an approximation of . In some error norms, only the interpolated function is integrated rather than . Finally, also the numerical integration based on Gauss quadrature introduces an integration error. Nevertheless, it was shown that is integrated and interpolated optimally on zerolevel sets implied by .
4 Decomposition
Once the interface elements are defined, i.e. the zerolevel set is meshed with higherorder interface elements, the cut background elements are to be decomposed into subelements depending on the topology of the cut situation, see Section 3.2. This is outlined in Fig. 4 for the two and threedimensional situation, respectively. The resulting subelements share the property that they feature one higherorder side coinciding with the reconstructed interface element. All other edges are straight. A key aspect for the decomposition is the mapping of element nodes into the special subelements on the two sides of the interface.
4.1 Decomposition in 2D



For the twodimensional case, a cut reference background element in the coordinate system falls into a triangular and a quadrilateral subelement with one higherorder side, see Fig. 4(a). In order to define the element nodes of the subelements, a mapping from a higherorder reference triangle or quadrilateral in coordinates to each subelement in coordinates is sought, see Fig. 16(a) and (b) for an example. Such a mapping is, as before in Section 3.3.2 in the context of reconstructions, not unique. However, the properties of this mapping largely influence the convergence properties of the resulting subelements.
This is outlined for the case of a triangular subelement where we introduce three different mappings . That is, the reference element in is mapped differently to the coordinate system as shown in Fig. 16(c). Note that the outer contour is the same for the three mappings, however, different for the inner region. Mapping 1 is defined in [52] and is outlined in the appendix 7.1 for the case of three curved edges in . The situation here is the reduced case where the mapping is to and only one of the edges is curved (the diagonal edge 2 of the reference triangle). Nevertheless, the same formulas given in the appendix may be employed, it follows from Eq. (7.1) that only for .
A second mapping is the original blending function mapping proposed in [21, 22, 58] naïvely extended to triangles. Starting point is the approach in the appendix restricted to one curved edge. We set and for . For , we define , , and follows from Eq. (7.1). The third mapping is the intersection mapping of [49, 19] and is motivated from geometrical considerations. In the reference element, auxiliary points on the element boundary are computed for any given inner point . These points are mapped to the subelement in which is not a problem because they are on the outer contour. They imply two straight lines and the intersection of the two lines is the mapped inner point .



Although Fig. 16(c) indicates that all three mappings successfully map to the linear triangle with one higherorder side, it is found that the second map yields suboptimal accuracy in the convergence studies. This could be traced back to the fact that the mapping is not as smooth as needed for higherorder accuracy. See Fig. 17 for visualizations of the Jacobians of the three mappings. It is seen that the second variant does not feature a smooth Jacobian near the node opposite to the higherorder side which explains its suboptimal performance. The intersection mapping and the map based on [52] and outlined in the appendix feature smooth Jacobians and perform optimal. However, the intersection mapping is not easily extended to the threedimensional situation (three straight lines do not necessarily intersect in one point). Therefore, in the following we restrict ourselves only to the mappings in two and three dimensions as described in [52]; they are general, smooth and lead to optimal element nodes in the subelements as shall be seen below.
4.2 Decomposition in 3D
In three dimensions, it was outlined in Section 3.2 that a tetrahedron may fall into two different topologies. For topology 1, the tetrahedron was cut into a subtetrahedron and a subprism where the curved face is triangular, see also Fig. 8(b). For topology 2, the cut tetrahedron is decomposed into two prisms where one of the quadrilateral faces is curved in each subprism, see Fig. 8(c). The curved faces coincide with the reconstructed higherorder interface elements as discussed in Section 3. Each subelement has straight edges except for those belonging to the curved face.
Again, the task is to define maps (of element nodes) from higherorder reference tetrahedra and prisms in coordinates to the cut reference background element in . For some possible situations, examples are shown in Fig. 18 for all types of subelements. For the subelements with one higherorder face, the corresponding mappings are outlined in appendix 7.2 for tetrahedra and in appendix 7.3 for prisms. The formulas are based on [52] and adapted to the present situation. The resulting formulas are quite involved wherefore they are moved to the appendix.



4.3 Map of integration points
Once the cut background elements are successfully decomposed into subelements in the reference domain, one may simply map the obtained element nodes to the physical background element using the isoparametric concept. It is then simple to map integration points according to any desired quadrature rule to the subelements. Depending on the particular context for which higherorder accurate integration points are sought, this may be preferred in the subelements of the cut reference or physical background element. For the integration studies performed herein, the integration points are generated in the physical background elements. However, in the context of FDM and XFEMrelated methods one would generate integration points first in the reference background elements, then evaluate the corresponding shape functions and finally map them to the physical background elements using the isoparametric concept.
4.4 Numerical results for decompositions in 2D and 3D
The numerical results presented here are a direct extension of the studies of Section 3.4 to the integration and interpolation of implicitly defined areas and volumes.
4.4.1 Preliminary examples
Let us first consider some examples of decompositions. Fig. 19 shows some resulting subelements for complex levelset data in 2D. According to Section 3.1, recursive refinements are necessary to obtain valid levelset data in the refined elements. Note the relation to Fig. 6. Examples for decompositions in 3D are shown in Fig. 20 which correspond to zerolevel sets visualized in Fig. 7.








4.4.2 Numerical results for decompositions in 2D
The setup is identical to the one described in Section 3.4.1. The same background meshes, levelset functions and integrands are considered here. Three different (relative) error norms are introduced for the integration in :
(4.1)  
(4.2)  
(4.3) 
where are 2D integration points in the special subelements and the corresponding weights.




Convergence results are shown in Fig. 21. Fig. 21(a) confirms that most variants of the search directions perform optimally. A very similar result was shown in Fig. 13(a) for the integration on the zerolevel sets and the same conclusions are drawn: Firstly, search directions should be based on the normal vectors on the intermediate reconstruction or on the gradient of the levelset function. Secondly, Hermite reconstructions are not superior to linear ones. Therefore, we suggest variant 13 (based on normal vectors on linear reconstructions) for its simplicity. Fig. 13(b) indicates that all three error norms lead to identical results which was confirmed for the whole parameter space tested. An exception was only found for the suboptimal second variant of the mapping to the subelements described in Section 4.1: There, it turned out that results are suboptimal in (convergence rates are bounded by 4). Optimal results for and were still possible depending on the applied quadrature rules (which is also not desirable). Therefore, as mentioned before, we shall only present results where the mappings based on [52] are used, see also the appendix. Finally, the full convergence data for variant 13 for and according to Eqs. (3.2) and (3.3) is shown in Fig. 21(c) and (d): Optimal convergence rates are obtained. Only data points are shown where no recursive refinements were needed so that there is a clear element length associated to every result.
4.4.3 Numerical results for decompositions in 3D
In 3D, the convergence studies follow Section 3.4.2. The same definitions of the error norms from Eqs. (4.1) to (4.3) are used. Results are shown in Fig. 22 where the same form of presentation is employed than in Section 3.4.2. Fig. 22(a) uses variant 13 on the faces of the tetrahedra and variants A, B, C for the inner nodes. The differences are negligible and we suggest to use variant A because it is based on the normal vector of the intermediate reconstruction just as variant 13. Of course, all other combinations of search directions for the inner and outer nodes of the interface elements have been investigated as well: Those variants for the outer nodes which performed optimal in Section 4.4.2 also perform optimal in this 3D context. The fact that all three error norms behave similar is confirmed in Fig. 22(b). Finally, the full convergence data (for variant A13) for and according to Eqs. (3.10) and (3.11) is shown in Fig. 22(c) and (d) and, again, optimal convergence rates are obtained. It is noted that a previous work of the authors [18] was not able to achieve optimal convergence rates in 3D because the intermediate reconstructions used there were not sufficiently smooth.




It is confirmed for all studies that properties of the reconstruction are inherited by the decomposition. In particular, suboptimal results for reconstructions can not be better for decompositions. Of course, optimal reconstructions may still lead to suboptimal decompositions e.g. if the maps to the subelements, as described above, are not sufficiently smooth.
5 Corners and edges
So far, the reconstruction and decomposition of background elements has been described with respect to one levelset function. It is possible to define a levelset function whose zerolevel set also involves corners and edges. Then, this function is necessarily continuous. However, the meshing procedures described in the previous sections rely on an interpolated levelset function . Provided that is only continuous, an accurate interpolation poses requirements on the background mesh that we wish to avoid here. For example, in 2D, for being able to capture a corner, one would have to impose the constraint that the corner is at least on the edge of the background mesh. Otherwise, there would be a dramatic loss in accuracy in the interpolation and a higherorder decomposition would not be possible.
Therefore, we suggest to use multiple levelset functions for the domain description. Then, several levelset values are present at the nodes of the background mesh, one for each smooth part of the interface. In two dimensions, two levelset functions define a corner. See Fig. 23 for an example. As seen in Fig. 23(a), the domain features an interface/boundary with kinks, i.e. corners. Fig. 23(b) shows the two zerolevel sets and of the levelset functions and , respectively. Assume the levelset functions are negative inside the circular zerolevel sets. It is seen in Fig. 23(c) how the two levelset functions imply 4 different subdomains based on their signs. The gray region in Fig. 23(a) is . It is thus seen, how two levelset functions are able to define innerelement corners in 2D. No additional requirements on the background mesh are needed to capture the corners accurately.




In three dimensions, two levelset functions are able to define an edge and three functions a corner. See Fig. 24 for a graphical representation. Three zerolevel sets of , , and are shown in Fig. 24(a). It is seen how the intersections of any two levelset functions implicitly define edges and how three levelset functions imply corners. Of course, based on the signs of the three levelset functions, one may distinguish subregions of the background mesh. Depending on the application, one may now want to compose the domain based on these subregions. Two different examples are shown in Figs. 24(b) and (c); both feature implicit edges and corners.



It is now clear how multiple levelset functions are able to define edges and corners. It remains to describe how the corresponding reconstructions and decompositions are realized. The task is to individually mesh the subregions based on the signs of the involved levelset functions. The resulting meshes are already visualized in Figs. 23(c) and (d) in 2D and Figs. 24(b) and (c) in 3D. Thereafter, it is simple to assemble the domains of interest (e.g. for integration, interpolation, or approximation) from these subregions.
Fortunately, the decomposition of background elements with respect to several levelset functions is straightforward. The key aspect is that the decomposition is always carried out in reference background elements. That is, a background element is first decomposed with respect to the first levelset function . Then, the remaining levelset functions are interpolated at the new element nodes of the subelements of the cut background element. The decomposition with respect to the next levelset function is carried out just as if the already decomposed mesh were the initial mesh. Therefore, it may be useful to further decompose quadrilateral subelements in 2D into triangles and prismatic subelements in 3D into tetrahedra, so that the newly generated mesh (conforming to the zerolevel sets of previous ) only consists of the same element type. Fig. 25 shows an example where a reference element is decomposed with respect to level set functions. The resulting decompositions for each new levelset functions are shown and the different subregions that may be identified by the signs of the levelset functions are colorcoded. Fig. 25(f) shows a subdomain which features 3 corners within this one element. An example for the situation in 3D is shown in Fig. 26.












Hence, the decomposition with respect to several levelset functions is a truly straightforward successive application of the meshing procedures described above. The ability to naturally consider an arbitrary number of levelset functions, and even the presence of multiple edges and corners in one element seems to be a unique feature of the proposed method.
6 Conclusions
Implicit geometries occur in many applications, in particular in the context of fictitious domain methods and XFEMrelated methods. The boundaries of a domain of interest or interfaces inside the domain are then easily defined based on the levelset method. Using a higheroder background mesh enables a straightforward path to a higherorder accurate implicit description. However, during an analysis, integration points are required, typically, for the integration of the weak form of a model and accurate quadrature rules in cut background elements become a critical issue. Herein, a consistent path to obtain higherorder accurate integration points in cut background elements is described. There are virtually no requirements on the background mesh, in particular, there is no need to use Cartesian meshes or elements with straight edges.
The key idea is to generate higherorder finite elements on the two sides of an interface which is called decomposition. Before, the zerolevel sets have to be identified and approximated by interface elements which is called reconstruction. In fact, the reconstruction and decomposition may be seen as two meshing steps. A number of different variants to realize the meshing have been studied. It is found that one has to carefully define the start values and search directions used for the nonlinear detection of the zerolevel sets. For the decomposition, mappings of element nodes into the special subelements with one higherorder side are needed and they have to be sufficiently smooth. It is noted that depending on the levelset data inside the elements one may have to recursively refine background elements until valid levelset data is obtained. Recursive refinements are also needed when the (initial) reconstruction or decomposition fails e.g. due to negative Jacobians of the resulting subelements. When the curvatures of the levelset functions are reasonably adjusted with the resolution of the background elements, recursive refinements are only needed in a very small number of elements. The number of generated subelements per background element is then small, thus also the number of additional integration points in the cut elements.
The proposed meshing procedures will be used in upcoming publications to approximate boundary values problems without using typical FDM or XFEMrelated methods. In fact, because the generated elements align with boundares and interfaces, the mesh is used in a classical FEM context. Major issues are then the quality of the automatically generated elements which are not necessarily wellshaped. Nevertheless, it shall be seen that higherorder convergence rates are possible with such meshes not only in pure integration and interpolation problems as shown here but also in the approximation of boundary value problems.
The implementation of the proposed meshing procedures is not without efforts, especially in three dimensions. Therefore, we plan to soon enable a download of the underlying software on the institute’s webpage at www.ifb.tugraz.at.
7 Appendix
Crucial ingredients of the proposed remeshing are the maps of element nodes (or start values) to the element interiors based on the situation on the outer contour. For example, for the subelements on the two sides of the zerolevel sets, the outer contour is given by linear elements with one higherorder, typically curved side. For the placements of start values for the reconstruction in tetrahedra, three or four higherorder edges of the sought interface element are prescribed. In the following, the maps needed herein are defined in detail for triangles, quadrilaterals, tetrahedra and prisms. Because the definitions of these maps are quite lengthy, we decided to move this to the appendix. A more general assessment of general transfinite mappings is given e.g. in [52].
7.1 Mappings for triangles and quadrilaterals
Assume that a triangular (or quadrilateral) element is implied by (or ) curved line elements defining the outer contour, i.e. the edges. The line elements are of the same order, i.e. the same number of nodes and corresponding shape functions is associated with each edge. The task is to define a smooth mapping for all points in a reference triangle (or quadrilateral) to the threedimensional situation in . Obviously, this defines the shape of the triangular (or quadrilateral) surface in 3D, see Fig. 27. Such a mapping is not unique and the definition outlined below follows the textbook [52] adapted to the present situation.




The nodal coordinates of the line elements are denoted as with and . They must build a closed contour, i.e. the start and end coordinates must match properly. The (bi)linear shape functions related to the reference triangle (or quadrilateral) are called , more precisely, for triangles
and for quadrilaterals
The higherorder shape functions of the line elements are denoted as and are defined in a onedimensional reference element with coordinates . Herein, these are standard Lagrange shape functions. Furthermore, the linear shape functions and are needed. Local coordinates along each of the edges of the triangle or quadrilateral are introduced, see Fig. 27(a) and (c). They have to be related to the coordinates . For triangles, this gives
and for quadrilaterals
Next, we define
(7.1) 
where the first term on the right hand side defines the curve in , and the other two subtract the linear interpolant. Furthermore, a ramp function is needed which is, for triangles,
(7.2) 
and, for quadrilaterals,
(7.3) 
The overall map , i.e. the shape of the resulting surface element which is implied by the curved line elements, is then defined as
The first term on the right hand side is the linear interpolant of the corner points of the triangle (or quadrilateral) and is extracted from the start and end points of the curved line elements. See Fig. 27 for a graphical representation.
7.2 Mapping for tetrahedra
Let there be a tetrahedron with one curved face given by a higherorder triangular element with nodes. The more general case where all four faces are curved is not needed herein. The curved face has three curved edges, which are higherorder line elements with nodes each. The other three edges are straight. The situation is depicted in Fig. 28. The task is to define a smooth mapping for all points in a reference tetrahedron to the threedimensional situation in .



The nodal coordinates of the triangular face element are denoted as with . The three corner nodes of the face are set as the nodes with of the tetrahedron. The position of the remaining node of the tetrahedron must be provided as well. It shall be seen that the resulting mapping may be decomposed into three parts, which may be associated to vertex, face and edge contributions, hence,
(7.4) 
The first part is simply the linear mapping of a tetrahdron, i.e.
(7.5) 
with the linear shape functions
For the edge contributions, we proceed similar as in Section 7.1. That is, local coordinates along each of the three edges belonging to the higherorder face are introduced, see Fig. 28(b). They are related to the coordinates as
It is then simple to define as above in Eq. (7.1). The ramp function of Eq. (7.2) is adapted as
(7.6) 
Finally, the overall edge contribution for Eq. (7.4) is defined as
(7.7) 
It remains to define the face contribution wherefore we need the local coordinates corresponding to the triangular face element, see also Fig. 28(b),
(7.8) 
One needs the higherorder shape functions and the linear shape functions of the triangular face. Together with the corresponding nodes and this leads to the map
(7.9) 
We need a bubble function which strictly lives inside the triangular face. Therefore, one needs to subtract the edge contributions from the face contribution. We have transformed coordinates to based on Eq. (7.8). Next, special points are generated based on and the linear shape functions of the triangle, i.e.
where , , are simply the coordinates of the corner nodes of the tetrahedron in the reference configuration as seen in Fig. 28(a). One may then evaluate the edge contributions (7.7) based on obtaining We are now ready to define the face contribution to Eq. (7.4) as
(7.10) 
with the ramp function
7.3 Mapping for prisms
The general situation where all faces of the prism are defined by (curved) higherorder elements is discussed in [52] but the implementation is rather tedious. Herein, only one of the faces of the prismatic element may be of higherorder resulting into two different situations: Case 1 arises from the situation where a tetrahedral background element is cut into a subtetrahedron and a subprism; the higherorder side of the prism is then a triangle, see Fig. 29(a). Case 2 results when the tetrahedron is split into two subprisms; the higherorder side of the prism is then a quadrilateral, see Fig. 29(b). The map for the first case is defined in a straightforward way using the original idea of [21, 22], which is detailed for the present situation in a previous work of the authors [18].
In the second case, we wish to avoid the general definition of the map and rather outline the procedure in an illustrative way: As we are only interested in placing inner element nodes, it is noted that one may cut the prism into slices according to the blue lines in Fig. 29(b). Then, each slice is a triangle with one higherorder side and the mapping defined in Section 7.1 applies. Doing so in each slice defines all nodes of the sought higherorder prism.


References
 [1] Abedian, A.; Parvizian, J.; Düster, A.; Khademyzadeh, H.; Rank, E.: Performance of different integration schmenes in facing discontinuities in the Finite Cell Method. International Journal of Computational Methods, 10, 1 – 24, 2013.
 [2] Abedian, A.; Parvizian, J.; Düster, A.; Rank, E.: The finite cell method for the J2 flow theory of plasticity. Finite Elem. Anal. Des., 69, 37 – 47, 2013.
 [3] Babuška, I.; Caloz, G.; Osborn, J.E.: Special finite element methods for a class of second order elliptic problems with rough coefficients. SIAM J. Numer. Anal., 31, 945 – 981, 1994.
 [4] Babuška, I.; Melenk, J.M.: The partition of unity method. Internat. J. Numer. Methods Engrg., 40, 727 – 758, 1997.
 [5] Bathe, K.J.: Finite Element Procedures. PrenticeHall, Englewood Cliffs, NJ, 1996.
 [6] Belytschko, T.; Black, T.: Elastic crack growth in finite elements with minimal remeshing. Internat. J. Numer. Methods Engrg., 45, 601 – 620, 1999.
 [7] Belytschko, T.; Liu, W.K.; Moran, B.: Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, Chichester, 2000.
 [8] Burman, E.; Claus, S.; Hansbo, P.; Larson, M.G.; Massing, A.: CutFEM: Discretizing geometry and partial differential equations. Internat. J. Numer. Methods Engrg., 104, 472 – 501, 2015.
 [9] Burman, E.; Hansbo, P.: Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method. Comp. Methods Appl. Mech. Engrg., 199, 2680 – 2686, 2010.
 [10] Burman, E.; Hansbo, P.: Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. Applied Numerical Mathematics, 62, 328 – 341, 2012.
 [11] Cessenat, O.; Despres, B.: Application of an ultra weak variational formulation of elliptic PDEs to the twodimensional Helmholtz problem. SIAM J. Numer. Anal., 35, 255 – 299, 1998.
 [12] Chen, G.; Ohnishi, Y.; Ito, T.: Development of Highorder Manifold Method. Internat. J. Numer. Methods Engrg., 43, 685 – 712, 1998.
 [13] Cheng, K.W.; Fries, T.P.: Higherorder XFEM for curved strong and weak discontinuities. Internat. J. Numer. Methods Engrg., 82, 564 – 590, 2010.
 [14] Dréau, K.; Chevaugeon, N.; Moës, N.: Studied XFEM enrichment to handle material interfaces with higher order finite element. Comp. Methods Appl. Mech. Engrg., 199, 1922 – 1936, 2010.
 [15] Düster, A.; J. Parvizian, Z. Yang; Rank, E.: The finite cell method for threedimensional problems of solid mechanics. Comp. Methods Appl. Mech. Engrg., 197, 3768 – 3782, 2008.
 [16] Fries, T.P.: Higherorder accurate integration for cut elements with ChenBabuška nodes. In Advances in Discretization Methods: Discontinuities, virtual elements, fictitious domain methods. (Ventura, G.; Benvenuti, E., Eds.), Vol. 0, SEMA SIMAI Springer Series, Springer, Berlin, accepted, 2016.
 [17] Fries, T.P.; Belytschko, T.: The extended/generalized finite element method: An overview of the method and its applications. Internat. J. Numer. Methods Engrg., 84, 253 – 304, 2010.
 [18] Fries, T.P.; Omerović, S.: Higherorder accurate integration of implicit geometries. Internat. J. Numer. Methods Engrg., 106, 323 – 371, 2016.
 [19] Gockenbach, M.S.: Understanding and Implementing the Finite Element Method. SIAM, Philadelphia, PA, 2006.
 [20] Goldstein, C.I.: The weak element method applied to Helmholtz type equations. Appl. Numer. Math., 2, 409 – 426, 1986.
 [21] Gordon, W.J.; Hall, C.A.: Construction of curvilinear coordinate systems and applications to mesh generation. Internat. J. Numer. Methods Engrg., 7, 461 – 477, 1973.
 [22] Gordon, W.J.; Hall, C.A.: Transfinite element methods: blending function interpolation over arbitrary curved element domains. Numer. Math., 21, 109 – 129, 1973.
 [23] Hansbo, A.; Hansbo, P.: An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems. Comp. Methods Appl. Mech. Engrg., 191, 5537 – 5552, 2002.
 [24] Joulaian, M.; Hubrich, S.; Düster, A.: Numerical integration of discontinuities on arbitrary domains based on moment fitting. Comput. Mech., 57, 979 – 999, 2016.
 [25] Kudela, L.; Zander, N.; Kollmannsberger, S.; Rank, E.: Smart octrees: Accurately integrating discontinuous functions in 3D. Comp. Methods Appl. Mech. Engrg., 306, 406 – 426, 2016.
 [26] Laborde, P.; Pommier, J.; Renard, Y.; Salaün, M.: Highorder extended finite element method for cracked domains. Internat. J. Numer. Methods Engrg., 64, 354 – 381, 2005.
 [27] Legay, A.; Wang, H.W.; Belytschko, T.: Strong and weak arbitrary discontinuities in spectral finite elements. Internat. J. Numer. Methods Engrg., 64, 991 – 1008, 2005.
 [28] Legrain, G.; Chevaugeon, N.; Dréau, K.: High order XFEM and levelsets for complex microstructures: Uncoupling geometry and approximation. Computer Methods in Applied Mechanics and Engineering, 241244, 172 – 189, 2012.
 [29] Lehrenfeld, C.: High order unfitted finite element methods on level set domains using isoparametric mappings. Comp. Methods Appl. Mech. Engrg., 300, 716 – 733, 2016.
 [30] Leveque, R.; Randall, J.; Li, Z.: Immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J. Numer. Anal., 31, 1019 – 1044, 1994.
 [31] Li, Z.: A fast iterative algorithm for elliptic interface problems. SIAM J. Numer. Anal., 35, 230 – 254, 1998.
 [32] Löhner, R.; Cebral, J.R.; Camelli, F.F.; Baum, J.D.; Mestreau, E.L.; Soto, O.A.: Adaptive embedded/immersed unstructured grid techniques. Archives Of Computational Methods In Engineering, 14, 279 – 301, 2007.
 [33] Melenk, J.M.; Babuška, I.: The partition of unity finite element method: basic theory and applications. Comp. Methods Appl. Mech. Engrg., 139, 289 – 314, 1996.
 [34] Moës, N.; Dolbow, J.; Belytschko, T.: A finite element method for crack growth without remeshing. Internat. J. Numer. Methods Engrg., 46, 131 – 150, 1999.
 [35] Monk, P.; Wang, DaQing: A leastsquares method for the Helmholtz equation. Comp. Methods Appl. Mech. Engrg., 175, 121 – 136, 1999.
 [36] Mote, C.D.: Globallocal finite element. Internat. J. Numer. Methods Engrg., 3, 565 – 574, 1971.
 [37] Moumnassi, M.; Belouettar, S.; Béchet, É.; Bordas, S.P.A.; Quoirin, D.; PotierFerry, M.: Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces. Comp. Methods Appl. Mech. Engrg., 200, 774 – 796, 2011.
 [38] Mousavi, S.E.; Sukumar, N.: Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons. Comput. Mech., 47, 535 – 554, 2011.
 [39] Müller, B.; Kummer, F.; Oberlack, M.: Highly accurate surface and volume integration on implicit domains by means of momentfitting. Internat. J. Numer. Methods Engrg., 96, 512 – 528, 2013.
 [40] Neittaanmäki, P.; Tiba, D.: An embedding of domains approach in free boundary problems and optimal design. SIAM J. Control Optim., 33, 1587 – 1602, 1995.
 [41] Noor, A.K.: Globallocal methodologies and their applications to nonlinear analysis. Finite Elem. Anal. Des., 2, 333 – 346, 1986.
 [42] Osher, S.; Fedkiw, R.P.: Level Set Methods and Dynamic Implicit Surfaces. Springer, Berlin, 2003.
 [43] Parvizian, J.; Düster, A.; Rank, E.: Finite cell method: h and pextension for embedded domain problems in solid mechanics. Comput. Mech., 41, 121 – 133, 2007.
 [44] Rose, M.E.: Weak element approximations to elliptic differential equations. Numer. Math., 24, 185 – 204, 1975.
 [45] Saiki, E.M.; Biringen, S.: Numerical Simulation of a Cylinder in Uniform Flow: Application of a Virtual Boundary Method. J. Comput. Phys., 123, 450 – 465, 1996.
 [46] SalaLardies, E.; FernándezMéndez, S.; Huerta, A.: Optimally convergent highorder XFEM for problems with voids and inclusion. Proceedings of the ECCOMAS 2012, Vienna, Austria, 2012.
 [47] Schillinger, D.; Düster, A.; Rank, E.: The hpdadaptive finite cell method for geometrically nonlinear problems of solid mechanics. Internat. J. Numer. Methods Engrg., 89, 1171 – 1202, 2012.
 [48] Schillinger, D.; Ruess, M.: The Finite Cell Method: A Review in the Context of HigherOrder Structural Analysis of CAD and ImageBased Geometric Models. Archive Comp. Mech. Engrg., 22, 391 – 455, 2015.
 [49] Scott, R.L.: Finite Element Techniques for Curved Boundaries. PhDDissertation, Massachusetts Institute of Technology, Cambridge, MA, USA, 1973.
 [50] Sethian, J.A.: Level Set Methods and Fast Marching Methods. Cambridge University Press, Cambridge, 2 edition, 1999.
 [51] Shi, G.: Manifold method of material analysis. Trans. 9th Army Conf. on Applied Mathematics and Computing, Minneapolis, Minnesota, 51 – 76, 1992.
 [52] Šolín, P.; Segeth, K.; Doležel, I.: Higherorder finite element methods. CRC Press, Boca Raton, FL, 2003.
 [53] Stavrev, A.; Nguyen, L.H.; Shen, R.; Varduhn, V.; Behr, M.; Elgeti, S.; Schillinger, D.: Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method. Comp. Methods Appl. Mech. Engrg., 310, 646 – 673, 2016.
 [54] Stazi, F.L.; Budyn, E.; Chessa, J.; Belytschko, T.: An extended finite element method with higherorder elements for curved cracks. Comput. Mech., 31, 38 – 48, 2003.
 [55] Strouboulis, T.; Babuška, I.; Copps, K.: The design and analysis of the generalized finite element method. Comp. Methods Appl. Mech. Engrg., 181, 43 – 69, 2000.
 [56] Strouboulis, T.; Copps, K.; Babuška, I.: The generalized finite element method: an example of its implementation and illustration of its performance. Internat. J. Numer. Methods Engrg., 47, 1401 – 1417, 2000.
 [57] Sukumar, N.; Chopp, D.L.; Moës, N.; Belytschko, T.: Modeling holes and inclusions by level sets in the extended finiteelement method. Comp. Methods Appl. Mech. Engrg., 190, 6183 – 6200, 2001.
 [58] Szabó, B.; Düster, A.; Rank, E.: The pVersion of the Finite Element Method, Chapter 5, 119 – 139. John Wiley & Sons, Chichester, 2004.
 [59] Uzgoren, E.; J. Sim, J.; Shyy, W.: Markerbased, 3D adaptive cartesian grid method for multiphase flow around irregular geometries. Comput. Phys. Comm., 5, 1 – 41, 2009.
 [60] Ventura, G.: On the elimination of quadrature subcells for discontinuous functions in the eXtended finite element method. Internat. J. Numer. Methods Engrg., 66, 761 – 795, 2006.
 [61] Ventura, G.; Gracie, R.; Belytschko, T.: Fast integration and weight function blending in the extended finite element method. Internat. J. Numer. Methods Engrg., 77, 1 – 29, 2009.
 [62] Ye, T.; Mittal, R.; Udaykumar, H.S.; Shyy, W.: An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries. J. Comput. Phys., 156, 209 – 240, 1999.
 [63] Zi, G.; Belytschko, T.: New cracktip elements for XFEM and applications to cohesive cracks. Internat. J. Numer. Methods Engrg., 57, 2221 – 2240, 2003.
 [64] Zienkiewicz, O.C.; Taylor, R.L.: The Finite Element Method, Vol. 1 – 3. ButterworthHeinemann, Oxford, 2000.