Numerical integration over implicitly defined domains for higher order unfitted finite element methodsthanks: Partially supported by NSF through the Division of Mathematical Sciences grants 1315993 and 1522252.

Maxim A. Olshanskii1    Danil Safin1
11Department of Mathematics, University of Houston, Houston, Texas 77204-3008 (molshan,dksafin)

The paper studies several approaches to numerical integration over a domain defined implicitly by an indicator function such as the level set function. The integration methods are based on subdivision, moment–fitting, local quasi-parametrization and Monte-Carlo techniques. As an application of these techniques, the paper addresses numerical solution of elliptic PDEs posed on domains and manifolds defined implicitly. A higher order unfitted finite element method (FEM) is assumed for the discretization. In such a method the underlying mesh is not fitted to the geometry, and hence the errors of numerical integration over curvilinear elements affect the accuracy of the finite element solution together with approximation errors. The paper studies the numerical complexity of the integration procedures and the performance of unfitted FEMs which employ these tools.

1 Introduction

Numerical approaches for solving PDEs that integrate the underlying geometric information, such as isogeometric analysis [14], are in the focus of research over the last decade. In the isogeometric analysis, the functions used for geometry representation are also employed to define functional spaces in the Galerkin method. Several other approaches try to make the grid generation and geometry description independent. Unfitted finite element methods, such as immersed boundary methods [16], extended FEM [4, 9], cut FEM [5], or trace FEM [19], use a sufficient regular background grid, but account for geometric details by modifying the spaces of test and trial functions or the right-hand side functional.

The accuracy of unfitted FEM depends on several factors. These are the approximation property of basic finite element space, the accuracy of underlying geometry recovering, and the error introduced by numerical integration. The present paper first reviews the analysis of two higher order finite element methods: one is an unfitted FEM for the Neumann problem in a bounded curvilinear domain, another one is a narrow-band FEM for an elliptic PDE posed on a closed smooth manifold. In both cases, the domain (a volume or a surface) is given implicitly by a discrete level-set function. Practical implementation of these methods (as well as many other unfitted FEM) leads to the following problem: Given a simplex , a smooth function defined on and a polynomial of degree such that on , evaluate the integral


The paper focuses on several numerical approaches to problem (1) in . In the context of numerical solution of PDEs, is typically a polynomial; is the order of geometry recovery.

We note that for one has to integrate over a polygon (or polyhedron). Then for a polynomial function an exact numerical integration is straightforward through subdividing the polygon (polyhedron) into a finite number of triangles (tetrahedra) and applying standard Gauss quadratures on each triangle. However, for the problem appears to be less trivial and building an exact quadrature rule for does not look feasible. This can be realised by considering the simple 1D example with and . Solving (1) becomes equivalent to finding the root of . The latter problem is resolved (only) in radicals for and is well known to have no general algebraic solution for by the Abel theorem.

The problems of building numerical quadratures for the implicitly defined volume integrals (1) and the implicitly defined surface integrals

However, for a general superposition of and the zero level set of , the smearing leads to significant integration errors which are hard to control. Another numerical integration technique is based on an approximation of by elementary shapes. Sub-triangulations or quadtree (octree) Cartesian meshes are commonly used for these purposes. On each elementary shape a standard quadrature rule is applied. The sub-triangulation is often adaptively refined towards the zero level of . The approach is popular in combination with higher order extended FEM for problems with interfaces, see, e.g., [1, 17, 8], and the level-set method [15, 13]. Although numerically stable, the numerical integration based on sub-partitioning may significantly increase the computational complexity of a higher order finite element method, since the number of function evaluations per scales with for some depending on the order of the FEM. In several recent papers [18, 21, 10] techniques for numerical integration over implicitly defined domains were devised that have optimal computational complexity. The moment–fitting method from  [18] uses polynomial divergence free basis of vector function to approximate the integrand and further reduce the volume integrals to boundary integrals to find the weights of a quadrature formula by a least-square fitting procedure. We recall the moment–fitting method in section 3. For the case when is a hyper-rectangle, the approach in [21] converts the implicitly given geometry into the graph of an implicitly defined height function. The approach leads to a recursive algorithm on the number of spatial dimensions which requires only one-dimensional root finding and one-dimensional Gaussian quadrature. In [10], a zero level-set of is approximated by higher order interface elements. These elements are extended inside so that is covered by regular and curvilinear simplexes. For curvilinear simplexes a mapping to the reference simplex is constructed. Further standard Gauss quadratures are applied.

In this paper, we develop an approach for (1) based on the local quasi-parametrization of the zero level set of . Similar to [21] the zero level set is treated as a graph of an implicitly given function. With the help of a 1D root finding procedure the integration over is reduced to the integration over regular triangles and the recursive application of 1D Gauss quadratures. The technique is also related to the method of local parametrization for higher order surface finite element method of [11]. We compare the developed method with several other approaches for the numerical integration in the context of solving partial differential equations in domains with curvilinear boundaries and over surfaces. For the comparison purpose we consider the method of sub-triangulation for , the moment–fitting method, and the Monte-Carlo method.

The rest of the paper is organized as follows. We first recall unfitted finite element method for solving an elliptic PDE in a domain with curvilinear boundary and an elliptic PDE posed on a surface. For the surface PDE we use the method from [20] of a regular extension to a narrow band around the surface. The error analysis of these unfitted FE methods is also reviewed. Further in section 3 we discuss the methods for numerical integration of (1), which further used to build the FEM stiffness matrices. Section 4 collects the result of numerical experiments. Section 5 concludes the paper with a few closing remarks.

2 Unfitted FEM

We assume that is an open bounded subset in , , with a boundary , which is a connected compact hypersurface in . In this paper we apply unfitted FE methods to elliptic equations posed in and on . As model problems, let us consider the Poisson and the Laplace–Beltrami problems:




with some strictly positive or , respectively.

2.1 Preliminaries

To define finite element methods, we need a formulation of the surface PDE (3) based on normal extension to a narrow band. First, we introduce some preliminaries. Denote by a domain consisting of all points within a distance from less than some :

Let be the signed distance function, for all . The surface is the zero level set of :

We may assume on the interior of and on the exterior. We define for all . Thus, is the normal vector on , and for all . The Hessian of is denoted by :

For , the non-zero eigenvalues of are the principal curvatures. Hence, one can choose such sufficiently small positive that is uniformly positive definite on . For denote by the closest point on . Assume that is sufficiently small such that the decomposition is unique for all . For a function on we define its extension to :

Thus, is the extension of along normals on .

We look for solving the following elliptic problem


with . The Neumann condition in (4) is the natural boundary condition. The following results about the well-posedness of (4) and its relation to the surface equations (3) have been proved in [20]:

(i) For , the problem (4) has the unique weak solution , which satisfies , with a constant depending only on and ;

(ii) For the solution to (4) the trace function is an element of and solves a weak formulation of the surface equation (3).

(iii) The solution to (4) satisfies . Using the notion of normal extension, this can be written as

(iv) Additionally assume , then and with a constant depending only on , and .

2.2 FEM formulations

Let , , be a polygonal (polyhedral) domain such that for problem and for problem . Assume we are given a family of regular triangulations of such that . For a triangle denote by the diameter of the inscribed circle. Denote


For the sake of presentation, we assume that triangulations of are quasi-uniform, i.e., is uniformly bounded in .

It is computationally convenient not to align (not to fit) the mesh to or . Thus, the computational domain approximates or and has a piecewise smooth boundary which is not fitted to the mesh .

Let be a continuous piecewise, with respect to , polynomial approximation of the surface distance function in the following sense:


with some (for problem (2) a generic level-set function can be considered, not necessary a signed distance function). Then one defines


with , . Note, that in some applications the surface may not be known explicitly and only a finite element approximation to the distance function is known. Otherwise, one may set , where is a suitable piecewise polynomial interpolation operator. Estimate (6) is reasonable, if is a polynomial of degree and . The latter is the case for -smooth .

The space of all continuous piecewise polynomial functions of a degree with respect to is our finite element space:


The finite element method for problem (2) reads: Find satisfying


where are suitable extensions of and to . For problem (3), the finite element method is based on the extended formulation (4) and consists of finding that satisfies


It should be clear that only those basis functions from contribute to the finite element formulations (9), (10) and are involved in computations that do not vanish everywhere on . Error estimates for the unfitted finite element methods (9), (10) depend on the geometry approximation and the order of finite elements. Let and assume solves the Neumann problem (2) and solves (9). Then it holds


where a constant is independent of , is the degree of the finite element polynomials, and the geometry approximation order defined in (6).

For the Neumann problem with and a compatibility condition on the data, the estimate in (11) is proved in [3] subject to an additional assumption on . In that paper it is assumed that and match on the edges of intersected by . This is not necessary the case for defined implicitly from the discrete level function . For implicitly defined domains, the convergence result in (11) follows from the more recent analysis in  [12], where a suitable mapping of on was constructed.

For the FE formulation (10) of the Laplace-Beltrami problem (3), approximations to and are required. If is given explicitly, one can compute and and set , and in (9). Otherwise, if the surface is known approximately as, for example, the zero level set of a finite element distance function , then, in general, and one has to define a discrete Hessian and also set . A discrete Hessian can be obtained from by a recovery method, see, e.g., [2, 23]. Assume that some is provided and denote by the approximation order for in the (scaled) -norm:


where denotes the area (volume) of .

The convergence of the finite element method (9) is summarized in the following result from [20]. Let , , and assume solves the surface problems (3) and solves (9). Then it holds


where a constant is independent of , and , , are the finite elements, Hessian recovery, and distance function approximation orders defined in (8), (6) and (12), respectively. Numerical experiments in [20] show that typically demonstrates a one order higher convergence rate than the norm of the error.

Note that both error estimates (11) and (13) assume exact numerical integration. In the introduction, we discussed that for the exact numerical integration is not feasible. The error of the numerical quadrature should be consistent with the finite element interpolation and geometric error to ensure that the FEM preserves the optimal accuracy.

Remark 1

For the purpose of improving algebraic properties of an unfitted finite element method a stabilization procedure was suggested in [6] for elliptic equations posed in volumetric domains and further extended to surface PDEs in [7]. The procedure consists of adding a special term penalizing the jump of the solution gradient over the edges (faces) of the triangles (tetrahedra) cut by .

Let ( is the set of all elements having non-empty intersection with the boundary of the numerical domain ). By denote the set of all edges in 2D or faces in 3D shared by any two elements from . Define the term

Here denotes the jump of the normal derivative of across ; are stabilization parameters (in our computations we set ). The edge-stabilized trace finite element reads: Find such that


for all . Here and are the bilinear forms and the right-hand side functional corresponding to the finite element methods (9) or (10).

For continuous bulk finite element methods on quasi-uniform regular tetrahedral meshes, the optimal orders of convergence for (14) were proved in [6]. In our numerical studies we tested the stabilized formulation (14) with higher order elements. We observed very similar convergence results for the formulations with and without term, including sub-optimal/irregular behaviour with moment-fitting quadratures and optimal with other integration techniques. For the systems of linear algebraic equations we use exact ‘backslash’ solves in either case. Therefore, we shall not report results for (14) in addition to finite element formulations (9) and (10).

3 Numerical integration

For a simplex , we are interested in computing integrals over with a certain accuracy , i.e. for a sufficiently smooth we look for a numerical quadrature such that


with a constant uniformly bounded over . We restrict ourselves with the two-dimensional case, .

As a pre-processing step we may compute a simple polygonal approximation to the curvilinear integration domain . To find the polygonal subdomain , we invoke a root finding procedure (several iterations of the secant method are used in our implementation) to find all intersection points of the zero level set of with the edges of . is defined as a convex hull of these intersection points and the vertices of lying in . The integration of a polynomial function over can be done exactly using Gaussian quadratures on a (macro) triangulation of . The integration over has to be done approximately. We can write

For the numerical integration over and , we consider several approaches. We start with the most straightforward technique, the Monte-Carlo method.

For the Monte-Carlo method we seed points using a uniform random distribution over a narrow rectangular strip containing . Further compute


Further we calculate the variance . If the variance exceeds a predefined threshold , then we increase the set of points used and update the sum on the right-hand side of (16). Noting , a conservative estimate gives .

In the Monte-Carlo method the number of function evaluations per a grid element intersected by is too far from being optimal. A sub-triangulation method below allows to decrease the number of function evaluations per cell and still constitutes a very robust approach.

In the sub-triangulation algorithm we construct a local triangulation of . This is done by finding points on the curvilinear boundary . In our implementation, these points are found as intersections of a uniform ray corn tailored to a basis point on , see Figure 1 for the example of how such local triangulations were constructed for cut triangles and (left plot) and a cut triangle (right plot) of a bulk triangulation (FE functions are integrated over the green area). Further, the integral over a cut element is computed as the sum of integrals over the resulting set of smaller triangles.

Figure 1: Local subdivisions of cut triangles.

The approach can be viewed as building a local piecewise linear approximation of with some , with , i.e. for elements and for elements.

The number of function evaluations per a triangle intersected by is O(), which is better than with the Monte-Carlo method, but still sub-optimal. The first integration method delivering optimal complexity we consider is the Moment-Fitting method from [18].

In the Moment-Fitting method, one first defines a set of points , , for a given cell intersected by . The choice of the points can be done for a reference triangle and is independent on how intersects : can be regularly spaced, come from a conventional quadrature scheme, or even randomly distributed. The boundary of the integration domain consists of straight edges and the curvilinear part , cf. Figure 2. Let , to be the unit outward normal vectors for each and respectively.

Figure 2: Integration domain for the moment–fitting.

The moment–fitting method calculates quadrature weights for the chosen points by finding the least-square solution to the system


for a given set of basis functions , . For example, for P2 FEM we use the basis


The number of points is recommended in [18] to exceed the number of basis functions, i.e., holds.

The integrals on the right-hand side of (17) are evaluated approximately by resorting to divergence-free basis functions. For divergence-free basis functions, one applies the divergence theorem to reduce dimensions of the integrals and relate integration over implicit surface to simple line integrals. In 2D, the div-free basis complementing (18) is given by

Since is a closed piecewise smooth curve, for it holds

Hence, we obtain

Integrals over any interval on the right side can be computed with a higher order Gauss quadrature rule. This allows to build a quadrature for the numerical integration over implicitly given curvilinear edge based on the interior nodes . To this end, one calculates the weights for the set of nodes by solving the following system:

For area integration, take the second set of functions related to as :

The divergence theorem gives

Now the right-hand side values in (17) are (approximately) computed using the above identity and the previously computed surface quadrature rule for the numerical integration over :

This algorithm can be extended to higher order quadratures by expanding the function sets , and ; or to integration over a higher dimensional curvilinear simplex by adding a further moment–fitting step.

The complexity of the moment–fitting integration is optimal, e.g. O(1) of function evaluations per triangle. However, the weights computed by the fitting procedure are not necessarily all non-negative. There is no formal prove of the resulting quadrature accuracy as in (15) or the consistency order. Furthermore, in experiments we observe that the finite element methods with stiffness matrices assembled using moment–fitting can be less stable compared to applying other numerical integration techniques discussed here. Below we consider another integration method of optimal complexity.

Consider the curvilinear remainder and denote by , the intersection points of with . If there are more than 2 such points, then the calculations below should be repeated for each of the simply-connected component of . Choose points on as a nodes of a Gaussian quadrature with weights . Consider the normal vector for the line passing through and . For each point , one finds the point on the boundary such that , , and , where is the -level value for , i.e. for the FEM (9) and for the FEM (10). Secant method finds up to machine precision within a few steps. Further we employ the same 1D quadrature rule to place points {} on each segment (, ). The integral over is computed through


The weights for integrating over may be combined in to write the final quadrature formula . We remark that for problem (10), the factor in (19) is replaced by or depending on the level set of .

Similar to the moment–fitting method, the complexity of the numerical integration is optimal, i.e. O(1) of function evaluations per triangle. All weights are positive, and the accuracy analysis of numerical integration of a smooth function over , which invokes the constructed quadrature for handling boundary terms, is straightforward and outlined below.

We estimate the error of integration of a sufficiently smooth over . The integration uses a conventional quadrature scheme for interior cells of and the quadrature (19) for the curvilinear remainders of cut cells. This composed numerical integral is denoted by . Assume the 2D quadrature over regular triangles has accuracy. By the triangle inequality, we have

We apply the co-area formula to estimate the first term

For the second term we estimate

where is a quadrature we use to integrate over . Decomposing the mesh into interior cells and cut cells (those intersected by ) and applying triangle inequalities, we obtain:

The interior cells are integrated with the error by a conventional method:

For , and polygonal defined earlier, we have

The estimate below assumes that a Gaussian 1D quadrature with nodes is used in the construction of (19). For each interval we introduce the local orthogonal coordinate system , where parameterizes the interval. The graph of the zero level of is the implicit function given by . We note the identity

We assume and to be smooth enough that


with a constant uniform over all and independent of . Applying standard estimates for the Gaussian quadratures, we get:


For a 2D boundary, the number of unfitted regions grows at an order of O(). Then the entire error for all unfitted regions will be .

To satisfy (15), it is sufficient to set . Therefore, for sufficiently smooth the assumption (20) is valid if with the -norm uniformly bounded over all cut triangles and independent of . The latter follows from our assumptions on and . Indeed, let be a suitable polynomial interpolant for on a cut triangle . By triangle inequality we have


Applying the finite element inverse inequality, condition (5), estimate (6) and approximation properties of polynomials we get for the first term on the right-hand side of (22):

with a constant independent of and . Estimating the second and the third terms on the right-hand side of (22) in an obvious way, we obtain


with a constant independent of and . The desired estimate on the -norm of follows from (23), the properties of implicit function and assumptions on .

4 Numerical examples

In this section we demonstrate the results of a few experiments using different numerical integration approaches described in the previous section.

4.1 Integral of a smooth function

We first experiment with computing integral of a smooth function over an implicitly defined domain in . For the domain we choose the annular region defined by the level set function

For given in polar coordinates by

the exact value is known and can be used to test the accuracy of different approaches.

Due to the symmetry, the computational domain is taken to be the square . Further, uniform triangulation with meshes of sizes , , are built to triangulate . To avoid extra geometric error and assess the accuracy of numerical integration, in these experiments we set . All four methods are set up to deliver the local error estimate (15) with or . This should lead to and global accuracy, respectively.

Tables 14 demonstrate that all methods demonstrate convergent results with local parametrization and sub-triangulation being somewhat more accurate in terms of absolute error values. At the same time, only moment–fitting and local parametrization approaches are optimal in terms of the computational complexity. Here we measure complexity in terms of the number of function evaluations. The total number of function evaluations to compute integrals over cut elements and interior elements is shown. Monte-Carlo method appears to be the most computationally expensive. Both Monte-Carlo and sub-triangulation methods become prohibitively expensive for fine meshes so that we make only 4 refining steps with those methods for and only 2 refining steps for . The actual CPU timings (not shown) depend on particular implementation. For a Matlab code we used, the moment–fitting was the fastest among the four tested for a given .

h MF rate MC rate ST rate LP rate
0.1000 7.41e+01 5.50e-01 1.02e-01 1.66e+00
0.0500 6.22e+00 3.57 4.78e-02 3.52 8.53e-03 3.58 7.06e-03 7.88
0.0250 9.14e-02 6.09 1.56e-03 4.94 4.80e-04 4.15 4.46e-03 0.66
0.0125 7.96e-03 3.52 7.93e-04 0.98 1.35e-05 5.15 5.91e-05 6.24
0.0062 1.28e-03 2.64 1.22e-04 2.70 1.17e-06 3.53 1.18e-05 2.32
0.0031 9.05e-06 7.14 2.62e-07 5.49
0.0016 1.11e-05 -0.29 3.99e-08 2.72
0.0008 4.17e-08 8.06 1.91e-09 4.38
0.0004 1.68e-07 -2.01 1.44e-10 3.73
Table 1: The global error and the error reduction rates for the numerical integration using moment–fitting (MF), Monte-Carlo (MC), sub-triangulation (ST), and local parametrization (LP) algorithms to handle cut elements. The table shows results with .
0.1000 1 352 25 802 4 136 423
0.0500 3 718 93 491 16 966 1 608
0.0250 11 726 346 005 68 942 12 750
0.0125 40 144 1 342 357 277 456 42 192
0.0062 146 926 40 645 578 1 113 070 151 022
0.0031 561 912 - 570 104
0.0016 2 194 452 - 2 210 836
0.0008 8 671 104 - 8 703 872
0.0004 34 471 164 - 34 536 700
Table 2: The number of function evaluations for the numerical integration using moment–fitting (MF), Monte-Carlo (MC), sub-triangulation (ST), and local parametrization (LP) algorithms to handle cut elements. The table shows results with .
h MF rate MC rate ST rate LP rate
0.1000 2.80e+01 6.49e-02 5.51e-02 4.16e-02
0.0500 2.79e-01 6.65 8.64e-04 6.23 7.60e-03 2.86 1.13e-04 8.52
0.0250 4.44e-03 5.97 2.25e-05 8.40 7.18e-06 3.98
0.0125 3.61e-04 3.62 9.75e-09 9.52
0.0062 2.11e-05 4.10 1.53e-09 2.67
0.0031 9.59e-07 4.46 7.73e-12 7.63
Table 3: The global error and the error reduction rates for numerical integration using moment–fitting (MF), Monte-Carlo (MC), sub-triangulation (ST), and local parametrization (LP) algorithms to handle cut elements. The table shows results with .
0.1000 1 352 2 341 229 38 696 1 928
0.0500 3 718 117 184 420 308 806 4 870
0.0250 11 726 2 465 102 14 030
0.0125 40 144 44 752
0.0062 146 926 156 142
0.0031 561 912 580 344
Table 4: The number of function evaluations for numerical integration using moment–fitting (MF), Monte-Carlo (MC), sub-triangulation (ST), and local parametrization (LP) algorithms to handle cut elements. The table shows results with .

4.2 Unfitted FEM

In the next series of experiments we solve the Poisson equation with Neumann’s boundary condition in the unit disc domain defined implicitly as with . We are interested in solving (2) with and the right-hand side given in polar coordinates by

The corresponding solution is . In experiments we set and .

Figure 3: Finite element method error for (left) and (right) elements. The error plots are shown for moment–fitting (MF) and local parametrization (LP) used to treat cut elements.
Figure 4: Dependence of the total CPU time to assemble stiffness matrices for (left) and (right) elements. Times are shown for moment–fitting (MF), local parametrization (LP), sub-triangulation (ST) and Monte-Carlo (MC) methods used to treat cut elements.

The bulk domain is triangulated using uniform meshes of sizes , . We experiment with and finite elements. The discrete level set function is the finite element interpolant to the distance function , . Hence, the estimate (6) holds with and , respectively. To be consistent with the geometric error and polynomial order, all four integration methods for cut cells were set up to deliver the local error estimate (15) with and , respectively. According to (11), we should expect and convergence in the energy norm and and convergence in the norm.

Figure 3 shows the error plots for the unfitted finite element method (9) for different mesh sizes. The results are shown with the moment–fitting and local parametrization algorithms used for the integration over cut triangles. The FE errors for Monte-Carlo and sub-triangulation were very similar to those obtained with the local parametrization quadratures and hence they are not shown. These results are in perfect agreement with the error estimate (11). Interesting that using the moment–fitting method for computing the stiffness matrix and the right-hand side leads to larger errors for the computed FE solution.

Further, Figure 4 shows the CPU times needed for the setup phase of the finite element method using different numerical integration tools. For element, the moment–fitting, the local parametrization and the sub-triangulation method show similar scaling with respect to since the complexity is dominated by the matrix assemble over internal triangles, while for elements local parametrization is superior in terms of final CPU times. As expected from the above analysis, the Monte–Carlo algorithm is non-optimal in either case.

4.3 Narrow-band unfitted FEM

In the final series of experiments we apply the narrow-band unfitted finite element method (10) to solve the Laplace-Beltrami equation (3) with on the implicitly defined surface with . The solution and the right-hand side are given in polar coordinates by

Figure 5: and errors for the narrow-band with (left) and (right) bulk elements. The error plots are shown for moment–fitting (MF) and local parametrization (LP) used to treat cut elements.

The bulk domain, triangulations and the discrete level set function are the same as used in the previous series of experiments in section 4.2. For the extended finite element formulation (10) we define the following narrow band domain: