Higher-order meshing of implicit geometries - part I: Integration and interpolation in cut elements

TL;DR

This paper introduces a higher-order meshing strategy for implicit geometries using level-set methods, enabling accurate integration and interpolation in cut elements without restrictions on background mesh structure.

Contribution

It presents an automatic meshing approach for cut elements with higher-order interface meshing and sub-element decomposition for precise integration.

Findings

Effective higher-order interface meshing in 2D and 3D.

Accurate integration within cut elements.

Handling of corners and edges in implicit geometries.

Abstract

An accurate implicit description of geometries is enabled by the level-set method. Level-set data is given at the nodes of a higher-order background mesh and the interpolated zero-level sets imply boundaries of the domain or interfaces within. The higher-order accurate integration of elements cut by the zero-level sets is described. The proposed strategy relies on an automatic meshing of the cut elements. Firstly, the zero-level sets are identified and meshed by higher-order interface elements. Secondly, the cut elements are decomposed into conforming sub-elements on the two sides of the zero-level sets. Any quadrature rule may then be employed within the sub-elements. 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.

Equations

$${\mathbf{r}}^{i+1}={\mathbf{r}}^i-\frac{{\varphi }^h\left({\mathbf{r}}^i\right)}{\nabla \varphi \left({\mathbf{r}}^i\right)\cdot \mathbf{N}}\cdot \mathbf{N}.$$

$${\varphi }_1\left(\mathbf{x}\right)$$

$${\varphi }_2\left(\mathbf{x}\right)$$

$$f\left(\mathbf{x}\right)=\text{nicefrac}12\cdot x+\text{nicefrac}14\cdot y+x^2+2\cdot y^3.$$

$${\epsilon }_1^{\Gamma }$$

$${\epsilon }_{\varphi }^{\Gamma }$$

$${\epsilon }_f^{\Gamma }$$

$${\epsilon }_{f_{1D}^h}^{\Gamma }$$

$${\epsilon }_{f_{2D}^h}^{\Gamma }$$

$${\varphi }_1\left(\mathbf{x}\right)$$

$${\varphi }_2\left(\mathbf{x}\right)$$

$$f\left(\mathbf{x}\right)=x^2+y^2+\text{nicefrac}12\cdot \cos (z)$$

$${\epsilon }_{f_{2D}^h}^{\Gamma }$$

$${\epsilon }_{f_{3D}^h}^{\Gamma }$$

$${\epsilon }_1^{\Omega }$$

$${\epsilon }_f^{\Omega }$$

$${\epsilon }_{f^h}^{\Omega }$$

$$N_1^{\star }\left(\mathbf{a}\right)=1-a-b,N_2^{\star }\left(\mathbf{a}\right)=a,N_3^{\star }\left(\mathbf{a}\right)=b$$

$$N_1^{\star }\left(\mathbf{a}\right)$$

$$N_3^{\star }\left(\mathbf{a}\right)$$

$$u_1=2a-1,u_2=b-a,u_3=1-2b$$

$$u_1=a,u_2=b,u_3=-a,u_4=-b.$$

$${\mathbf{r}}^{\mathrm{edge}k}\left(u_k\right)=\sum _{i=1}^{n_p}{N}_i^{\mathrm{HO}}\left(u_k\right)\cdot {\mathbf{r}}_i^{\mathrm{edge}k}-N_1^{\mathrm{lin}}\left(u_k\right)\cdot {\mathbf{r}}_1^{\mathrm{edge}k}-N_2^{\mathrm{lin}}\left(u_k\right)\cdot {\mathbf{r}}_{n_p}^{\mathrm{edge}k}$$

$$R_1=\frac{N_1^{\star }\left(\mathbf{a}\right)\cdot N_2^{\star }\left(\mathbf{a}\right)}{N_1^{\mathrm{lin}}(u_1)\cdot N_2^{\mathrm{lin}}(u_1)},R_2=\frac{N_2^{\star }\left(\mathbf{a}\right)\cdot N_3^{\star }\left(\mathbf{a}\right)}{N_1^{\mathrm{lin}}(u_2)\cdot N_2^{\mathrm{lin}}(u_2)},R_3=\frac{N_3^{\star }\left(\mathbf{a}\right)\cdot N_1^{\star }\left(\mathbf{a}\right)}{N_1^{\mathrm{lin}}(u_3)\cdot N_2^{\mathrm{lin}}(u_3)}$$

$$R_1=N_1^{\star }\left(\mathbf{a}\right)+N_2^{\star }\left(\mathbf{a}\right),R_2=N_2^{\star }+N_3^{\star },R_3=N_3^{\star }+N_4^{\star },R_4=N_4^{\star }+N_1^{\star }.$$

$$\mathbf{r}\left(\mathbf{a}\right)=\sum _{i=1}^{n_e}{N}_i^{\star }\left(\mathbf{a}\right)\cdot {\mathbf{r}}_i^{\star }+\sum _{k=1}^{n_e}{R}_k\cdot {\mathbf{r}}^{\mathrm{edge}k}\left(u_k(\mathbf{a})\right).$$

$$\mathbf{r}\left(\mathbf{a}\right)={\mathbf{r}}_{\mathrm{vertex}}\left(\mathbf{a}\right)+{\mathbf{r}}_{\mathrm{edge}}\left(\mathbf{a}\right)+{\mathbf{r}}_{\mathrm{face}}\left(\mathbf{a}\right).$$

$${\mathbf{r}}_{\mathrm{vertex}}\left(\mathbf{a}\right)=\sum _{i=1}^4{N}_i^{♢}\left(\mathbf{a}\right)\cdot {\mathbf{r}}_i^{\mathrm{tetra}}$$

$$N_1^{♢}\left(\mathbf{a}\right)=1-a-b-c,N_2^{♢}\left(\mathbf{a}\right)=a,N_3^{♢}\left(\mathbf{a}\right)=b,N_4^{♢}\left(\mathbf{a}\right)=c.$$

$$w_1=b-a,w_2=c-b,w_3=a-c.$$

$$R_1=\frac{N_2^{♢}\left(\mathbf{a}\right)\cdot N_3^{♢}\left(\mathbf{a}\right)}{N_1^{\mathrm{lin}}(w_1)\cdot N_2^{\mathrm{lin}}(w_1)},R_2=\frac{N_3^{♢}\left(\mathbf{a}\right)\cdot N_4^{♢}\left(\mathbf{a}\right)}{N_1^{\mathrm{lin}}(w_2)\cdot N_2^{\mathrm{lin}}(w_2)},R_3=\frac{N_4^{♢}\left(\mathbf{a}\right)\cdot N_2^{♢}\left(\mathbf{a}\right)}{N_1^{\mathrm{lin}}(w_3)\cdot N_2^{\mathrm{lin}}(w_3)}.$$

$${\mathbf{r}}_{\mathrm{edge}}\left(\mathbf{a}\right)=\sum _{k=1}^3{R}_k\cdot {\mathbf{r}}^{\mathrm{edge}k}\left(w_k(\mathbf{a})\right).$$

$$u\left(\mathbf{a}\right)=b+\text{nicefrac}13(1-a-b-c),v\left(\mathbf{a}\right)=c+\text{nicefrac}13(1-a-b-c).$$

$${\mathbf{r}}^{face234}\left(\mathbf{u}\right)=\sum _{i=1}^{n_q}{N}_i^{\mathrm{HO}}\left(\mathbf{u}\right)\cdot {\mathbf{r}}_i^{\mathrm{face}}-\sum _{i=1}^3{N}_i^{\mathrm{lin}}\left(\mathbf{u}\right)\cdot {\mathbf{r}}_{i+1}^{\mathrm{tetra}}.$$

$${\mathbf{a}}^{\star }\left(\mathbf{u}\left(\mathbf{a}\right)\right)=\sum _{i=1}^3{N}_i^{\mathrm{lin}}\left(\mathbf{u}\right)\cdot {\mathbf{a}}_{i+1}^{\mathrm{tetra}}$$

Full text

Keywords: Numerical integration, level-set method, fictitious domain method, XFEM, GFEM, interface capturing

Institute of Structural Analysis

Graz University of Technology

Lessingstr. 25/II, 8010 Graz, Austria

www.ifb.tugraz.at

fries@tugraz.at

1. 1 Introduction

2. 2 Background mesh and level-set data

3. 3 Reconstruction

4. 3.1 Valid level-set data and recursion

5. 3.2 Topological cases

6. 3.3 Root search and interface elements

7. 3.3.1 Start values and search directions in 2D

8. 3.3.2 Start values and search directions in 3D

9. 3.3.3 Iterative procedure

10. 3.4 Numerical results for reconstructions in 2D and 3D

11. 3.4.1 Numerical results for zero-isolines in 2D

12. 3.4.2 Numerical results for zero-isosurfaces in 3D

13. 4 Decomposition

14. 4.1 Decomposition in 2D

15. 4.2 Decomposition in 3D

16. 4.3 Map of integration points

17. 4.4 Numerical results for decompositions in 2D and 3D

18. 4.4.1 Preliminary examples

19. 4.4.2 Numerical results for decompositions in 2D

20. 4.4.3 Numerical results for decompositions in 3D

21. 5 Corners and edges

22. 6 Conclusions

23. 7 Appendix

24. 7.1 Mappings for triangles and quadrilaterals

25. 7.2 Mapping for tetrahedra

26. 7.3 Mapping for prisms

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 higher-order 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 $p$-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 higher-order 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 higher-order 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). Non-standard 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 inner-element 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 non-smooth solutions involve the weak element method [20, 44], the ultra weak variational method [11], the least-squares method as proposed in [35], and the global-local FEM as described in [36, 41]. Finally, we mention the manifold method [12, 51] which may be used for inner-element discontinuities.

A shared property of FDMs and XFEM-related 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 sub-elements or not. For the first class, the standard approach is to recursively decompose a cut element into polygonal sub-cells until the desired accuracy is obtained [1, 37, 14]. This typically leads to a very large number of integration points in the context of higher-order approximations [54, 63, 28, 26, 14]. It was already noted in [27, 13, 18, 46] that a decomposition into sub-elements with curved, higher-order edges or faces is a strategy which consistently enables the generation of higher-order 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, higher-order accuracy, general integrands, and to the presence of corners and edges often poses problems in this class of approaches.

The higher-order accurate numerical integration in three-dimensional, 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 level-set 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, higher-order sub-elements. Let us assume the abstract domain represented by Fig. 2(a) featuring an internal interface $\Gamma_{0}$ implicitly defined by the level-set method. It does not matter here whether $\Gamma_{0}$ is an interface between two materials in $\Omega^{+}$ and $\Omega^{-}$ or a boundary if e.g. $\Omega^{-}$ is considered to be a void region. A higher-order background mesh is introduced as shown in Fig. 2(b) and the level-set data is only given at the nodes. The implied zero-level 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 higher-order accuracy. Therefore, the zero-level sets are approximated by higher-order interface elements (reconstruction) and higher-order sub-elements 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 sub-elements and used in FDMs or XFEM-related 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 level-set 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 zero-level sets, i.e. on manifolds. Therefore, only the mesh of reconstructed interface elements will be needed. In the third part, BVPs with inner-element boundaries and interfaces are considered using meshes consisting of regular and cut background elements, the latter decomposed into conforming sub-elements as described herein. Further studies will then outline how fully automatic, higher-order 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 higher-order background meshes used herein and the level-set method. In Section 3, the identification and meshing of the zero-level set is described in detail (reconstruction). The decomposition into higher-order, conforming sub-elements is detailed in Section 4. The proper consideration of corners and edges with several level-set 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 level-set data

A two or three-dimensional domain of interest, $\Omega\in\mathbb{R}^{d}$, $d=\{2,3\}$ is defined implicitly based on the level-set method [42, 50]. That is, a scalar function $\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$, $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\in\Omega$ is given, called level-set function, which is positive on one side of the boundary/interface and negative on the other. The function is continuous and its zero-level set defines the position of a smooth boundary or interface. Note that $\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ 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 level-set functions for the description.

The domain is fully immersed in a background mesh composed by higher-order 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 level-set data is only given at the nodes of the background mesh, $\phi_{i}=\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}\right)$, and is interpolated by standard finite element shape functions $N_{i}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ in-between, hence, $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)=\sum_{i}N_{i}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)\cdot\phi_{i}$. That is, once the nodal level-set data is given, no (further) use of any analytic knowledge of the level-set functions is made. Obviously, the background mesh enables a higher-order 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 higher-order accurate integration in the implicit geometry. Therefore, the background elements are decomposed into sub-elements 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 sub-element. The decomposition is successively applied in each element for the involved level-set functions. This is a two-step procedure for each zero-level set that cuts an element: (1) Reconstruction, which is the meshing of the zero-level set with higher-order interface elements. (2) Decomposition of cut elements into higher-order sub-elements 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 three-dimensional 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 level-set functions. Thereafter, (3) the sub-elements are mapped to the corresponding physical background element and (4) any desired quadrature points are mapped to these sub-elements, 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 zero-level set of a level-set function. The nodal values of the level-set function at the element nodes are taken from the physical element, i.e. $\phi(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i})=\phi(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i})$. The task is to approximate the zero-level set of $\phi^{h}(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}})$ by means of a higher-order interface element. That is, roots have to be identified on the zero-level set which are then the element nodes of the interface element. In-between, the zero-level 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 root-finding 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 level-set data and recursion

It is obvious that higher-order elements may feature very complicated topologies of zero-level sets. For example, element edges may be cut several times, or the zero-level 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 level-set 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 level-set data at the sample points is considered for these checks.

If the level-set 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 level-set 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 sub-element after the decomposition features a negative Jacobian, a recursive refinement is still triggered although the initial level-set data was possibly valid in the above sense. This may, for example, be the case for strongly curved zero-level sets, see Fig. 6(d). A level-set function must not be zero right at the corner nodes of the background element; if that is the case, the level-set value is simply perturbed in the range of $10^{-13}$ without any noticeable change in the results.

3.2 Topological cases

From now on, it is assumed that valid level-set data is present in a (refined) reference element. Then, the zero-level 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 sub-triangle and one sub-quadrilateral; the latter may be further decomposed into sub-triangles. A cut tetrahedron is decomposed either into one sub-tetrahedron and one sub-prism (topology 1), see Fig. 8(b), or into two sub-prisms (topology 2), see Fig. 8(c). Obviously, the sub-prisms may be further decomposed into tetrahedra. There is no need to avoid hanging nodes in the decomposition into sub-elements 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 level-set function $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\right)$ which is obtained by solving a non-linear 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 zero-level 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 zero-level 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). 2. 2.

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

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

In direction of the gradient of the level-set function, i.e. $\nabla\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\right)$, see Fig. 9(d) and (h).

For later reference in the numerical studies, the 2D variants are characterized by two-digit numbers: The first digit is either $1$ for a linear or $2$ 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 higher-order 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 zero-isosurface based on the mapping defined in [52]. This mapping defines a surface implied by the three or four curved edges (higher-order line elements) of the triangular or quadrilateral zero-isosurface within the cut tetrahedron, see Fig. 10. The definition of the mapping is outlined in the appendix 7.1 of this work. Next, a higher-order 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 level-set function.

It is thus seen that also in 3D, the 2D-variants 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 zero-level 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 $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\right)$ 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 zero-level set.

Mathematically, the algorithm is described as follows. The task is to find the root of $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\right)$ in the reference element, i.e. some position on the zero-level set. The starting point of the iterative procedure is labeled $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{0}$. The Newton-Raphson-type algorithm for all approaches considered here is based on the following iteration:

[TABLE]

$\mathchoice{\displaystyle\boldsymbol{N}}{\textstyle\boldsymbol{N}}{\scriptstyle\boldsymbol{N}}{\scriptscriptstyle\boldsymbol{N}}$ 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 $\mathchoice{\displaystyle\boldsymbol{N}}{\textstyle\boldsymbol{N}}{\scriptstyle\boldsymbol{N}}{\scriptscriptstyle\boldsymbol{N}}=\nabla\phi\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{0}\right)$. Otherwise, it may also be the current gradient $\mathchoice{\displaystyle\boldsymbol{N}}{\textstyle\boldsymbol{N}}{\scriptstyle\boldsymbol{N}}{\scriptscriptstyle\boldsymbol{N}}=\nabla\phi\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{i}\right)$.

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 zero-level 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 sub-optimal.

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 sub-elements 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 zero-isolines in 2D

We start with the two-dimensional situation. Following [16], two different level-set functions are considered:

[TABLE]

with $r=0.7123$, $R\left(\theta\right)=0.5+0.1\cdot\sin\left(8\theta\right)$ and $\theta\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)=\mathrm{atan}\left(y/x\right)$. See Fig. 12(a) for a visualization of the zero-level sets of $\phi_{1}$ and $\phi_{2}$. The zero-level set of $\phi_{1}$ is a circle with radius $r$ and is frequently used in the literature, e.g. [18], and $\phi_{2}$ is similar to [31]. The function $f\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ is defined as

[TABLE]

and is later integrated on the zero-level sets.

In the convergence studies, $\left\{6,10,20,30,50,70,100,150,200,300\right\}$ elements are employed per dimension. We use standard Gauss rules for the integration on the zero-level set $\Gamma_{0}$ with a rather high order of $11$ 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 $f\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ on $\Gamma_{0}$ [16]:

[TABLE]

where $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}$ are the integration points in the interface elements with integration weights $w_{i}$. It is noted that, in contrast to the nodes of the interface elements, the integration points at $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}$ are only approximately on the zero-level set of $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$.

The first error, $\varepsilon_{1}^{\Gamma}$ 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 $\phi_{1}$, this is the circumference of a circle). The second error, $\varepsilon_{\phi}^{\Gamma}$, integrates the level-set function on the zero-level set, which, in the ideal case, would be zero. Nevertheless, because the nodes of the interface elements are on the zero-level set of $\phi^{h}$ rather than $\phi$ and, furthermore, the integration points are only approximately on the zero-level set of $\phi^{h}$, $\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}\right)$ is not exactly zero. The other error norms integrate either the exact or the interpolated function $f\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ as given in Eq. (3.4). For $\varepsilon_{f}^{\Gamma}$, the function $f\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ is evaluated exactly at the integration points $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}$. Instead, $\varepsilon_{f_{1D}^{h}}^{\Gamma}$ 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, $\varepsilon_{f_{2D}^{h}}^{\Gamma}$ 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 $14$ and $24$, where the search directions are based on the gradient of the level-set function are split into two sub-versions 14a/14b and 24a/24b depending on whether the search path is constant, $\nabla\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{0}\right)$, or changes during the iteration, $\nabla\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{i}\right)$. A few variants actually underperform, in particular variant $11$ and $12$. 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 level-set 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 $1$) compared to the Hermite reconstruction (first digit is $2$). Because it is found that this also applies to the three-dimensional 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, $\phi_{1}$ is chosen and an order of $5$ 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 level-set 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 $\phi_{1}$ and $\phi_{2}$, 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 $h$ and complicates the proper placement of the result in the convergence plots).

3.4.2 Numerical results for zero-isosurfaces in 3D

Numerical results are now presented for the integration on zero-isosurfaces $\Gamma_{0}$ 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 $13$ 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 $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{0}$ and the second based on the gradient of the level-set function at $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{0}$. Variant C is based on the gradient as well, but depending on the points $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{i}$ during the iterative procedure.

The following two level-set functions are considered in 3D:

[TABLE]

with $r=0.7123$. The zero-isosurfaces of these two level-set functions are visualized in Fig. 14. There, also a coarse tetrahedral background mesh is shown. For the integrand, the function

[TABLE]

is chosen. In the convergence studies, $\left\{6,10,14,20,30,50,70,100\right\}$ 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

[TABLE]

Therein, $f_{2D}^{h}$ is the integrand interpolated by the shape functions of the reconstructed interface elements and $f_{3D}^{h}$ 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 zero-isosurfaces of $\phi_{1}$ and $\phi_{2}$.

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 zero-level set of the interpolated level-set $\phi^{h}$ which, in fact, is only an approximation of $\phi$. In some error norms, only the interpolated function $f^{h}$ is integrated rather than $f$. Finally, also the numerical integration based on Gauss quadrature introduces an integration error. Nevertheless, it was shown that $f\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ is integrated and interpolated optimally on zero-level sets implied by $\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$.

4 Decomposition

Once the interface elements are defined, i.e. the zero-level set is meshed with higher-order interface elements, the cut background elements are to be decomposed into sub-elements depending on the topology of the cut situation, see Section 3.2. This is outlined in Fig. 4 for the two- and three-dimensional situation, respectively. The resulting sub-elements share the property that they feature one higher-order 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 sub-elements on the two sides of the interface.

4.1 Decomposition in 2D

For the two-dimensional case, a cut reference background element in the coordinate system $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}$ falls into a triangular and a quadrilateral sub-element with one higher-order side, see Fig. 4(a). In order to define the element nodes of the sub-elements, a mapping from a higher-order reference triangle or quadrilateral in coordinates $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}$ to each sub-element in coordinates $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}$ 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 sub-elements.

This is outlined for the case of a triangular sub-element where we introduce three different mappings $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\left(\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}\right)$. That is, the reference element in $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}$ is mapped differently to the coordinate system $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}$ 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 $\mathbb{R}^{3}$. The situation here is the reduced case where the mapping is to $\mathbb{R}^{2}$ 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 $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{\mathrm{edge}\,k}\left(u_{k}\right)\neq\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}$ only for $k=2$.

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 $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{\mathrm{edge}\,k}\left(u_{k}\right)=\mathchoice{\displaystyle\boldsymbol{0}}{\textstyle\boldsymbol{0}}{\scriptstyle\boldsymbol{0}}{\scriptscriptstyle\boldsymbol{0}}$ and $R_{k}=0$ for $k=1,3$. For $k=2$, we define $u_{2}=(b-a)/(a+b)$, $R_{2}=(a+b)$, and $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{\mathrm{edge}\,2}\left(u_{2}\right)$ 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 $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}_{i}$. These points are mapped to the sub-element in $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}$ 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 $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}$.

Although Fig. 16(c) indicates that all three mappings successfully map to the linear triangle with one higher-order side, it is found that the second map yields sub-optimal accuracy in the convergence studies. This could be traced back to the fact that the mapping is not as smooth as needed for higher-order 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 higher-order side which explains its sub-optimal 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 three-dimensional 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 sub-elements as shall be seen below.

The procedure for quadrilaterals with one higher order side, is again an adaption of the more general map outlined in the appendix 7.1. In this case, the mapping coincides with the blending function mapping of [21, 22, 58].

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 sub-tetrahedron and a sub-prism 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 sub-prism, see Fig. 8(c). The curved faces coincide with the reconstructed higher-order interface elements as discussed in Section 3. Each sub-element has straight edges except for those belonging to the curved face.

Again, the task is to define maps (of element nodes) from higher-order reference tetrahedra and prisms in coordinates $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}$ to the cut reference background element in $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}$. For some possible situations, examples are shown in Fig. 18 for all types of sub-elements. For the sub-elements with one higher-order 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 sub-elements 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 sub-elements. Depending on the particular context for which higher-order accurate integration points are sought, this may be preferred in the sub-elements 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 XFEM-related 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 sub-elements for complex level-set data in 2D. According to Section 3.1, recursive refinements are necessary to obtain valid level-set data in the refined elements. Note the relation to Fig. 6. Examples for decompositions in 3D are shown in Fig. 20 which correspond to zero-level 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, level-set functions and integrands are considered here. Three different (relative) error norms are introduced for the integration in $\Omega^{-}$:

[TABLE]

where $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}$ are 2D integration points in the special sub-elements and $w_{i}$ 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 zero-level 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 level-set 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 sub-optimal second variant of the mapping to the sub-elements described in Section 4.1: There, it turned out that results are sub-optimal in $\varepsilon_{f^{h}}^{\Omega}$ (convergence rates are bounded by 4). Optimal results for $\varepsilon_{1}^{\Omega}$ and $\varepsilon_{f}^{\Omega}$ 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 $\phi_{1}$ and $\phi_{2}$ 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 $h$ 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 $\phi_{1}$ and $\phi_{2}$ 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, sub-optimal results for reconstructions can not be better for decompositions. Of course, optimal reconstructions may still lead to sub-optimal decompositions e.g. if the maps to the sub-elements, 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 level-set function. It is possible to define a level-set function $\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ whose zero-level set also involves corners and edges. Then, this function is necessarily $C_{0}$-continuous. However, the meshing procedures described in the previous sections rely on an interpolated level-set function $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$. Provided that $\phi\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ is only $C_{0}$-continuous, an accurate interpolation $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ 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 $\phi^{h}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)$ and a higher-order decomposition would not be possible.

Therefore, we suggest to use multiple level-set functions for the domain description. Then, several level-set values are present at the nodes of the background mesh, one for each smooth part of the interface. In two dimensions, two level-set 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 zero-level sets $\Gamma_{1,0}$ and $\Gamma_{2,0}$ of the level-set functions $\phi_{1}$ and $\phi_{2}$, respectively. Assume the level-set functions are negative inside the circular zero-level sets. It is seen in Fig. 23(c) how the two level-set functions imply 4 different sub-domains based on their signs. The gray region in Fig. 23(a) is $\Omega^{\star}=\left\{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}:\phi_{1}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)>0\;\mathrm{and}\;\phi_{2}\left(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\right)>0\right\}$. It is thus seen, how two level-set functions are able to define inner-element corners in 2D. No additional requirements on the background mesh are needed to capture the corners accurately.

In three dimensions, two level-set functions are able to define an edge and three functions a corner. See Fig. 24 for a graphical representation. Three zero-level sets of $\phi_{1}$, $\phi_{2}$, and $\phi_{3}$ are shown in Fig. 24(a). It is seen how the intersections of any two level-set functions implicitly define edges and how three level-set functions imply corners. Of course, based on the signs of the three level-set functions, one may distinguish $2^{3}=8$ sub-regions of the background mesh. Depending on the application, one may now want to compose the domain based on these sub-regions. Two different examples are shown in Figs. 24(b) and (c); both feature implicit edges and corners.

It is now clear how multiple level-set 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 sub-regions based on the signs of the involved level-set 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 sub-regions.

Fortunately, the decomposition of background elements with respect to several level-set 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 level-set function $\phi_{1}$. Then, the remaining level-set functions are interpolated at the new element nodes of the sub-elements of the cut background element. The decomposition with respect to the next level-set function is carried out just as if the already decomposed mesh were the initial mesh. Therefore, it may be useful to further decompose quadrilateral sub-elements in 2D into triangles and prismatic sub-elements in 3D into tetrahedra, so that the newly generated mesh (conforming to the zero-level sets of previous $\phi_{i}$) only consists of the same element type. Fig. 25 shows an example where a reference element is decomposed with respect to $4$ level set functions. The resulting decompositions for each new level-set functions are shown and the different sub-regions that may be identified by the signs of the level-set functions are color-coded. Fig. 25(f) shows a sub-domain 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 level-set functions is a truly straightforward successive application of the meshing procedures described above. The ability to naturally consider an arbitrary number of level-set 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 XFEM-related methods. The boundaries of a domain of interest or interfaces inside the domain are then easily defined based on the level-set method. Using a higher-oder background mesh enables a straightforward path to a higher-order 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 higher-order 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 higher-order finite elements on the two sides of an interface which is called decomposition. Before, the zero-level 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 non-linear detection of the zero-level sets. For the decomposition, mappings of element nodes into the special sub-elements with one higher-order side are needed and they have to be sufficiently smooth. It is noted that depending on the level-set data inside the elements one may have to recursively refine background elements until valid level-set data is obtained. Recursive refinements are also needed when the (initial) reconstruction or decomposition fails e.g. due to negative Jacobians of the resulting sub-elements. When the curvatures of the level-set 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 sub-elements 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 XFEM-related 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 well-shaped. Nevertheless, it shall be seen that higher-order 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 sub-elements on the two sides of the zero-level sets, the outer contour is given by linear elements with one higher-order, typically curved side. For the placements of start values for the reconstruction in tetrahedra, three or four higher-order 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 $n_{e}=3$ (or $n_{e}=4$) 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 $n_{p}$ and corresponding shape functions is associated with each edge. The task is to define a smooth mapping for all points $\left(a,b\right)$ in a reference triangle (or quadrilateral) to the three-dimensional situation in $\left(r,s,t\right)$. 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 $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}^{\mathrm{edge}\,k}\in\mathbb{R}^{3}$ with $i=1,\dots,n_{p}$ and $k=1,\dots,n_{e}$. 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 $N_{k}^{\star}\left(\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}\right)$, more precisely, for triangles

[TABLE]

and for quadrilaterals

[TABLE]

The higher-order shape functions of the line elements are denoted as $N_{i}^{\mathrm{HO}}\left(u\right)$ and are defined in a one-dimensional reference element with coordinates $u\in[-1,1]$. Herein, these are standard Lagrange shape functions. Furthermore, the linear shape functions $N_{1}^{\mathrm{lin}}(u)=\nicefrac{{1}}{{2}}(1-u)$ and $N_{2}^{\mathrm{lin}}(u)=\nicefrac{{1}}{{2}}(1+u)$ are needed. Local coordinates $u_{k}\in\left[-1,1\right]$ 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 $\left(a,b\right)$. For triangles, this gives

[TABLE]

and for quadrilaterals

[TABLE]

Next, we define

[TABLE]

where the first term on the right hand side defines the curve in $\mathbb{R}^{3}$, and the other two subtract the linear interpolant. Furthermore, a ramp function is needed which is, for triangles,

[TABLE]

and, for quadrilaterals,

[TABLE]

The overall map $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\left(\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}\right)$, i.e. the shape of the resulting surface element which is implied by the curved line elements, is then defined as

[TABLE]

The first term on the right hand side is the linear interpolant of the corner points of the triangle (or quadrilateral) and $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}^{\star}$ 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 higher-order triangular element with $n_{q}$ nodes. The more general case where all four faces are curved is not needed herein. The curved face has three curved edges, which are higher-order line elements with $n_{p}$ 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 $\left(a,b,c\right)$ in a reference tetrahedron to the three-dimensional situation in $\left(r,s,t\right)$.

The nodal coordinates of the triangular face element are denoted as $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}^{\mathrm{face}}\in\mathbb{R}^{3}$ with $i=1,\dots,n_{q}$. The three corner nodes of the face are set as the nodes $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}^{\mathrm{tetra}}$ with $i=2,3,4$ of the tetrahedron. The position of the remaining node of the tetrahedron $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{1}^{\mathrm{tetra}}$ 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,

[TABLE]

The first part is simply the linear mapping of a tetrahdron, i.e.

[TABLE]

with the linear shape functions

[TABLE]

For the edge contributions, we proceed similar as in Section 7.1. That is, local coordinates $w_{k}\in\left[-1,1\right]$ along each of the three edges belonging to the higher-order face are introduced, see Fig. 28(b). They are related to the coordinates $\left(a,b,c\right)$ as

[TABLE]

It is then simple to define $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}^{\mathrm{edge}\,k}\left(w_{k}\right)$ as above in Eq. (7.1). The ramp function of Eq. (7.2) is adapted as

[TABLE]

Finally, the overall edge contribution for Eq. (7.4) is defined as

[TABLE]

It remains to define the face contribution wherefore we need the local coordinates $(u,v)$ corresponding to the triangular face element, see also Fig. 28(b),

[TABLE]

One needs the higher-order shape functions $N_{i}^{\mathrm{HO}}\left(u,v\right)$ and the linear shape functions $N_{i}^{\mathrm{lin}}\left(u,v\right)$ of the triangular face. Together with the corresponding nodes $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}^{\mathrm{face}}$ and $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{i}^{\mathrm{tetra}}$ this leads to the map

[TABLE]

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 $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}$ to $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}$ based on Eq. (7.8). Next, special points $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}^{\star}$ are generated based on $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}$ and the linear shape functions $N_{i}^{\mathrm{lin}}\left(\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}\right)$ of the triangle, i.e.

[TABLE]

where $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}_{k}^{\mathrm{tetra}}$, $k=2,3,4$, 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 $\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}^{\star}$ obtaining $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}_{\mathrm{edge}}\left(\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}^{\star}\right).$ We are now ready to define the face contribution to Eq. (7.4) as

[TABLE]

with the ramp function

[TABLE]

7.3 Mapping for prisms

The general situation where all faces of the prism are defined by (curved) higher-order elements is discussed in [52] but the implementation is rather tedious. Herein, only one of the faces of the prismatic element may be of higher-order resulting into two different situations: Case 1 arises from the situation where a tetrahedral background element is cut into a sub-tetrahedron and a sub-prism; the higher-order side of the prism is then a triangle, see Fig. 29(a). Case 2 results when the tetrahedron is split into two sub-prisms; the higher-order 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 $\mathchoice{\displaystyle\boldsymbol{r}}{\textstyle\boldsymbol{r}}{\scriptstyle\boldsymbol{r}}{\scriptscriptstyle\boldsymbol{r}}\left(\mathchoice{\displaystyle\boldsymbol{a}}{\textstyle\boldsymbol{a}}{\scriptstyle\boldsymbol{a}}{\scriptscriptstyle\boldsymbol{a}}\right)$ 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 higher-order side and the mapping defined in Section 7.1 applies. Doing so in each slice defines all nodes of the sought higher-order prism.