Dual weighted residual error estimation for the finite cell method

3.1 Error representation

The error err_ex from equation (3.3) ought to be computable except for minor perturbations: We assume that there exists a sufficiently precise approximation J^(n+k) of J for a fixed k ∈ ℕ, such that the resulting perturbation error is negligibly small. Moreover, the representation should allow for a separation of two error sources, i.e., the discretization and the quadrature. We derive such a representation by adapting the proof stated for the standard FEM case in [33, Prop. 3.1], where a representation of the error with respect to any perturbation v_h ∈ V_h of the discrete solution is provided. In the FCM setting, we show that the error is composed of the sum of a discretization-related error term, a quadrature-related error term, and some terms which are assumed to be negligibly small (e.g., of higher order). The representation requires computable replacements u⁺, z⁺ for the unknown solutions u, z, as well as improvements G^(n+k) of the approximations G⁽ⁿ⁾ for each G ∈ {A, A′, Â, Â′, J, J′} with the property that |G^(n+k)(⋅) − G(⋅)| is negligibly small. Recall that we assume that a quadrature scheme allowing for these improvements is available. For instance, when the space-tree is used, the approximations with index n + k may be obtained from approximations with index n by refining the space-tree k times.

Before stating the representation, we briefly explain the occurring terms. The residuals ρ, ρ^* are defined as

Here, we abbreviate by v₀ the extension of v|_Ω by zero onto Ω̂. The residuals ρ, ρ^* occur in the definitions of the terms eD(n+k) (u⁺|_Ω, z⁺|_Ω), e_NS,D (u⁺|_Ω, z⁺|_Ω) related to the discretization error. These terms are defined for v, w ∈ V by

We assume the approximations u⁺ and z⁺ to be such that the residuals in the errors u₀ − u0+ and z₀ − z0+ are of higher order compared to those in u0−(uh(n))0 and z0−(zh(n))0.

Abbreviating the errors e:=u−uh(n)|Ω  and  e∗:=z−zh(n)|Ω, we obtain the quadrature-related error term eQ(n+k) and additional contributions e_NS,Q, e_NS,J, and e_NS,𝓛:

These additional terms are assumed to be negligibly small for k sufficiently large due to the fact that they consist of evaluations of differences G − G^(n+k). To see that eQ(n+k) can be regarded as a quadrature error term, assume that exact integration is available, i.e., Â^(n+k) = Â. Then, we may abbreviate u_h := uh(n), z_h := zh(n), and

Therefore, the term eQ(n+k) vanishes if all operators are exact. If only approximations to the operators are available or a perturbed solution is inserted, the term eQ(n+k) will be nonzero in general and, thus, may be regarded as a perturbation error. This error may be caused by, e.g., numerical quadrature or by an iterative method as in [33].

Finally, the term e_ε describes the error incurred by the FCM approximation in Ω̂ ∖ Ω, which is assumed to be negligibly small for ε sufficiently small. It is defined as

Typically, the value of e_ε is 𝓞(ε), see also (2.3).

We assume that several error terms in (3.4), (3.5), and (3.6) are negligibly small. In Section 5, we will see that the adaptive algorithm presented in the next section can be configured in such a way that these terms are indeed negligibly small, which justifies our assumption. Ignoring these terms, we arrive at the following approximate error representation.

The effectivity index (or overestimation index) is defined as

3.2 Approximation of the solutions of the continuous problems

The unknown quantities u and z are approximated by computable functions u⁺ and z⁺. To this end, several methods have been proposed in the literature. The first method computes approximations by solving the discrete problems in a finite-element space of higher polynomial degree, e.g., by doubling each local polynomial degree [6]. However, this is too expensive except for simple test problems. Usually, patched meshes are employed, i.e., for an element of 𝓣_h with its parent element P in the mesh refinement history, all 2^d children of P are elements of 𝓣_h. This implies that whenever an element is refined, all its siblings are refined as well. In this case, the patches of the mesh can be joined to form a finite-element space V_2h,2p of double mesh width and double polynomial degree. The computation of u⁺, z⁺ is then approximately as expensive as the computation of uh(n),zh(n).

Another method requiring patched meshes uses local higher-order interpolation to compute more accurate approximations by, again, viewing each patch as a single element of a coarser mesh and doubling the polynomial degree [2]. This eliminates the need of computing additional discrete solutions. Under certain regularity assumptions, it can be shown that the error incurred is of higher order (see [2, Section 5.2]). However, one has to take care that the resulting functions are elements of V by ensuring continuity of the interpolation on the boundary of the patches, which is difficult when hanging nodes are present [30].

3.3 Localization

To perform the finite-cell mesh adaptation, the discretization error eD(n+k) has to be localized to nonnegative elementwise contributions ηD,T(n+k). To this end, several methods are available, for all of which good effectivity indices have been demonstrated for many practical problems (see, e.g., the references in the introduction). A first method applies elementwise partial integration leading to an inner residual and a boundary residual involving integrals over the boundary of each element. An obvious disadvantage of this method is the possibly costly computation of strong residuals and jump terms. Also, the strong adjoint residual formulation may not even be available [35]. In the context of the FCM, another disadvantage consists in the necessity of determining intersections between the boundary of the domain and the element boundaries, as these intersections are not required for the application of the FCM. Also, these intersections have to be determined with great precision in order not to introduce additional errors. A second method, which uses the variational formulation directly, has been proposed by Richter and Wick [35]. Here, a partition of unity for the nodes of the finite-element mesh is inserted into the error representation.

A third method known as the algebraic filtering approach has been described by Braack and Ern [9] which is also based on the variational formulation. The method relies on patched meshes and the associated canonical finite-element spaces V_2h,p and V_2h,2p formed by treating each patch as an element of degree p and 2p, respectively. In order to reconstruct higher-order solutions, interpolation and filtering operators are defined, which we briefly describe in the case of meshes without hanging nodes similar to [35]: As the spaces V_h,p and V_2h,2p have the same numbers of unknowns in the same Lagrange points, one can define an interpolation operator i^* : V_h,p → V_2h,2p by assigning to v_h ∈ V_h,p the element i^*v_h of V_2h,2p uniquely determined by the values of v_h in those Lagrange points. The filtering operator is defined by π_2h := id − i_2h for a finite-element interpolation operator i_2h : V → V_2h,p. The name of the method stems from the observation that π_2hv_h is a strictly local algebraic process acting on the coefficient vector v ∈ ℝ^N of v_h ∈ V_h,p, since π_2hv_h = v_h − i_2hv_h = ∑j=1Nv_j (φ_h,j − i_2hφ_h,j) = : ∑_j (π_2hv)_jφ_h,j for the basis (φh,j)j=1N of V_h,p (see [35, eq. (40)]). Finally, the localization of eD(n+k) to nodal contributions is as follows:

The node-wise error contributions ηD,j(n+k):=|η~D,j(n+k)| are then used for the marking step in the adaptive procedure, e.g., by refining all elements touching node j or by explicitly assigning an elementwise indicator ηD,T(n+k) based on the node-wise contributions and performing the usual elementwise refinement.

To measure the quality of the localization, we define the indicator index

similar to [35, eq. (26)] which takes into account the overestimation of the discretization error caused by taking the absolute value of possibly negative local values η~D,T(n+k).

The quadrature error estimator value eQ(n+k) is used to check if the quadrature mesh resolves the computational domain Ω sufficiently well in order to obtain a meaningful discrete solution. For its localization, we split it canonically into element-wise contributions

and utilize the same localization techniques as for the discretization error.

5.1 Quarter disk

The first domain is the quarter disk Ω := B₁(0) ∩ (0, 1)² with Dirichlet boundary part Γ_D := ([0, 1] × {0}) ∪ ({0} × [0, 1]) and Neumann boundary part Γ_N := ∂ Ω ∖ Γ_D. We face the difficulty that the circular domain cannot be represented exactly by quadrilateral finite elements. Thus, we embed Ω into the rectangle Ω̂ := (0, 1)². We solve the non-linear diffusion–reaction equation − Δu + u³ = f, u|_ΓD = 0, ∂_n u|_ΓN = g_N. Its weak form reads

with Gateaux derivative

The functions f and g_N are chosen such that the solution u(x, y) := sin(πx) sin(πy) is obtained. The quantity of interest J is selected according to the analytic dual solution

given in polar coordinates.

First, we investigate the dependence of the a posteriori error estimator on the number of additional quadrature refinements denoted by k. We fix ε = 10⁻¹⁰ and the number of mesh elements in 𝓣 to 262 144 and use a constant quadrature refinement of α(T) := n for each T ∈ 𝓣. In Fig. 3, the different contributions to the error identity (3.7) are plotted. We use a quadrature with a piecewise quadratic boundary approximation to evaluate the integrals almost exactly, if needed. We consider the piecewise constant and piecewise linear boundary approximation in the numerical quadrature described in Section 2 and set n := 4 for the piecewise constant and n := 1 for the piecewise linear scheme. For both quadrature techniques, we generally obtain the same results. The terms comp, e_D, e_Q, and e_NS,D are approximately constant as expected, if k is sufficiently large. For a piecewise constant boundary approximation, we need k ⩾ 4 and in the linear case k ⩾ 1. The terms e_NS,Q, e_NS,J, and e_{NS, 𝓛}, which we want to assume as negligibly small, are decreasing with increasing k. We denote the mesh width of the quadrature mesh by h_Q and note that an incrementation of k by 1 results in halving h_Q. For the piecewise constant boundary approximation the terms e_NS,Q, e_NS,J, and e_{NS, 𝓛} are of order 𝓞(h_Q) and for the piecewise linear one we find O(hQ2).

Fig. 3

Single contributions to the error identity for increasing number of additional quadrature refinements.

Second, we examine the quadrature level n. To this end, we choose ε = 10⁻¹⁰ and the number of mesh elements in 𝓣 again as 262 144. Furthermore, we set k = 4 for the piecewise constant boundary approximation and k = 2 for the linear one. The results are presented in Fig. 4. The terms comp and e_D are approximately constant from the beginning in the case of linear boundary approximation. The term e_NS,D becomes constant for n ⩾ 4 as expected. Also, all other terms are of order O(hQ2). In the case of a piecewise constant boundary approximation in the quadrature, we find that the terms comp and e_D become approximately constant for n ⩾ 10, while we do not see this for e_NS,D here. The reason is that the values of e_NS,D are polluted by the quadrature error even for n = 12. All other terms are of order 𝓞(h^Q).

Fig. 4

Single contributions to the error identity for increasing number of quadrature refinements.

Third, we consider varying mesh sizes. The results are depicted in Fig. 5 where ε = 10⁻¹⁰ and a piecewise linear boundary approximation is used. We observe the optimal convergence in comp of order 𝓞(h²). The boundary approximation is accurate enough such that no further refinements of the quadrature are needed. We see higher order convergence for e_NS,L, while the terms e_Q, e_NS,J, and e_NS,Q are of order 𝓞(h²) but are essentially smaller than the error in J. The term in which we are mainly interested here is e_NS,D measuring the error with respect to u⁺ and z⁺. This term is also of order 𝓞(h²) and is considerably smaller.

Fig. 5

Single contributions to the error identity for varying mesh size using a linear boundary approximation with n = k = 2.

Fourth, we take a look at the dependence of the a posteriori error estimator on the regularization parameter ε. We fix the number of mesh elements in 𝓣 to 262 144 and approximate the boundary in the quadrature by piecewise linear functions with n = k = 2. The results are illustrated in Fig. 6. We observe that for ε < 10⁻⁴ the error is approximately constant as well as all other terms except for e_Q, e_NS,D and e_ε. The numerical error estimator becomes constant for ε < 10⁻⁶, which illustrates its stronger dependence on ε. This is also observed for e_NS,D, which is constant for ε < 10⁻⁷. The term e_ε is of order 𝓞(ε) as expected. For large ε, we find the same behaviour also for the error and error estimator. Consequently, choosing ε = 10⁻¹⁰ is a sufficiently small value to ensure that the error due to the regularization is negligibly small.

Fig. 6

Single contributions to the error identity for decreasing ε

Finally, we test the adaptive algorithm in this example, where we set the minimum quadrature refinement l = 1 and start with α = 0 everywhere. We work with k = 4 for the piecewise constant boundary approximation and with k = 2 for the piecewise linear one. We expect uniform mesh refinements in the original domain because of the smoothness of the primal and dual solutions. Outside of the original domain, additional refinements should not occur. This expected structure of the adaptive meshes is produced by the adaptive algorithm, cf. Fig. 7. We find the optimal convergence order for adaptive and uniform refinement as depicted in Fig. 8, where N denotes the total number of degrees of freedom. The detailed results of the adaptive algorithms for piecewise constant and linear boundary approximations are listed for ρ = 1 in Table 1 and 3 respectively. For the piecewise linear boundary approximation, we find effectivity indices, which are very close to 1. It should be remarked, that no additional quadrature refinements are carried out. The behaviour in the case of the piecewise constant boundary approximation is completely different. Here, we need additional refinements of the quadrature mesh on each refinement level. The reason is the reduced convergence order of the quadrature rule. Furthermore, the numerical error is significantly higher than in the case of the piecewise linear boundary approximation and pollutes the a posteriori error estimate leading to effectivity indices which are oscillating and which are not close to one. The higher-order reconstruction used to evaluate the error identity numerically relies on asymptotic properties of the finite element solution in the nodes, cf. [2, Section 5.2.ii]. However, the relatively large numerical error disturbs these properties. In Table 2, we list the effectivity indices for a fixed mesh under quadrature refinement. We observe that the estimated numerical error should be a factor about 10 smaller than the estimated discretization error to obtain good effectivity indices. However, although the effectivity indices are not good, the adaptive algorithm performs optimally as shown in Fig. 8. All in all, we save a large amount of computational work in the numerical quadrature compared to choices like ρ = 0.1 or ρ = 0.01. However, there is a loss of accuracy in the error estimate.

Fig. 7

Adaptive meshes for the quarter disk example with ρ = 1.

Fig. 8

Convergence results for the quarter disk example with ρ = 1.

For the piecewise constant boundary approximation, we plot the distribution of the quadrature level in Fig. 9. The quadrature level is low at the ends of the fictitious boundary and high in the interior. Also, it can be seen that the quadrature level is higher on elements with a larger diameter. This is due to the fact that a certain quadrature level is required to ensure coercivity and the refinement of the finite-cell mesh induces a refinement of the quadrature.

Fig. 9

Quarter disk: Distribution of the quadrature level in the 14^th iteration of adaptive algorithm for the piecewise constant boundary approximation.

5.2 Circular domain with reentrant corner

In a second series of experiments, we consider the circular domain with reentrant corner Ω := B¹(0) ∖ ([0, 1] × [−1, 0]). The Dirichlet boundary part is Γ_D := ([0, 1] × {0}) ∪ ({0} × [−1, 0]) and the Neumann boundary part is Γ_N := ∂ Ω ∖ Γ_D. For the discretization via the finite-cell method, we embed Ω into the L-shaped domain Ω̂ := (−1, 1)² ∖ ([0, 1] × [−1, 0]). The initial finite-cell mesh 𝓣₀ consists of 3 ⋅ 16 square elements of degree 1. Again, we consider the diffusion–reaction equation −Δu + u³ = f, u|_ΓD = 0, ∂_n u|_ΓN = g_N. The functions f and g_N are chosen such that u(r, φ) = r^2/3sin23φ is the solution given in polar coordinates. Except for the circular arc on the Neumann boundary, this problem closely resembles the classical L-shaped domain problem. In particular, it features a corner singularity in (0, 0) leading to a reduced regularity. Consequently, we cannot expect that the finite cell method based on uniform refinement yields optimal algebraic convergence rates.

We aim to control the error of |∇ u|² at the point p=(cos54π,sin54π)∈∂Ω, which we approximate by J(v) := ∫_Ω̂Ĩ_S|∇ v|², where Ĩ_S with S = B_0.01(p) ∩ Ω is a smoothed indicator function given on S by

where w(r) = −10⁶ ⋅ r³ + 3 ⋅ 10⁴ ⋅ r² − 3 ⋅ 10² ⋅ r + 1.0. Since the exact solution u is known, the functional J can be evaluated up to machine precision on S. Here, we compute J(u) ≈ 0.0017523852871125355.

The adaptive meshes generated in the adaptive algorithm with ρ = 0.01 are depicted in Fig. 10. We find strong refinements close to the point p and in the reentrant corner as well as between these two regions, which is expected. The convergence is examined in Fig. 11. After some quadrature refinements in the beginning, we find optimal convergence of order 𝓞(N⁻¹) for the adaptive approach. Using uniform refinements, the convergence order is reduced as expected. The results for piecewise constant and linear boundary approximations in the quadrature are almost identical. However, the computation with the piecewise constant boundary approximation stops earlier due to the huge number of quadrature points. The two approaches deliver similar results which are given in detail in the Tables 4 and 5. However, we need approximately 690 times more quadrature points when we use the piecewise constant boundary approximation to achieve the same accuracy as with the linear approach. This is also substantiated by the significantly higher quadrature levels for the piecewise constant boundary approximation. In Fig. 12, the distribution of the quadrature level is shown for the piecewise linear boundary approximation. We find a high level near to S if the cells are comparatively coarse. Furthermore, in the upper left corner, we find high quadrature levels. The effectivity indices lie in an acceptable range between 0.8 and 1.3 for higher N even though the problem is of low regularity, cf. Tables 4 and 4.

Fig. 10

Adaptive meshes for the L-domain example.

Fig. 11

Convergence results for the L-domain example.

Fig. 12

L-domain: Distribution of the quadrature level in the 13^th iteration of adaptive algorithm for the piecewise linear boundary approximation.

5.3 B-domain

In this section, we consider a domain, which we call B-domain. It is illustrated together with the initial mesh in Fig. 13a, where the used notation is also defined. The domain has two holes and a strong singularity in the point (0, 0), which is located on the fictitious boundary and not in a node of the mesh. The underlying problem formulation is model plasticity with linear isotropic hardening, we refer to [38] for a complete description. It is given by the PDE

Fig. 13

Sketch of the B-domain example, adaptive mesh, and distribution of the quadrature level in the 8^th iteration of adaptive algorithm.

Here, σ₀ > 0 denotes the yield stress and 0 ⩽ζ ≪ 1 the hardening parameter. We choose σ₀ = 1 and ζ = 0.01. We assume homogeneous Neumann boundary conditions on the boundary of the holes and on two small parts of the outer boundary, Γ_N = Γ_I ∪ Γ_U, and homogeneous Dirichlet boundary conditions on the outer boundary, Γ_D = Γ_O ∪ Γ_F. Since the homogeneous Dirichlet boundary conditions on Γ_F cannot be realized in the usual finite element way in the finite cell method, we approximate them by a penalty approach with a large penalty parameter γ ≫ 0. The initial γ is set to 10⁻¹⁰ as well as ε = 10⁻¹⁰. The right hand side f is given by a constant function with f≡ 2. The weak form of equation (5.1) reads

It should be remarked that a need not be three times directional differentiable. Nevertheless, our approach works as shown in the following results. We are interested in the quantity of interest J given by

where ω_D(x, y) = 10⁸ ⋅ (x + 0.3)²(x + 0.1)²(y + 0.1)²(y − 0.1)² and ω_F(φ, r) = φ²(φ − 0.5 π)².

We have seen in the last example that the results for a piecewise linear and a piecewise constant boundary approximation are almost the same except for the number of quadrature points needed. To shorten the presentation we focus here on the linear boundary approximation with ρ = 0.01 and k = 2. In Fig. 13b, the mesh in the 14^th iteration of the adaptive algorithm is depicted. We find strong refinements around the singularity in (0, 0) as well as around S and Γ_J and around the edges of the holes, which are expected. Since the exact solution of this example is unknown, we cannot compute J(u) exactly. Hence, we approximate J(u) using extrapolation techniques based on the results of the adaptive algorithm and obtain

This numerically computed reference value is used in the calculation of the error and the effectivity index. We compare the convergence properties of the uniform and the adaptive approach in Fig. 14. The adaptive approach outperforms the uniform approach by far. The results are detailed in Table 6. For finer meshes the effectivity indices are in the range of 0.8. This is a very good value for such a complex example especially in the presence of the nonsmooth differential operator. The indicator indices are smaller than 1.85 in most iterations, which is also a good result. We find that the overall quadrature refinement levels are lower than for the L-domain example. The distribution of the quadrature level in the 8^th iteration of the adaptive algorithm for the piecewise linear boundary approximation is depicted in Fig. 13c. The quadrature levels are higher on coarser mesh cells and smaller on fine mesh cells similar to the other examples. It should be pointed out that the quadrature is exact on the straight lines of the holes such that no quadrature refinements at all are carried out. The level is one because the quadrature level is adjusted to the minimal quadrature level during the refinement process of the mesh cells.

Fig. 14

Convergence results for the B-domain example.

5.4 Implementation details and analysis of computing time

Finally, we give some remarks on the implementation and take a look at the computing time of the adaptive algorithm. We have implemented the numerical quadrature in parallel on shared memory computers, so that all integration operations as the assembling of vectors and matrices or the evaluation of J and the error estimator are performed in parallel with a good strong scaling. However, all other parts of the code especially the linear solvers are not parallelized. The nonlinear systems are solved with a damped Newton iteration. In this example, we need up to 40 Newton iterations with strong damping in the first iterations, when the mesh is refined. We refer to [38] for more information concerning the solution algorithm used for this class of problems. If a quadrature refinement is performed, the quadratic convergence of the Newton scheme is realized and only two or three Newton iterations are needed. The linear system in each Newton step is solved with a conjugate gradients (CG) method preconditioned by a symmetric SOR scheme.

In Fig. 15, we compare the computing time of the different parts of the adaptive algorithm. We summarize all operations concerning memory allocation, quadrature and parallelization preparation and so on under the designation ‘setup’. By ‘assemble’ all assembling operations during the solution of the primal problem are collected, while ‘solve’ is the cumulated computing time for solving all linear systems in this part. The counterparts concerning the dual problem are named by ‘dual assemble’ and ‘dual solve’. The evaluation of the quantity of interest is given by ‘evaluate J′. The process of the calculation of the error estimator η is summarized in ‘evaluate η′, where we note that several operations especially concerning the localization are not parallelised. The adaptive refinement routines for the mesh and the quadrature are considered under the designation ‘adaptivity’. The graphical output of the mesh and the output of the data for the tables and graphs are collected in ‘post processing’. The computations are carried out on a Sun Fire compute server with 8 AMD Quad-Core Opteron 8356 CPUs (2.3 GHz) and 64GB RAM, where we use 16 cores for the calculations concerning this example. Figure 15 shows that the computing time for a large number of unknowns is dominated by the solution process of the linear systems. In the last step nearly 92% of the computing time is needed for the solution of the linear systems. Approximately 1.2% is used for the assembling and 5.9 for the solution of the dual problem. All other parts are less than 1%. Note that fewer Newton iterations are required during the quadrature refinement steps. Thus the computational amount of solving the dual problem increases and the computational amount of solving the primal problem decreases. It is an interesting question for further research how the distribution of the computing time changes if the solver is also parallelized and if more sophisticated preconditioners, for instance, based on algebraic multigrid adapted to the special challenges of the finite cell method are used. However, such research is out of the scope of the article at hand.

Fig. 15

Percentage of the computing time in each adaptive iteration.