# Dual weighted residual error estimation for the finite cell method

• Paolo Di Stolfo , Andreas Rademacher and Andreas Schröder

## Abstract

The paper presents a goal-oriented error control based on the dual weighted residual method (DWR) for the finite cell method (FCM), which is characterized by an enclosing domain covering the domain of the problem. The error identity derived by the DWR method allows for a combined treatment of the discretization and quadrature error introduced by the FCM. We present an adaptive strategy with the aim to balance these two error contributions. Its performance is demonstrated for several two-dimensional examples.

Classification: 65N15; 65N85

## 1 Introduction

The finite cell method (FCM) is a well-established variant of the general fictitious domain approach [17, 18, 36] and was developed by Parvizian, Düster, and Rank [15, 27]. It has been applied to a vast number of both linear and nonlinear problems, including linear elasticity in 2D and 3D [15], shell problems [32], biomechanical problems [41, 42], wave propagation [21], elastoplasticity [1], and topology optimization in structural mechanics [28].

The basic idea of the FCM is to replace the possibly complicated domain of the problem by an enclosing domain of a geometrically simple shape, e.g., a (paraxial) quadrilateral in 2D or (paraxial) hexahedron in 3D. As the enclosing domain can be trivially subdivided into (paraxial) quadrilateral or hexahedral cells, mesh generation is simplified substantially. The finite element space is constructed on these cells, from which the name of the method is derived. To recover the geometry of the original problem, the integrals in the variational formulation of the problem are approximated by quadratures defined on the covering mesh of finite cells. To this end, an approximation of the original domain of sufficient quality has to be available, which is typically provided by a separate quadrature mesh. However, this approximation introduces a quadrature error which is assumed to be lower than the discretization error. A first mathematically rigorous investigation of the FCM for exact integration and certain boundary conditions as well as numerical experiments for inexact integration are provided in [12].

While it has become standard for modern finite-element techniques to include a posteriori error control and adaptivity, error estimators have neither been derived nor applied to the FCM to this day. In this work, we focus on the dual-weighted residual error (DWR) estimation method, which has become one of the most popular a posteriori techniques for standard finite elements in the last two decades. It is based on the preliminary work by Eriksson et al. [16] and was developed by Becker and Rannacher [5]. The DWR method allows for goal-oriented error estimation and, thus, supports more general, user-defined, possibly nonlinear expressions to be estimated, such as norms, point values, averages, or lift and drag coefficients, see [2] for an overview. The method relies on representing the error in terms of the solution of a dual problem, which is typical as duality arguments are the basis of many techniques in so-called goal-oriented error control [26, 29]. The DWR method has been applied to many practical problems including fluid mechanics, chemically reactive flows, and fluid--structure interaction (see, e.g., [3, 8, 19, 34, 39, 40]), as well as simplified Signorini and (frictional) contact problems (see [7, 31, 37]).

A posteriori error estimates are well-developed with respect to exact discrete solutions, i.e., solutions determined with no computational error incurred by, e.g., iterative methods or inexact integration. However, there are only a few publications dealing with a posteriori error estimates for inexact discrete solutions determined by an iterative process, such as the multigrid method [4] or Newton’s method in the context of nonlinear problems [33]. A common idea of these approaches is to apply a stopping criterion that is based on balancing the discretization error with the iteration error.

In this article, we discuss the derivation and application of the DWR method in the FCM context in order to estimate both the discretization error and the quadrature error with respect to a goal functional, along with an adaptive strategy with the aim to balance these two contributions by either refining the finite-cell mesh or its associated quadrature mesh. The balancing is to be understood in the sense that, instead of striving for a very small quadrature error as it is usually the case for standard finite elements, the quadrature mesh is refined so that the quadrature error is only small enough to provide a useful discrete solution. We utilize the localization strategy for the DWR method by Braack and Ern [9] which does not require jumps over element facets and, thus, is well-suited for the FCM.

The outline of the paper is as follows. In Section 2, we give an overview of the FCM. Section 3 discusses the DWR method for goal-oriented a posteriori error estimation and provides the error identity containing terms representing the discretization and the quadrature error. An adaptive strategy realizing both discretization and quadrature mesh refinements is discussed in Section 4. In Section 5, numerical experiments for several 2D examples with different characteristics are presented. Finally, conclusions are drawn in Section 6.

## 2 Abstract framework for the finite cell method

In this section, we present a general nonlinear setting for the finite-cell method (FCM). For this purpose, let Ω ⊂ ℝd be a bounded domain and Ω̂Ω be a paraxial d-dimensional interval (i.e., a rectangle for d = 2 or a cuboid for d = 3), and let ΓDΩ be the Dirichlet boundary part.

Given a Hilbert space V of functions defined on Ω with its dual space V* and an operator A : VV*, we aim to find a solution uV such that

A(u)(φ)=0φV(2.1)

where we assume that (2.1) is uniquely solvable. Furthermore, we assume there exists a space of functions defined on Ω̂ extending V, i.e., |ΩV. Also, we assume there exists an operator Â : * such that

A(v)(φ)=A^(v0)(φ0)v,φV

where w0 denotes the extension by zero onto Ω̂ of a function w defined on Ω.

In the discrete setting, a triangulation 𝓣h of Ω̂ into intervals and a finite-element space Vh on 𝓣h can be constructed easily due to the simple form of Ω̂. In the framework used in the following, we assume Vh|ΩV, i.e., the space of restrictions is conforming. The discrete problem is to find a solution uhVh such that

A^(uh)(φh)=0φhVh.(2.2)

It is assumed that the contributions |Â(uh)(φh) − A(uh|Ω)(φh|Ω)| are sufficiently small so that the model error can be neglected. To illustrate the spaces and operators, we consider a 2D example based on the Poisson model problem

A(u)(φ):=ΩuφΩfφ=0φV

where Ω : = (0, 1)2B1(0) is the quarter disk, fL2(Ω), the space V := HΓD1(Ω) = {vH1(Ω); v|ΓD = 0}, and V* = HΓD1(Ω)* for the Dirichlet part ΓD := [0, 1] × {0} ⊆ ∂ Ω. As the operator Â, we may choose

A^(v)(φ):=Ωvφ+εΩ^ΩvφΩfφv,φV^

where Ω̂ := (0, 1)2, ε ≈ 0 (e.g., ε = 10−12) is a positive parameter large enough to secure coercivity, and

V^:={v:vL2(Ω^),v|ΩV,v|Ω^ΩH1(Ω^Ω)}.

For the discrete setting, we introduce a triangulation 𝓣h of Ω̂ consisting of four square elements and define Vh to be the H1(Ω̂)-conforming finite-element space of degree 1 on 𝓣h respecting the Dirichlet boundary condition on ΓD, implying Vh|ΩV and Vh. Since, here, the boundary of Ω matches a union of facets, Dirichlet boundary conditions can be applied in a strong manner. However, in general, the Dirichlet boundary is non-matching with Ω̂, i.e., it does not equal the union of some facets in 𝓣h. In this case, Dirichlet boundary conditions may be applied weakly by, e.g., Nitsche’s method [25, 43]. The contributions

|A^(vh)(φh)A(vh|Ω)(φh|Ω)|=ε|Ω^Ωvhφh|(2.3)

result in a model error of O(ε) in the energy norm (see [12, p. 1047]).

While the operator Â is defined on Ω̂, the domains of the involved integrals may depend on Ω. Therefore, for the computation of Â in the discrete setting, numerical integration has to be performed. In the context of the FCM, this usually involves an approximation of Ω by geometrically simple objects. These approximations result in approximate operators Â(∼) and perturbed discrete problems

A^()(uh())(φh)=0φhVh(2.4)

yielding perturbed discrete solutions uh()Vh. A geometrically simple replacement for Ω used in practice is the spacetree with its specializations to two and three dimensions commonly referred to as quadtree and octree, respectively [13]. Here, to each element T ∈ 𝓣h, a set of intervals QT is assigned via a number α(T) ∈ ℕ0 indicating the number of recursive refinements of T towards the boundary ∂ Ω. The set QT is generated by the following recursive procedure:

1. Set i := 0, QT(i) := {T}.

2. If i = α(T), then QT := QT(i), exit. Otherwise, replace each interval in QT(i) that is intersected nontrivially by ∂ Ω by 2d sub-intervals, yielding QT(i+1). Increase i by 1 and go to step 2.

Finally, as an approximation Ω(∼) of the domain of integration Ω, one may use the union of all intervals in any QT having non-trivial intersection with Ω. The union of the remaining intervals is then an approximation to Ω̂Ω. Similarly, an approximation to Ω̂ can be obtained. For the approximation of the integrals involved in the operators Â(∼), the usual quadrature rules used in the finite-element context are applied on each interval.

The result of the procedure is visualized in Fig. 1 for a finite-cell mesh for the quarter disk, where the unit square is subdivided into four equally sized elements, along with the assigned number of recursive refinements.

Fig. 1

The presented procedure for establishing the space-tree only makes use of a point-in-domain test, which is typically applied to sample points of each interval (e.g., the four vertices of a rectangle). Thus, the condition whether an interval intersects ∂ Ω nontrivially is checked only approximately. Due to its simplicity, the space-tree can be easily applied to complicated geometries, e.g., generated in the context of constructive solid geometry. An obvious disadvantage is the fact that it offers only a piecewise constant approximation to ∂ Ω. Therefore, a high number of recursive refinements may be required to approximate the domain sufficiently well. For domains with smooth boundaries, higher-order approximations of the boundary may be used. In the context of the FCM, several improvements over the space-tree have been developed (see, e.g., [11, 20, 22, 23, 24]). However, implementing these improvements is usually rather involved. Whenever the boundary is of a particular simple shape (such as a circle), one may apply the following improvements: Firstly, the quadtree refinement on a given cell can be performed so that the intersection points of the boundary with the cell coincide with vertices of the quadtree, whenever possible (see Fig. 2a). Secondly, one may derive a piecewise linear approximation of the boundary by subdividing the elements by horizontal or vertical line segments, computing intersection points of the boundary with the segments, and replacing the boundary by line segments, so that a subdivision into triangles and quadrilaterals is produced (see Fig. 2b).

Fig. 2

Visualization of the described improvements applicable in the case of the quarter disk

In the following, we do not confine our treatment to the space-tree, but permit a general quadrature scheme which generates an approximation Ω(∼) that allows for refinements, i. e., given an initial approximation Ω(0) := Ω(∼) for Ω, a sequence of domains Ω(n), n ⩾ 1, may be generated fulfilling Ω(n+1)ΩΩ(n)Ω, so that each approximation is an improvement of the previous one. The operators with integrals defined on these approximate domains are indicated with superscript (n), such as A(n). Also, we make the natural assumption that the refinement of the quadrature occurs on the level of the finite cells, i.e., there is a function αn : 𝓣h → ℕ0 indicating the number of refinements of the quadrature scheme used to approximate TΩ in Ω(n). An improvement Ω(n+1) can be generated by assigning αn+1(T) := αn(T) + 1 for each T ∈ 𝓣h. Alternatively, one may only improve a subset 𝓢 ⊂ 𝓣h by

αn+1(T):=αn(T)+1,TS,αn+1|ThS:=αn|ThS.

In particular, note that the theoretical considerations presented in Sect. 3 are valid for any domain approximation scheme that generates approximate operators, each of which is an improvement of the previous one.

## 3 The dual weighted residual method

In this section, we tailor the dual weighted residual (DWR) method for a (possibly nonlinear) problem and a (possibly nonlinear) goal functional to the FCM setting. To this end, we consider the spaces and operators introduced in Section 2 and, in addition, assume A to be three times Gateaux-differentiable. Moreover, let J : V → ℝ be a three times Gateaux-differentiable goal functional. For the kth order Gateaux derivative of a function g: XY for Banach spaces X, Y in a point xX, we adopt the usual identification g(k)(x)(ψ1, …, ψk) = g(k)(x)(ψ1)⋯(ψk), indicating that g(k)(x) is k-linear.

For the unique solution uV of the problem in (2.1), we may formulate the following trivial optimization problem which connects the problem with the goal functional:

u=argminφVJ(φ) subject to A(u)(φ)=0φV.

Introducing the Lagrangian 𝓛 : V × V → ℝ with 𝓛(v, w) := J(v) − A(v)(w), we seek a Lagrangian multiplier zV such that (u, z) is a stationary point of 𝓛, yielding

L(u,z)(φ,ψ)=(J(u)(φ)A(u)(φ,z),A(u)(ψ))=0φ,ψV.

Thus, in addition to seeking a solution u in (2.1), we seek a function z that solves the dual problem

A(u)(φ,z)=J(u)(φ)φV.(3.1)

In the FCM setting, the nonconformity VhV and the approximation of operators and functionals have to be taken into account in the computation of the discrete solutions. Let Ω(n) be an approximation for Ω that allows for refinement. Similar to the case of A and Â(n) described in Section 2, we assume that approximations J(n), J′(n) for J, J′ exist, where all integrals on Ω occurring in the definition are replaced by a quadrature rule on Ω(n). Instead of the discrete dual problem, its perturbation

A^(n)(uh(n))(φh,zh(n))=J(n)(uh(n))|Ω(φh|Ω)φhVh(3.2)

is solved. We are interested in a representation of the exact error

errex:=J(u)J(uh(n)|Ω).(3.3)

In the next section, we derive an error representation for errex. It includes non-computable terms, which we assume to be negligibly small, and higher-order approximations to continuous solutions. The question of how to compute these approximations is addressed in Section 3.2. Also, localization techniques for the discretization and the quadrature error are discussed in Section 3.3.

### 3.1 Error representation

The error errex 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 vhVh 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(n) 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

ρ():=A^(n+k)uh(n)0()ρ():=J(n+k)uh(n))|Ω(|Ω)A^(n+k)uh(n)0,zh(n)0.

Here, we abbreviate by v0 the extension of v|Ω by zero onto Ω̂. The residuals ρ, ρ* occur in the definitions of the terms eD(n+k) (u+|Ω, z+|Ω), eNS,D (u+|Ω, z+|Ω) related to the discretization error. These terms are defined for v, wV by

eD(n+k)(v,w):=12(ρ(w0(zh(n))0)+ρ(v0(uh(n))0))eNS,D(v,w):=12(ρ(z0w0)+ρ(u0v0)).(3.4)

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

Abbreviating the errors e:=uuh(n)|Ωande:=zzh(n)|Ω, we obtain the quadrature-related error term eQ(n+k) and additional contributions eNS,Q, eNS,J, and eNS,𝓛:

eQ(n+k):=A^(n+k)uh(n)zh(n)eNS,Q:=A^(n+k)uh(n)0zh(n)0Auh(n))|Ωzh(n)|ΩeNS,J:=JJ(n+k)uh(n))|ΩeNS,L:=12JJ(n+k)uh(n)|ΩeAuh(n)|Ωe,zh(n)|ΩA^(n+k)uh(n)0e0,zh(n)0Auh(n)|Ω(e)A^(n+k)uh(n)0((e)0)).(3.5)

These additional terms are assumed to be negligibly small for k sufficiently large due to the fact that they consist of evaluations of differences GG(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 uh := uh(n), zh := zh(n), and

eQ:=eQ(n+k)=A^(uh)(zh)=0.

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

eϵ:=A^(n+k)uh(n)zh(n)A^(n+k)uh(n)0zh(n)0.(3.6)

### Proposition 3.1

Let u resp. z be the solution of the primal resp. dual problem in (2.1) resp. (3.1) with approximations u+, z+V+, where V+Vh. Then, for the perturbed discrete solutions uh(n) resp. zh(n) of (2.4) resp. (3.2), it holds that

J(u)J(n+k)uh(n))|Ω=eD(n+k)(u+|Ω,z+|Ω)+eQ(n+k)+eNS,Q+eNS,J+eNS,L+eNS,D(u+|Ω,z+|Ω)+eε+Rh(3)(3.7)

where the higher-order remainder Rh(3) is given by

Rh(3):=01Juh(n)|Ω+te(e,e,e)A(uh|Ω+te)(e,e,e,zh(n)|Ω+te)3Auh(n)|Ω+te(e,e,e)t(t1)dt.

### Proof

Let : ℝ → ℝ be defined as

(t):=L(γu(t),γz(t))=Luh(n)|Ω),zh(n)|Ω+t(e,e)

where γu(t):=uh(n)|Ω+te,γz(t):=zh(n)|Ω+te. Note that γu(t)=e,γz(t)=e. The derivative of is

(t)=L(γu(t),γz(t))(e,e)(1,1)T.

Applying the definition of 𝓛, we get

(t)=J(γu(t))(e)A(γu(t))(e,γz(t))A(γu(t))(e).

Applying differentiation twice more yields

(t)=J(γu(t))(e,e,e)A(γu(t))(e,e,e,γz(t))3A(γu(t))(e,e,e).(3.8)

The error using the exact functional J can be written in the following form:

J(u)Juh(n))|Ω=L(u,z)+A(u)(z)L(uh(n))|Ω,zh(n)|ΩAuh(n))|Ωzh(n)|Ω=L(u,z)L(uh(n))|Ω,zh(n)|ΩAuh(n))|Ωzh(n)|Ω=(1)(0)Auh(n))|Ωzh(n)|Ω=01(t)dtAuh(n))|Ωzh(n)|Ω.

It follows that

J(u)J(n+k)uh(n))|Ω=01(t)dt+eQ(n+k)+eε+eNS,Q+eNS,J.

We use the error representation of the trapezoidal rule to obtain

01(t)dt=12(0)+(1)+1201(t)t(t1)dt.

Since ″′(t) has been determined in (3.8), it remains to inspect the terms ′(0) and ′(1). We use that γu(0) = uh(n)|Ω,γz(0)=zh(n)|Ω, and γu(1) = u, γz(1) = z, so that

(1)=J(u)(e)A(u)(e,z)A(u)(e)=0

since (u, z) is a stationary point of 𝓛. We see that

(0)=Juh(n))|Ω(e)Auh(n))|Ωe,zh(n)|ΩAuh(n))|Ω(e)=ρ(e0)+ρ((e)0)+2eNS,L.

Thus, (0)=2eD(n+k)(u,z)+2eNS,L. The definition Rh(3)=1201(t)t(t1)dt and the calculation

eD(n+k)(u,z)=12ρz0zh(n)0+ρu0uh(n)0=12(ρ(z0(z+)0)+ρ(u0(u+)0))+12ρ(z+)0zh(n)0+ρ(u+)0uh(n)0=eNS,D(u+|Ω,z+|Ω)+eD(n+k)(u+|Ω,z+|Ω)

imply the proposed error identity. □

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.

### Corollary 3.1

Omitting all error terms which can be ensured to be negligibly small, we have the approximate error representation

errexcomp(n+k):=J(u)J(n+k)uh(n))|ΩeD(n+k)(u+|Ω,z+|Ω)+eQ(n+k)=:η(n+k).(3.9)

The effectivity index (or overestimation index) is defined as

eff(n+k):=η(n+k)comp(n+k).(3.10)

### 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 2d 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 V2h,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 V2h,p and V2h,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 Vh,p and V2h,2p have the same numbers of unknowns in the same Lagrange points, one can define an interpolation operator i* : Vh,pV2h,2p by assigning to vhVh,p the element i*vh of V2h,2p uniquely determined by the values of vh in those Lagrange points. The filtering operator is defined by π2h := id − i2h for a finite-element interpolation operator i2h : VV2h,p. The name of the method stems from the observation that π2hvh is a strictly local algebraic process acting on the coefficient vector v ∈ ℝN of vhVh,p, since π2hvh = vhi2hvh = j=1Nvj (φh,ji2hφh,j) = : ∑j (π2hv)jφh,j for the basis (φh,j)j=1N of Vh,p (see [35, eq. (40)]). Finally, the localization of eD(n+k) to nodal contributions is as follows:

eD(n+k)j=1N12A^(n+k)uh(n)0(iid)(ϕh,j)0(π2hz)j+J(n+k)uh(n)|Ω(iid)(ϕh,j)0(π2hu)j|ΩA^(n+k)uh(n)0(iid)(ϕh,j)0(π2hu)j,zh(n)0=:j=1Nη~D,j(3.11)

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

indD(n+k):=TThηD,T(n+k)|eD(n+k)|(3.12)

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

eQ(n+k)=TThηQ,T(n+k)

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

## 4 Refinement strategy

As identity (3.7) allows for a separation of the error into a term representing the discretization error eD(n+k) and a term representing the quadrature error eQ(n+k), we may perform finite-cell and quadrature mesh adaptation using these terms.

In the case of finite elements with exact quadrature, each step in the Solve–Estimate–Mark–Refine (SEMR) loop for adaptivity is well-examined at least for linear problems [10, 14]. However, for the finite-cell method, there are no strategies available for choosing the accuracy of the quadrature mesh. In practice, when e.g. the space-tree is used, a fixed number of recursive refinements throughout the computation is chosen. However, the chosen depth might be too high and a coarser integration mesh might be sufficient.

An issue that requires attention is the fact that, whenever an element T of the finite-cell mesh is refined and α(T) = 0, the quadrature mesh is refined in T as well, so that the quadrature error might decrease when only a decrease in the discretization error is intended. Hence, if α(T) > 0, it is reasonable to set α(T′) := α(T) − 1 for any child T′ of T in order not to introduce an additional improvement of the quadrature mesh where only an improvement in the finite-cell mesh is indicated. Furthermore, it is useful to introduce a lower bound l ∈ ℕ0 on α, e.g., αl := 0, when derefinements are not supported.

The adaptive strategy is summarized in the following steps. We emphasize that the error terms from (3.9) as well as the finite-cell mesh and the number of quadrature mesh refinements α now depend on the iteration index i of the SEMR loop. Also, to indicate the dependence of the approximate terms on α when the space-tree is used, we replace the generic number n indicating a sequence of approximations to the exact operators and functionals by the number of quadrature mesh refinements per element, given by the function α. Hence, we write α + k instead of n + k indicating that the quadrature mesh defined by α is refined k times globally.

1. Set i := 0. Initialize the finite-cell mesh 𝓣i. Choose an initial depth d ∈ ℕ0 and set the number of refinements of the quadrature scheme αi(T) := d for each T ∈ 𝓣i. Set l to be the minimum possible depth. Choose 0 < ρ ⩽ 1, e. g., ρ := 0.01. Choose a stopping criterion, e.g., stop if the maximum number of degrees of freedom is reached or if a prescribed error tolerance is met.

2. Construct the quadrature mesh for 𝓣i associated to αi.

3. Solve: Compute solutions uh,i(αi),zh,i(αi) of the perturbed discrete problems from (2.4), (3.2). Compute approximations ui+,zi+ to u, z.

4. Estimate: Choose k ∈ ℕ and construct a quadrature mesh on 𝓣i associated to αi + k to compute estimators eD,i(αi+k),eQ,i(αi+k) and indicators ηD,T,i(αi+k),ηQ,T,i(αi+k) for each T ∈ 𝓣i. If the stopping criterion is fulfilled, stop. If |eQ,i(αi+k)|ρ|eD,i(αi+k)|, localize the quadrature error contributions, mark the elements T with highest error contributions using a marking strategy (see step 5), set αi(T) := αi(T) + 1 for such elements, and go to step 2.

5. Mark: Choose an appropriate marking strategy, such as fixed-fraction marking or maximum marking [5], and mark elements with respect to the local discretization error ηD,T,i(n+k) for finite-cell mesh refinement.

6. Refine: Refine each marked element in 𝓣i to obtain 𝓣i+1. For the quadrature mesh, let αi+1 : 𝓣i+1 → ℕ0 and set αi+1|𝓣i∩𝓣i+1 := αi except for children T′ of any marked element T, where αi+1(T′) := αi(T) − 1 unless αi(T) − 1 < l.

7. Increase i by 1 and go to step 2.

## 5 Numerical results

In this section, we study the properties of the a posteriori error estimator and the adaptive algorithm by means of three examples which are characterized by some nonlinearities. In the first example, the primal and the dual solutions are known and smooth. Then we consider an example with a known but non-smooth solution. In the third example, the setting is more complex with respect to the geometry as well as the problem formulation. In this section, we omit indices such as (n + k) on the variables to improve the readability.

### 5.1 Quarter disk

The first domain is the quarter disk Ω := B1(0) ∩ (0, 1)2 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)2. We solve the non-linear diffusion–reaction equation − Δu + u3 = f, u|ΓD = 0, n u|ΓN = gN. Its weak form reads

a(u)(v)=Ωuv+u3v=F(v)=Ωfv+ΓNgv

with Gateaux derivative

a(u)(v,w)=Ωvw+3u2vw.

The functions f and gN 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

z(r,φ):=(r2+2r)16π2φ2+8πφ

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−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, eD, eQ, and eNS,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 eNS,Q, eNS,J, and eNS, 𝓛, which we want to assume as negligibly small, are decreasing with increasing k. We denote the mesh width of the quadrature mesh by hQ and note that an incrementation of k by 1 results in halving hQ. For the piecewise constant boundary approximation the terms eNS,Q, eNS,J, and eNS, 𝓛 are of order 𝓞(hQ) 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−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 eD are approximately constant from the beginning in the case of linear boundary approximation. The term eNS,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 eD become approximately constant for n ⩾ 10, while we do not see this for eNS,D here. The reason is that the values of eNS,D are polluted by the quadrature error even for n = 12. All other terms are of order 𝓞(hQ).

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−10 and a piecewise linear boundary approximation is used. We observe the optimal convergence in comp of order 𝓞(h2). The boundary approximation is accurate enough such that no further refinements of the quadrature are needed. We see higher order convergence for eNS,L, while the terms eQ, eNS,J, and eNS,Q are of order 𝓞(h2) but are essentially smaller than the error in J. The term in which we are mainly interested here is eNS,D measuring the error with respect to u+ and z+. This term is also of order 𝓞(h2) 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−4 the error is approximately constant as well as all other terms except for eQ, eNS,D and eε. The numerical error estimator becomes constant for ε < 10−6, which illustrates its stronger dependence on ε. This is also observed for eNS,D, which is constant for ε < 10−7. 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−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 ε

Fig. 7

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

Fig. 8

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

Tab. 1

Quarter disk: Adaptive iteration, number of degrees of freedom, range of quadrature levels, number of quadrature points and associated computational error, estimated discretization and quadrature error, and effectivity and indicator indices for a piecewise constant boundary approximation.

rDOFαQP|comp||eD||eQ|effindD
18103986.98⋅ 10−33.16⋅ 10−31.00⋅ 10−21.891.03
4810–38181.97⋅ 10−32.70⋅ 10−31.74⋅ 10−30.481.02
52531–31 5625.83⋅ 10−45.00⋅ 10−51.51⋅ 10−32.671.05
72531–53 1706.58⋅ 10−47.28⋅ 10−43.73⋅ 10−40.541.02
88211–45 1981.41⋅ 10−41.87⋅ 10−54.00⋅ 10−42.701.04
108212–610 8861.93⋅ 10−42.23⋅ 10−49.52⋅ 10−50.661.04
112 8331–519 3342.52⋅ 10−51.67⋅ 10−59.48⋅ 10−53.101.07
132 8333–741 6785.17⋅ 10−56.15⋅ 10−52.47⋅ 10−50.711.05
143 2412–741 9823.94⋅ 10−55.04⋅ 10−52.37⋅ 10−50.681.03
1511 9931–676 2481.06⋅ 10−51.84⋅ 10−62.48⋅ 10−52.161.05
1711 9933–8151 9687.14⋅ 10−61.12⋅ 10−58.72⋅ 10−60.351.04
1846 1892–7281 8305.71⋅ 10−61.24⋅ 10−69.00⋅ 10−61.791.05
2146 1895–10890 5662.22⋅ 10−63.29⋅ 10−61.82⋅ 10−60.661.04
2248 2055–10904 4782.18⋅ 10−62.99⋅ 10−61.71⋅ 10−60.591.02
23185 5334–91 422 8587.87⋅ 10−71.55⋅ 10−81.64⋅ 10−62.071.05
25185 5336–113 116 9704.40⋅ 10−77.50⋅ 10−75.32⋅ 10−70.501.04
26190 9895–113 100 0304.37⋅ 10−76.90⋅ 10−75.22⋅ 10−70.381.03
27728 2914–105 155 9762.98⋅ 10−74.78⋅ 10−85.05⋅ 10−71.851.06
29728 2916–1211 001 7049.24⋅ 10−81.83⋅ 10−71.51⋅ 10−70.341.06
30744 8995–1210 984 6429.22⋅ 10−81.66⋅ 10−71.48⋅ 10−70.191.02

Tab. 2

Quarter disk: Adaptive iteration, range of quadrature levels, number of quadrature points and associated computational error, estimated discretization and quadrature error, and effectivity and indicator indices for a piecewise constant boundary approximation for a fixed uniformly refined mesh with 16 641 degrees of freedom.

rαQP|comp||eD||eQ|effindD
1068 0787.57⋅ 10−44.45⋅ 10−47.04⋅ 10−41.51921.32
20–170 3823.70⋅ 10−42.12⋅ 10−43.51⋅ 10−41.52231.20
30–270 4421.77⋅ 10−49.64⋅ 10−51.75⋅ 10−41.53321.12
40–374 8348.60⋅ 10−54.20⋅ 10−59.21⋅ 10−51.55921.08
50–483 6784.22⋅ 10−51.58⋅ 10−55.22⋅ 10−51.61331.06
61–5101 4741.66⋅ 10−51.84⋅ 10−62.89⋅ 10−51.85451.04
72–6135 4582.68⋅ 10−65.59⋅ 10−61.64⋅ 10−54.01381.02
83–7204 7945.95⋅ 10−61.03⋅ 10−58.55⋅ 10−60.28501.02
94–8342 7461.05⋅ 10−51.32⋅ 10−54.48⋅ 10−60.83491.02
105–9618 9381.29⋅ 10−51.43⋅ 10−52.32⋅ 10−60.92971.02
116–101 167 0981.40⋅ 10−51.48⋅ 10−51.26⋅ 10−60.96211.02
127–112 259 1941.46⋅ 10−51.50⋅ 10−57.03⋅ 10−70.97821.02
138–124 447 2261.50⋅ 10−51.52⋅ 10−53.73⋅ 10−70.98801.02
149–138 823 2901.52⋅ 10−51.53⋅ 10−51.86⋅ 10−70.99331.02
1510–1417 556 9861.53⋅ 10−51.53⋅ 10−59.32⋅ 10−80.99571.02
1611–1535 024 3781.54⋅ 10−51.54⋅ 10−54.66⋅ 10−80.99691.02

Tab. 3

Quarter disk: Adaptive iteration, number of degrees of freedom, range of quadrature levels, number of quadrature points and associated computational error, estimated discretization and quadrature error, and effectivity and indicator indices for a piecewise linear boundary approximation.

rDOFαQP|comp||eD||eQ|effindD
18104015.08⋅ 10−34.65⋅ 10−31.09⋅ 10−31.131.03
22530–11 3231.84⋅ 10−31.45⋅ 10−37.03⋅ 10−41.171.01
382113 6413.44⋅ 10−43.12⋅ 10−44.21⋅ 10−51.031.03
42 833112 1099.48⋅ 10−58.54⋅ 10−51.43⋅ 10−51.051.06
53 229114 1777.55⋅ 10−57.06⋅ 10−58.66⋅ 10−61.051.02
611 917149 6672.00⋅ 10−51.85⋅ 10−52.89⋅ 10−61.071.04
712 389152 2611.82⋅ 10−51.73⋅ 10−51.76⋅ 10−61.051.01
846 9971192 4514.73⋅ 10−64.48⋅ 10−65.49⋅ 10−71.061.04
948 5331200 2714.39⋅ 10−64.22⋅ 10−63.17⋅ 10−71.031.02
10184 5491748 6111.12⋅ 10−61.08⋅ 10−68.07⋅ 10−81.041.05
11189 9571774 6191.06⋅ 10−61.04⋅ 10−64.06⋅ 10−81.011.02
12729 84512 943 0232.72⋅ 10−72.67⋅ 10−71.19⋅ 10−81.031.06
13747 53713 024 9912.63⋅ 10−72.60⋅ 10−77.48⋅ 10−91.021.02

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 14th 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 Ω := B1(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)2 ∖ ([0, 1] × [−1, 0]). The initial finite-cell mesh 𝓣0 consists of 3 ⋅ 16 square elements of degree 1. Again, we consider the diffusion–reaction equation −Δu + u3 = f, u|ΓD = 0, n u|ΓN = gN. The functions f and gN are chosen such that u(r, φ) = r2/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|2 at the point p=(cos54π,sin54π)Ω, which we approximate by J(v) := ∫Ω̂ĨS|∇ v|2, where ĨS with S = B0.01(p) ∩ Ω is a smoothed indicator function given on S by

I~S(x)=w(|xp|),xS,0,xΩ^S

where w(r) = −106r3 + 3 ⋅ 104r2 − 3 ⋅ 102r + 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−1) 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 13th iteration of adaptive algorithm for the piecewise linear boundary approximation.

Tab. 4

L-domain: Adaptive iteration, number of degrees of freedom, range of quadrature levels, number of quadrature points and associated computational error, estimated discretization and quadrature error, and effectivity and indicator indices for a piecewise constant boundary approximation.

rDOFαQP|comp||eD||eQ|effindD
122501 1947.15⋅ 10−52.40⋅ 10−41.33⋅ 10−4−1.506.74
112256–10161 5501.25⋅ 10−51.57⋅ 10−51.54⋅ 10−71.252.43
124276–10162 2421.44⋅ 10−62.68⋅ 10−82.32⋅ 10−70.1424.67
1942713–1718 777 4101.23⋅ 10−61.95⋅ 10−71.45⋅ 10−90.1632.46
2069713–1718 483 4341.39⋅ 10−61.05⋅ 10−61.70⋅ 10−90.7525.05
211 47312–1718 019 3289.20⋅ 10−77.60⋅ 10−79.48⋅ 10−100.835.33
224 31311–1718 300 4203.06⋅ 10−72.21⋅ 10−73.69⋅ 10−90.734.38
234 31312–1832 161 2843.04⋅ 10−72.21⋅ 10−71.75⋅ 10−90.734.39
246 34312–1831 603 5601.70⋅ 10−71.52⋅ 10−71.50⋅ 10−90.914.34
2521 54711–1729 032 9285.52⋅ 10−84.41⋅ 10−81.33⋅ 10−90.823.34
2621 54712–1852 183 5205.46⋅ 10−84.43⋅ 10−86.48⋅ 10−100.823.34
2721 54713–1999 172 8325.43⋅ 10−84.46⋅ 10−83.22⋅ 10−100.833.34
2825 69713–1894 223 8783.26⋅ 10−82.91⋅ 10−82.99⋅ 10−100.904.02

Tab. 5

L-domain: Adaptive iteration, number of degrees of freedom, range of quadrature levels, number of quadrature points and associated computational error, estimated discretization and quadrature error, and effectivity and indicator indices for a piecewise linear boundary approximation.

rDOFαQP|comp||eD||eQ|effindD
122501 2033.12⋅ 10−53.08⋅ 10−51.15⋅ 10−5−0.621.60
52250–44 6833.21⋅ 10−71.61⋅ 10−53.53⋅ 10−850.332.41
64270–45 4771.43⋅ 10−62.42⋅ 10−85.53⋅ 10−8−0.0241.57
84272–614 6291.48⋅ 10−62.19⋅ 10−77.16⋅ 10−100.1537.26
96972–615 5711.42⋅ 10−61.04⋅ 10−69.47⋅ 10−90.7224.04
101 4732–618 2019.20⋅ 10−77.36⋅ 10−76.91⋅ 10−100.805.33
114 3131–629 3072.78⋅ 10−71.95⋅ 10−72.48⋅ 10−80.614.48
134 3133–856 4753.02⋅ 10−72.20⋅ 10−79.35⋅ 10−100.724.41
146 3432–763 2871.68⋅ 10−71.52⋅ 10−75.25⋅ 10−100.904.34
1521 5471–7121 2895.37⋅ 10−84.40⋅ 10−83.26⋅ 10−100.813.34
1625 6691–7136 8433.25⋅ 10−82.88⋅ 10−81.19⋅ 10−100.893.98
1793 8491–6410 4031.04⋅ 10−87.87⋅ 10−96.97⋅ 10−100.824.02
1993 8493–8496 3719.70⋅ 10−98.46⋅ 10−95.27⋅ 10−110.884.02
20105 9632–8544 1736.67⋅ 10−96.08⋅ 10−95.11⋅ 10−110.924.38
21389 4191–71 667 9692.01⋅ 10−91.68⋅ 10−95.93⋅ 10−110.864.84
23389 4193–92 014 4491.95⋅ 10−91.64⋅ 10−95.03⋅ 10−120.844.84
24416 8912–82 115 6551.51⋅ 10−91.50⋅ 10−95.08⋅ 10−121.004.63
25495 9391–82 431 6351.12⋅ 10−91.14⋅ 10−95.19⋅ 10−121.025.41
261 596 0131–76 814 2673.64⋅ 10−104.56⋅ 10−105.14⋅ 10−121.265.26
271 596 0132–87 204 0593.62⋅ 10−103.81⋅ 10−102.82⋅ 10−121.065.26
281 736 8411–87 762 0073.14⋅ 10−103.39⋅ 10−102.79⋅ 10−121.096.60

### 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 8th iteration of adaptive algorithm.

div(C(u))=f,C(u)=u,|u|<σ0ζu+1ζ|u|u,else.(5.1)

Here, σ0 > 0 denotes the yield stress and 0 ⩽ζ ≪ 1 the hardening parameter. We choose σ0 = 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−10 as well as ε = 10−10. The right hand side f is given by a constant function with f≡ 2. The weak form of equation (5.1) reads

a(u)(v)=(C(u),v)+γΓFu3v=F(v)=(f,v)vV.

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

J(v)=SωDv2+ΓJωFv,S=(0.3,0.1)×(0.1,0.1)

where ωD(x, y) = 108 ⋅ (x + 0.3)2(x + 0.1)2(y + 0.1)2(y − 0.1)2 and ωF(φ, r) = φ2(φ − 0.5 π)2.

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 14th 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

J(u)4.3067544444371864104.

Fig. 14

Convergence results for the B-domain example.

Tab. 6

B-domain: Adaptive iteration, number of degrees of freedom, number of quadrature points, range of quadrature levels, estimated discretization and quadrature error, as well as effectivity and indicator indices for a piecewise linear boundary approximation.

rDOFαQP|comp||eD||eQ|effindD
115301 0281.99⋅ 10−64.91⋅ 10−63.96⋅ 10−72.261.14
21530–11 3642.19⋅ 10−64.98⋅ 10−62.42⋅ 10−72.161.15
31530–22 0282.39⋅ 10−65.10⋅ 10−66.11⋅ 10−82.111.15
41530–33 3242.43⋅ 10−65.13⋅ 10−61.58⋅ 10−82.101.15
53130–33 9522.10⋅ 10−61.09⋅ 10−61.82⋅ 10−80.512.09
63130–46 5042.12⋅ 10−61.11⋅ 10−64.62⋅ 10−90.522.09
78210–48 7502.02⋅ 10−66.87⋅ 10−74.42⋅ 10−90.341.49
81 7910–412 6861.05⋅ 10−63.93⋅ 10−71.87⋅ 10−90.371.43
93 1051–318 1424.73⋅ 10−71.55⋅ 10−73.97⋅ 10−90.321.49
103 1051–423 9504.76⋅ 10−71.58⋅ 10−71.00⋅ 10−90.331.49
116 1751–436 2002.45⋅ 10−79.03⋅ 10−89.68⋅ 10−100.361.43
126 1751–547 6242.46⋅ 10−79.09⋅ 10−82.82⋅ 10−100.371.43
1310 5691–564 9541.01⋅ 10−74.46⋅ 10−83.03⋅ 10−100.441.53
1420 4731–4104 2584.35⋅ 10−82.31⋅ 10−83.33⋅ 10−100.521.45
1520 4731–5127 6984.38⋅ 10−82.34⋅ 10−88.36⋅ 10−110.531.45
1636 0871–5189 1721.90⋅ 10−81.20⋅ 10−89.03⋅ 10−110.631.54
1754 0891–5260 4161.00⋅ 10−87.55⋅ 10−98.78⋅ 10−110.751.51
1854 0892–6313 2001.01⋅ 10−87.59⋅ 10−92.20⋅ 10−110.751.51
1992 3671–5462 9945.05⋅ 10−94.15⋅ 10−92.12⋅ 10−110.821.51
20173 8431–5783 6622.37⋅ 10−91.89⋅ 10−92.23⋅ 10−110.781.71
21173 8432–6883 0542.39⋅ 10−91.90⋅ 10−95.65⋅ 10−120.791.70
22347 2091–51 568 2601.15⋅ 10−99.77⋅ 10−105.77⋅ 10−120.841.59
23565 9291–52 432 7406.34⋅ 10−105.49⋅ 10−105.62⋅ 10−120.861.71
24565 9292–62 615 1566.39⋅ 10−105.52⋅ 10−101.54⋅ 10−120.861.71
251 230 7651–65 249 3062.96⋅ 10−102.49⋅ 10−101.51⋅ 10−120.841.72
262 046 3231–58 504 9681.76⋅ 10−101.39⋅ 10−101.52⋅ 10−120.781.81
272 046 3232–68 891 9441.77⋅ 10−101.40⋅ 10−104.10⋅ 10−130.791.81
282 661 7371–611 318 6821.37⋅ 10−101.10⋅ 10−104.15⋅ 10−130.801.72

### 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.

## 6 Conclusion

In this article, we present a dual weighted residual (DWR) error estimator for the finite cell method (FCM). The DWR method allows for goal-oriented error control and incorporates the information of a user-defined quantity of interest into the solution of a dual problem that has to be solved alongside the primal problem.

1. funding: The first and the last author gratefully acknowledge support by the Deutsche Forschungsgemeinschaft in the Priority Programme 1748 ‘Reliable simulation techniques in solid mechanics. Development of non-standard discretization methods, mechanical and mathematical analysis’ under the project ‘High-order immersed-boundary methods in solid mechanics for structures generated by additive processes', Grant SCHR 1244/4-1.

## References

[1] A. Abedian, J. Parvizian, A. Düster, and E. Rank, Finite cell method compared to h-version finite element method for elasto- plastic problems, Appl. Math. Mech. 35 (2014), No. 10, 1239-1248.10.1007/s10483-014-1861-9Search in Google Scholar

[2] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhäuser, 2013.Search in Google Scholar

[3] R. Becker, V. Heuveline, and R. Rannacher, An optimal control approach to adaptivity in computational fluid mechanics, Int. J. Numer. Methods Fluids40 (2002), No. 1-2, 105–120.10.1002/fld.269Search in Google Scholar

[4] R. Becker, C. Johnson, and R. Rannacher, Adaptive error control for multigrid finite element, Computing55 (1995), No. 4, 271–288.10.1007/BF02238483Search in Google Scholar

[5] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numerica10 (2001), 1–102.10.1017/S0962492901000010Search in Google Scholar

[6] H. Blum, A. Schröder, and F.-T. Suttmeier, A Posteriori Estimates for FE-Solutions of Variational Inequalities, pp. 669–680, Springer Milan, Milano, 2003.10.1007/978-88-470-2089-4_60Search in Google Scholar

[7] H. Blum and F.-T. Suttmeier, An adaptive finite element discretisation for a simplified Signorini problem, Calcolo37 (2000), No. 2, 65–77.10.1007/s100920070008Search in Google Scholar

[8] M. Braack and A. Ern, Adaptive Computation of Reactive Flows with Local Mesh Refinement and Model Adaptation, Numerical Mathematics and Advanced Applications, Springer, 2004, pp. 159–168.10.1007/978-3-642-18775-9_13Search in Google Scholar

[9] M. Braack and A. Ern, Coupling multimodelling with local mesh refinement for the numerical computation of laminar flames, Combustion Theory Model. 8 (2004), No. 4, 771–788.10.1088/1364-7830/8/4/006Search in Google Scholar

[10] D. Braess, Finite Elemente: Theorie, Schnelle Löser Und Anwendungen in Der Elastizitätstheorie, Springer-Verlag, 2013.10.1007/978-3-642-34797-9Search in Google Scholar

[11] A. Byfut and A. Schröder, A fictitious domain method for the simulation of thermoelastic deformations in NC-milling processes, Int. J. Numer. Methods Engrg. 113 (2017), No. 2, 208–229.10.1002/nme.5609Search in Google Scholar

[12] M. Dauge, A. Düster, and E. Rank, Theoretical and numerical investigation of the finite cell method, J. Sci. Comp. 65 (2015), No. 3, 1039–1064.10.1007/s10915-015-9997-3Search in Google Scholar

[13] M. De Berg, M. Van Kreveld, M. Overmars, and O. Schwarzkopf, Computational Geometry, Springer, 2000.10.1007/978-3-662-04245-8Search in Google Scholar

[14] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996), No. 3, 1106–1124.10.1137/0733054Search in Google Scholar

[15] A. Düster, J. Parvizian, Z. Yang, and E. Rank, The finite cell method for three-dimensional problems of solid mechanics, Comput. Methods Appl. Mech. Engrg. 197 (2008), No. 45, 3768–3782.10.1016/j.cma.2008.02.036Search in Google Scholar

[16] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Introduction to adaptive methods for differential equations, Acta Numerica4 (1995), 105–158.10.1017/S0962492900002531Search in Google Scholar

[17] R. Glowinski and Y. Kuznetsov, Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems, Comput. Methods Appl. Mech. Engrg. 196 (2007), No. 8, 1498–1506.10.1016/j.cma.2006.05.013Search in Google Scholar

[18] R. Glowinski, T.-W. Pan, and J. Periaux, A fictitious domain method for Dirichlet problem and applications, Comput. Methods Appl. Mech. Engrg. 111 (1994), No. 3-4, 283–303.10.1016/0045-7825(94)90135-XSearch in Google Scholar

[19] T. Grätsch and K.-J. Bathe, Goal-oriented error estimation in the analysis of fluid flows with structural interactions, Comput. Methods Appl. Mech. Engrg. 195 (2006), No. 41, 5673–5684.10.1016/j.cma.2005.10.020Search in Google Scholar

[20] S. Hubrich, P. Di Stolfo, L. Kudela, S. Kollmannsberger, E. Rank, A. Schröder, and A. Düster, Numerical integration of discontinuous functions: moment fitting and smart octree, Comp. Mech. 60 (2017), No. 5, 863–881.10.1007/s00466-017-1441-0Search in Google Scholar

[21] M. Joulaian, S. Duczek, U. Gabbert, and A. Düster, Finite and spectral cell method for wave propagation in heterogeneous materials, Comput. Mech. 54 (2014), No. 3, 661–675.10.1007/s00466-014-1019-zSearch in Google Scholar

[22] M. Joulaian, S. Hubrich, and A. Düster, Numerical integration of discontinuities on arbitrary domains based on moment fitting, Comp. Mech. 57 (2016), No. 6, 979–999.10.1007/s00466-016-1273-3Search in Google Scholar

[23] L. Kudela, N. Zander, T. Bog, S. Kollmannsberger, and E. Rank, Eflcient and accurate numerical quadrature for immersed boundary methods, Adv. Modeling Simul. Engrg. Sci. 2 (2015), No. 1, 10.10.1186/s40323-015-0031-ySearch in Google Scholar

[24] L. Kudela, N. Zander, S. Kollmannsberger, and E. Rank, Smart octrees: Accurately integrating discontinuous functions in 3D, Comput. Methods Appl. Mech. Engrg. 306 (2016), 406–426.10.1016/j.cma.2016.04.006Search in Google Scholar

[25] J. Nitsche, Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg36 (1971), No. 1, 9–15.10.1007/BF02995904Search in Google Scholar

[26] M. Paraschivoiu, J. Peraire, and A. Patera, A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations, Comput. Methods Appl. Mech. Engrg. 150 (1997), No. 1-4, 289–312.10.1016/S0045-7825(97)00086-8Search in Google Scholar

[27] J. Parvizian, A. Düster, and E. Rank, Finite cell method, Comput. Mech. 41 (2007), No. 1, 121–133.10.1007/s00466-007-0173-ySearch in Google Scholar

[28] J. Parvizian, A. Düster, and E. Rank, Topology optimization using the finite cell method, Optim. Engrg. 13 (2012), No. 1, 57–78.10.1007/s11081-011-9159-xSearch in Google Scholar

[29] S. Prudhomme and J. T. Oden, On goal-oriented error estimation for elliptic problems: application to the control of point- wise errors, Comput. Methods Appl. Mech. Engrg. 176 (1999), No. 1-4, 313–331.10.1016/S0045-7825(98)00343-0Search in Google Scholar

[30] A. Rademacher, Adaptive Finite Element Methods for Nonlinear Hyperbolic Problems of Second Order, Ph.D. thesis, TU Dortmund, 2010.Search in Google Scholar

[31] A. Rademacher and A. Schröder, Dual weighted residual error control for frictional contact problems., Comput. Meth. Appl. Math. 15 (2015), No. 3, 391–413.10.1515/cmam-2015-0014Search in Google Scholar

[32] E. Rank, S. Kollmannsberger, C. Sorger, and A. Düster, Shell finite cell method: a high order fictitious domain approach for thin-walled structures, Comput. Methods Appl. Mech. Engrg. 200 (2011), No. 45, 3200–3209.10.1016/j.cma.2011.06.005Search in Google Scholar

[33] R. Rannacher and J. Vihharev, Adaptive finite element analysis of nonlinear problems: balancing of discretization and iteration errors, J. Numer. Math. 21 (2013), No. 1, 23–62.10.1515/jnum-2013-0002Search in Google Scholar

[34] T. Richter, Goal-oriented error estimation for fluid-structure interaction problems, Comput. Methods Appl. Mech. Engrg. 223 (2012), 28–42.10.1016/j.cma.2012.02.014Search in Google Scholar

[35] T. Richter and T. Wick, Variational localizations of the dual weighted residual estimator, J. Comp. Appl. Math. 279 (2015), 192–208.10.1016/j.cam.2014.11.008Search in Google Scholar

[36] V. Sauľev, On the solution of some boundary value problems on high performance computers by fictitious domain method, Siberian Math. J. 4 (1963), No. 4, 912–925.Search in Google Scholar

[37] A. Schröder and A. Rademacher, Goal-oriented error control in adaptive mixed FEM for Signorini’s problem, Comput. Methods Appl. Mech. Engrg. 200 (2011), No. 1, 345–355.10.1016/j.cma.2010.08.015Search in Google Scholar

[38] F.T. Suttmeier, On plasticity with hardening: an adaptive finite element discretisation, Int. Math. Forum5 (2010), No. 52, 2591–2601.Search in Google Scholar

[39] K. van der Zee, E. van Brummelen, I. Akkerman, and R. de Borst, Goal-oriented error estimation and adaptivity for fluid- structure interaction using exact linearized adjoints, Comput. Methods Appl. Mech. Engrg. 200 (2011), No. 37, 2738–2757.10.1016/j.cma.2010.12.010Search in Google Scholar

[40] T. Wick, Goal-oriented mesh adaptivity for fluid-structure interaction with application to heart-valve settings, Arch. Mech. Engrg. 59 (2012), No. 1, 73–99.10.2478/v10180-012-0005-2Search in Google Scholar

[41] Z. Yang, S. Kollmannsberger, A. Düster, M. Ruess, E. Garcia, R. Burgkart, and E. Rank, Non-standard bone simulation: inter- active numerical analysis by computational steering, Comput. Visualization Sci. 14 (2011), No. 5, 207–216.10.1007/s00791-012-0175-ySearch in Google Scholar

[42] Z. Yang, M. Ruess, S. Kollmannsberger, A. Düster, and E. Rank, An eflcient integration technique for the voxel-based finite cell method, Int. J. Numer. Methods Engrg. 91 (2012), No. 5, 457–471.10.1002/nme.4269Search in Google Scholar

[43] N. Zander, T. Bog, M. Elhaddad, R. Espinoza, H. Hu, A. Joly, C. Wu, P. Zerbe, A. Düster, S. Kollmannsberger, J. Parvizian, M. Ruess, D. Schillinger, and E. Rank, FCMLab: A finite cell research toolbox for MATLAB, Adv. Engrg. Software74 (2014), 49–63.10.1016/j.advengsoft.2014.04.004Search in Google Scholar