Monique Dauge

Alexander Düster

Ernst Rank

Theoretical and numerical investigation of the finite cell method

We present a detailed analysis of the convergence properties of the finite cell method which is a fictitious domain approach based on high order finite elements. It is proved that exponential type of convergence can be obtained by the finite cell method for Laplace and Lamé problems in one, two as well as three dimensions. Several numerical examples in one and two dimensions including a well-known benchmark problem from linear elasticity confirm the results of the mathematical analysis of the finite cell method.

Introduction

The finite cell method (FCM) [START_REF] Parvizian | Finite cell method -h-and p-extension for embedded domain problems in solid mechanics[END_REF][START_REF] Üster | The finite cell method for three-dimensional problems of solid mechanics[END_REF][START_REF] Rank | Die Finite Cell Methode -Eine Fictitious Domain Methode mit Finiten-Element Ansätzen hoher Ordnung[END_REF] is a combination of a fictitious domain approach [START_REF] Saul'ev | A method for automatization of the solution of boundary value problems on high performance computers[END_REF][START_REF]On solution of some boundary value problems on high performance computers by fictitious domain method[END_REF][START_REF] Glowinski | A least squares/fictitious domain method for mixed problems and Neumann problems[END_REF][START_REF]A fictitious domain method for Dirichlet problem and applications[END_REF][START_REF] Glowinski | Distributed Lagrange multipliers based on fictitious domain method for second order elliptic problems[END_REF] with finite elements of high order [START_REF] Szab | Finite element analysis[END_REF][START_REF] Babu | The p and h-p versions of the finite element method, basic principles and properties[END_REF][START_REF]Approximation properties of the h-p version of the finite element method[END_REF][START_REF] Schwab | p-and hp-finite element methods, theory and applications in solid and fluid mechanics[END_REF][START_REF] Üster | The p-version of the finite element method for three-dimensional curved thin walled structures[END_REF] The main idea is to embed the domain of the problem to be solved into a bigger domain that has a simple geometric shape and can therefore be readily meshed. Thanks to the simple shape of the embedding or fictitious domain, mesh generation is dramatically simplified. The geometry of the problem is considered during the integration of the cell matrices, i.e. when computing the stiffness and mass matrices.

In this paper we will consider Neumann boundary conditions at the transition from the physical to the fictitious domain. Our fictitious domain approach relies on the introduction of an indicator function α which is equal to 1 inside the domain and 0 outside of the domain. In order to avoid conditioning problems, α is set to a very small value ε close to zero outside the domain. In this way the variational formulation is stabilized and the energy contribution of the fictitious domain is weakly penalized, shifting the effort of meshing towards the numerical integration of the cell matrices. Since the quality of the finite cell approximation strongly depends on the accuracy of the numerical integration, an adaptive quadrature scheme is applied to compute the stiffness and mass matrices of cells that are cut by the boundary of the domain or include holes. The adaptive integration can be carried out very generally by applying quadtree (in 2D) and octree (in 3D) space partitioning schemes in a fully automatic, errorcontrolled fashion [START_REF]Performance of different integration schemes in facing discontinuites in the finite cell method[END_REF]. Müller et al. [START_REF] Üller | Highly accurate surface and volume integration on implicit domains by means of moment-fitting[END_REF] proposed a promising approach that enables the numerical integration of functions that are at least partly defined by the zero iso-contour of a level set function. To this end a solution of a small linear system based on a simplified variant of the moment fitting equations 1 [START_REF] Wall | Quadrature schemes for arbitrary convex/concave volumes and integration of weak form in enriched partition of unity methods[END_REF] has to be performed to find the modified weights of Gaussian quadrature scheme. Summarizing, the finite cell method is based on three important ingredients: a fictitious domain approach, high order shape functions and an adaptive integration of the cell matrices. Combining these ingredients allows to achieve an exponential type of convergence when performing a p-extension of the trial and test functions of the cells.

The FCM has been applied to several problems like linear elasticity in 2D [START_REF] Parvizian | Finite cell method -h-and p-extension for embedded domain problems in solid mechanics[END_REF] and 3D [START_REF] Üster | The finite cell method for three-dimensional problems of solid mechanics[END_REF], to shell problems [START_REF] Rank | Shell Finite Cell Method: A High Order Fictitious Domain Approach for Thin-Walled Structures[END_REF] as well as to problems in biomechanics [START_REF] Üster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF][START_REF] Yang | Non-standard bone simulation: interactive numerical analysis by computational steering[END_REF][START_REF] Yang | An efficient integration technique for the voxel-based Finite Cell Method[END_REF] or wave propagation [START_REF] Duczek | Simulation of Lamb waves using the spectral cell method[END_REF][START_REF] Joulaian | Finite and spectral cell method for wave propagation in heterogeneous materials[END_REF][START_REF] Duczek | Numerical analysis of Lamb waves using the finite and spectral cell method[END_REF]. Nonlinear problems such as geometrically nonlinearity [START_REF] Schillinger | Small and large deformation analysis with the p-and B-spline versions of the finite cell method[END_REF] or elastoplasticity [START_REF] Abedian | The Finite Cell Method for Elasto-Plastic Problems[END_REF][START_REF] Abedian | The finite cell method for the J 2 flow theory of plasticity[END_REF] have been addressed as well. The FCM has also been successfully applied to the numerical homogenization of materials with complicated microstructures [START_REF] Üster | Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method[END_REF] or to topology optimization [START_REF] Üster | Topology optimization based on the finite cell method[END_REF][START_REF]Topology optimization using the finite cell method[END_REF] in structural mechanics. Instead of classical hierarchic shape functions [START_REF] Szab | Finite element analysis[END_REF] NURBS, which have become very popular thanks to the isogeometric analysis [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], can also be successfully used within the FCM, see [START_REF] Rank | An unfitted hp adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry[END_REF][START_REF] Rank | Geometric modeling, Isogeometric Analysis and the Finite Cell Method[END_REF]. Local refinement strategies have been also developed for the FCM and it turned out that the hp-d method [START_REF] Rank | Adaptive remeshing and h-p domain decomposition[END_REF][START_REF] Üster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF] presents a general framework for local improvement of accuracy within the FCM, see [START_REF] Schillinger | The hp-d-adaptive finite cell method for geometrically nonlinear problems of solid mechanics[END_REF][START_REF] Üster | Local enrichment of the finite cell method for problems with material interfaces[END_REF].

Despite the fact that the FCM has been numerically demonstrated to yield high convergence rates in many different problems, there is still a lack of a thoroughly mathematical analysis of its convergence properties. Thus, this paper is devoted to the analysis of the FCM, proving its capabilities of achieving exponential convergence under conditions, which are similar to those for the p-version of the finite element method. Here, as already mentioned, we will focus on Neumann boundary conditions at the transition from the physical to the fictitious domain. Dirichlet boundary conditions will be studied in future work.

The layout of the paper is as follows: In Section 2 the setting of the problem is defined and the convergence of the discrete and continuous problem with respect to the penalization parameter ε is presented. In Section 3 Céa and Strang lemmas are revisited to check that strict positiveness of the bilinear forms is not necessary. In Section 4 the convergence of the p-version of finite elements is addressed which is one of the main ingredients of the FCM. Numerical examples for 1D problems are presented in Section 5 and the observed exponential convergence is proved. In Section 6 a two-dimensional benchmark of linear elasticity is studied and it is demonstrated that the exponential convergence can be also obtained in 2D. Finally, a conclusion is drawn in Section 7.

Neumann condition obtained by penalization 2.1 The setting

We consider a connected bounded domain Ω ⊂ R n , n = 1, 2 or 3 with reasonable smoothness assumption (Ω can be regular or polyhedral with extension property across its boundary), see Figure 1. We assume that the boundary of Ω has at least two connected components, in such a way that Ω has a finite number of holes: The complement domain Ω ′ := R n \ Ω has one unbounded component Ω ′ 0 and a finite number of bounded connected components Ω ′ j , j = 1, . . . , J. We denote by H the hole

H = ∪ J j=1 Ω ′ j , (2.1)
by Γ be the boundary of Ω ′ 0 , by Σ the boundary of H and by D the domain with holes removed More specifically, we are given two differential bilinear forms b 0 and b 1 of degree 1 which we assume for simplicity to be real symmetric with constant coefficients (at this point we may consider systems -Lamé -as well). The simplest, typical, example is given by the gradient bilinear form

D = Ω ∪ H. (2.2) 2 Note that the boundary of D is Γ, the boundary of Ω is Σ ∪ Γ. Ω Ω ′ 1 Ω ′ 2 Ω ′ 3 Σ Σ Σ Γ Γ D = Ω ∪ H
b 0 (u, v)(x) = ∇u(x) • ∇v(x) + k 2 u(x)v(x) (k ≥ 0) and b 1 (u, v)(x) = ∇u(x) • ∇v(x), (2.3) while Lamé system corresponds to b 0 (u, v) = b 1 (u, v) = 2µǫ(u)(x) : ǫ(v)(x) + λ div u(x) div v(x).
Let us define variational spaces:

V (Ω) = {u ∈ H 1 (Ω) : u = 0 on Γ} and V (D) = H 1 0 (D), (2.4)
and the variational forms

a 0 (u, v) = Ω b 0 (u, v)(x) dx and a 1 (u, v) = H b 1 (u, v)(x) dx. (2.5) Let f ∈ L 2 (Ω).
We want to solve the variational mixed problem

Find u ∈ V (Ω), ∀v ∈ V (Ω), a 0 (u, v) = Ω f (x)v(x) dx. (2.6)
Instead, we solve for small ε > 0 the following variational problem on D

Find u ε ∈ V (D), ∀v ∈ V (D), a 0 (u ε , v) + ε a 1 (u ε , v) = D f (x)v(x) dx, (2.7)
where the right hand side f is extended by 0 in the hole H. Note that

a 0 (u, v) + ε a 1 (u, v) = Ω b 0 (u, v)(x) dx + ε H b 1 (u, v)(x) dx.
For our typical example (2.3), we have

a 0 (u, v) + ε a 1 (u, v) = D α(x) ∇u(x) • ∇v(x) + α M (x) u(x)v(x) dx, (2.8)
where the functions α and α M are defined as follows

α(x) = 1 and α M (x) = k 2 in Ω, α(x) = ε and α M (x) = 0 in H. (2.9)
At the discrete level, we replace a FEM discretization of problem (2.6) by a discretization of (2.7): With a finite dimensional approximation V ap (D) of V (D), and a numerical integration D|ap over D, the discrete problem with parameter ε is

Find u ap ε ∈ V ap (D), ∀v ∈ V ap (D), a ap 0 (u ap ε , v) + ε a ap 1 (u ap ε , v) = D|ap f (x)v(x) dx, (2.10)
where

a ap 0 (u, v) = D|ap 1 Ω (x) b 0 (u, v)(x) dx and a ap 1 (u, v) = D|ap 1 H (x) b 1 (u, v)(x) dx. (2.11)
We are going to prove that solutions of problems (2.7) and (2.10) are the sums of convergent power series in ε. This helps to understand the structure of these solutions. However, our analysis of the FCM does not rely on such power series expansions, but on an extended Strang lemma stated in Section 3.

Convergence of discrete problems with respect to the penalization parameter

Let us consider the finite dimensional space V = V ap (D) and the numerical integration D|ap as fixed, and let ε tend to 0. Later on, we will assume that the numerical integration is sufficiently accurate. The influence of the integration error will be studied in Section 5.4.5. However, for the following abstract setting, we do not need any special assumption on the numerical integration.

We prove in Lemma 2.1 that, under a simple assumption, problem (2.10) converges to a limit as ε → 0.

We introduce the kernel of a ap 0

K 0 = {v ∈ V : ∀u ∈ V, a ap 0 (u, v) = 0}, (2.12)
and its orthogonal space

K ⊥ 0 = {ϕ ∈ V ′ : ∀v ∈ K 0 , ϕ, v = 0}. (2.13)
Here ϕ, v denotes the duality pairing between V ′ and V .

We define the operators Āℓ for ℓ = 0, 1:

Āℓ : V -→ V ′ u -→ V ∋ v → a ap ℓ (u, v) , (2.14)
and introduce their restrictions

A 0 : V -→ K ⊥ 0 u -→ V ∋ v → a ap 0 (u, v) , (2.15)
and

A 1 : V -→ K ′ 0 u -→ K 0 ∋ v → a ap 1 (u, v) , (2.16)
and, finally, the operator A A :

V -→ K ⊥ 0 × K ′ 0 u -→ (A 0 u, A 1 u).
(2.17)

Lemma 2.1 In the finite dimensional framework, if A 1 K 0 is bijective, then A is bijective. Let f ∈ K ⊥ 0 .
Then for ε small enough, problem (2.10) has a unique solution u ε . Let u 0 be defined as A -1 (f, 0). Then u ε tends to u 0 in V , and in any fixed norm

• on V u ε -u 0 ≤ C A ε f . Proof: A 0 sends V into K ⊥ 0 , and A sends V into K ⊥ 0 × K ′ 0 .
These two latter spaces have the same dimension. Since A 1 K 0 is bijective, the kernel of A is reduced to {0}, hence the bijectivity of A.

We define recursively: u (0) := u 0 = A -1 (f, 0) and u (j) = -A -1 ( Ā1 u (j-1) , 0), j = 1, 2, . . . This makes sense, since by definition A 1 u (j-1) = 0, hence Ā1 u (j-1) belongs to K ⊥ 0 . Then, for ε ≤ ε 0 with ε 0 small enough, the series j≥0 ε j u (j) converges in V and is a solution of problem (2.10).

More generally, if the right-hand side is any element ϕ of V ′ , we solve the problem

Find u ∈ V, ∀v ∈ V, a ap 0 (u, v) + ε a ap 1 (u, v) = ϕ, v , (2.18)
by the series j≥-1

ε j u (j)
with u (-1) ∈ K 0 the unique solution of 1) , 0), and u (j) for j ≥ 1 defined as above.

∀v ∈ K 0 , a ap 1 (u, v) = ϕ, v , u (0) := A -1 (f -Ā1 u (-
This proves that the operator

V ∋ u -→ v → a ap 0 (u, v) + ε a ap 1 (u, v) ∈ V ′ is
onto as soon as ε ≤ ε 0 . Therefore, it is injective. This ends the proof.

Convergence of the continuous problem with respect to the penalization parameter

Let us recall that a real symmetric bilinear form a is said to be coercive on a space V endowed with a norm • if for some positive constant c, there holds a(u, u) ≥ c u 2 for all u ∈ V .

We assume that a 0 is coercive on V (Ω). Then the kernel

K 0 := {v ∈ V (D) : ∀u ∈ V (D), a 0 (u, v) = 0}, is by definition the space of v ∈ H 1 0 (D) such that for all u ∈ H 1 0 (D), Ω b 0 (u, v)(x) dx = 0.
Since any function u ∈ V (Ω) can be extended in a function ū ∈ H1 0 (D), and, conversely, the restriction of any v ∈ H 1 0 (D) to Ω is an element of V (Ω), the coercivity property of a 0 implies that

K 0 = {v ∈ H 1 0 (D) : v Ω = 0}
. Hence K 0 is the space of the extensions by zero to Ω of the elements of H 1 0 (H). Thus the orthogonal space is

K ⊥ 0 = {ϕ ∈ H -1 (D) : ∀v ∈ K 0 , ϕ, v = 0}.
It is the space of the extensions by zero 1 to H of the elements of V (Ω) ′ .

We assume that a 1 is coercive on H 1 0 (H) and define the operator A like in (2.15)-(2.17

). Let f ∈ K ⊥ 0 . The function u (0) = u 0 := A -1 (f, 0) is by definition the solution of Find u ∈ H 1 0 (D) such that ∀v ∈ H 1 0 (D), a 0 (u, v) = f, v and ∀v ∈ H 1 0 (H), a 1 (u, v) = 0.
(2.19) Let for k = 0, 1 the interior and boundary operators L k and B k be such that

a 0 (u, v) = - Ω L 0 u v dx + Σ B 0 u v dσ, ∀u, v ∈ V (Ω) such that L 0 u ∈ L 2 (Ω), (2.20) a 1 (u, v) = - H L 1 u v dx + Σ B 1 u v dσ, ∀u, v ∈ H 1 (H) such that L 1 u ∈ L 2 (H). (2.21)
We can see that u 0 Ω is the solution u + 0 of the mixed Dirichlet (on Γ) Neumann (on Σ) problem associated with a 0 on Ω, with right-hand side f , and u 0 H is the solution u - 0 of the Dirichlet problem

L 1 u - 0 = 0 in H and u - 0 Σ = u + 0 Σ .
Note that the next term u (1) := -A -1 ( Ā1 u (0) , 0) as defined in the proof of Lemma 2.1 has the following structure. Let u + 1 and u - 1 its restriction to Ω and H, respectively. Then u + 1 is solution of the mixed problem

L 0 u + 1 = 0 in Ω, u + 1 Γ = 0, and B 0 u + 1 Σ = B 1 u + 0 Σ , and u - 1 is the solution of the Dirichlet problem L 1 u - 1 = 0 in H and u - 1 Σ = u + 1 Σ .
We have a statement similar to Lemma 2.1

Lemma 2.2 In the continuous framework, let a 0 be coercive on V (Ω) and a 1 be coercive on

H 1 0 (H). Let f ∈ V (Ω) ′ . Then for ε small enough, the problem Find u ε ∈ V (D), ∀v ∈ V (D), a 0 (u ε , v) + ε a 1 (u ε , v) = f, v
has a unique solution. Let u 0 be the solution of (2.19). Then u ε tends to u 0 in H 1 (D), and

u ε -u 0 H 1 (D) ≤ Cε f V (Ω) ′ .
Proof: It follows the same lines as the proof of Lemma 2.1. Since

ϕ → A -1 (ϕ, 0) is continuous from V (Ω) ′ into H 1 (D), and u (j-1) → Ā1 u (j-1) is continuous from H 1 (D) into V (Ω) ′ , we have an estimate u (j) H 1 (D) ≤ C u (j-1) H 1 (D) , ∀j ≥ 1.
We deduce the convergence in H 1 (D) of the series j≥0 ε j u (j) and the estimate of the Lemma.

Céa and Strang lemmas

Céa and Strang lemmas are classical cornerstones of FEM analysis, see the monograph [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. But in general, to our knowledge, they use as an assumption that the bilinear forms involved are coercive. Here we show that we may in a certain amount, replace the assumption of positiveness by a non-negativeness assumption. Though the principles of these proofs are standard, we did not find them in such a general framework in the literature, so we prefer to give them for the sake of completeness.

Lemma 3.1 (Céa) Let a be a real symmetric bilinear form, non-negative on the space V :

∀u ∈ V, a(u, u) ≥ 0,

defining the semi-norm |u| a := a(u, u) 1 2
.

Let V ap be a subspace of V . Let f ∈ V ′ . We assume that u ∈ V and u ap ∈ V ap satisfy a(u, v) = f, v ∀v ∈ V and a(u ap , v ap ) = f, v ap ∀v ap ∈ V ap . Then |u -u ap | a ≤ |u -v ap | a ∀v ap ∈ V ap . (3.1) Proof: We have for all v ap ∈ V ap a(u -u ap , u -u ap ) = a(u -u ap , u -v ap + v ap -u ap ) = a(u -u ap , u -v ap ).

Inequality (3.1) then follows by Cauchy-Schwarz inequality.

The lemma that we present here is known as first Strang lemma [10, Thm. 4.1.1] when assorted with the usual coercivity assumptions. We may also refer to original papers by Strang himself [START_REF] Strang | The mathematical foundations of the finite element method with applications to partial differential equations[END_REF][START_REF]Piecewise polynomials and the finite element method[END_REF] where the bases of convergence analysis are laid when variational crimes are committed, see also [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]Chap. 10].

Lemma 3.2 (Strang)

Let a be a real symmetric bilinear form, non-negative on the space V . Let V ap be a subspace of V and let a ap be a real symmetric bilinear form, non-negative on V ap . Let d(f, v) be a duality pairing between V ′ and V , and d ap be a duality pairing between V ′ and V ap . We define the semi-norms |u| a := a(u, u)

1 2
and |u| a ap := a ap (u, u)

1 2 .
We assume that there exists a positive constant C ap such that

|v ap | a ≤ C ap |v ap | a ap ∀v ap ∈ V ap . (3.2)
We assume that u ∈ V and u ap ∈ V ap satisfy, for a given

f ∈ V ′ a(u, v) = d(f, v) ∀v ∈ V and a ap (u ap , v ap ) = d ap (f, v ap ) ∀v ap ∈ V ap .
Then for all v ap ∈ V ap the following two inequalities hold:

|u -u ap | a ≤ (1 + C 2 ap )|u -v ap | a + C 2 ap sup w ap ∈V ap (a -a ap )(v ap , w ap ) |w ap | a + sup w ap ∈V ap (d -d ap )(f, w ap ) |w ap | a (3.3) and |u -u ap | a ≤ (1 + C 2 ap )|u -v ap | a + C ap sup w ap ∈V ap (a -a ap )(v ap , w ap ) |w ap | a ap + sup w ap ∈V ap (d -d ap )(f, w ap ) |w ap | a ap (3.4)
Proof: Let us choose v ap ∈ V ap . We write

|u -u ap | a ≤ |u -v ap | a + |v ap -u ap | a . (3.5)
We set w ap = u apv ap . Then we evaluate |v apu ap | 2 a ap = a ap (v apu ap , v apu ap ):

|v ap -u ap | 2 a ap = a ap (u ap , w ap ) -a(v ap , w ap ) + (a -a ap )(v ap , w ap ). (3.6)
We note that a ap (u ap , w ap ) = d ap (f,

w ap ) = d(f, w ap ) -(d -d ap )(f, w ap ). thus a ap (u ap , w ap ) = a(u, w ap ) -(d -d ap )(f, w ap ). (3.7)
Combining (3.6) and (3.7):

|v ap -u ap | 2 a ap = a(u, w ap ) -a(v ap , w ap ) -(d -d ap )(f, w ap ) + (a -a ap )(v ap , w ap ). (3.8)
Identity (3.8) implies the inequality

|v ap -u ap | a ap |w ap | a ap ≤ |u -v ap | a |w ap | a + |(d -d ap )(f, w ap )| + |(a -a ap )(v ap , w ap )| . (3.9)
Combining (3.9) with (3.2): We can use Lemma 3.2 with a = a 0 or a = a 0 + εa 1 , and also take numerical integration into account in a ap and d ap .

|v ap -u ap | a |w ap | a ≤ C
Assuming exact integration, we can also use the lemma with a = a 0 , a ap = a 0 + εa 1 and d = d ap .

In this case, the third term in the right-hand side of (3.4) is zero and the second one is the sup for

w ap ∈ V ap of (a -a ap )(v ap , w ap ) |w ap | a ap = ε a 1 (v ap , w ap ) |w ap | a ap ≤ ε |v ap | a 1 |w ap | a 1 |w ap | 2 a 0 + ε|w ap | 2 a 1 1/2 ≤ √ ε |v ap | a 1 .
In this case, C ap = 1 and (3.4) yields:

Corollary 3.
3 Under the conditions of Lemma 3.2 with a = a 0 and a ap = a 0 + εa 1 , we assume moreover exact integration (d = d ap ). Then we have the estimate

|u -u ap ε | a 0 ≤ 2|u -v ap | a 0 + √ ε |v ap | a 1 ∀v ap ∈ V ap . (3.11)
Remark 3.4 Under the conditions of Lemma 3.1 with a = a 0 and Lemma 2.1 with a ap 0 = a 0 and a ap 1 = a 1 (i.e., assuming exact integration), we can write for all ε

|u -u ap ε | a 0 ≤ |u -u ap 0 | a 0 + |u ap 0 -u ap ε | a 0 . Lemma 3.

p-version of finite elements

Let T be a fixed mesh of the domain D by segments, quadrilateral or hexahedral elements according to the space dimension n. Let V p (D) be the p-extension over the mesh T with the boundary condition v = 0 on Γ and C 0 conformity between elements. By "p-extension" we mean that on each element K ∈ T , the discrete functions are mapped polynomials of partial degree p on the reference element K. We denote by

Ω ap = interior K∈T | K∩Ω = ∅ K and H ap = D \ Ω ap . (4.1)
We note that Γ is contained in ∂Ω ap and we denote the common boundary ∂Ω ap ∩ ∂H ap by Γ ap .

We denote by A(U) the space of analytic functions up to the boundary of the domain U .

Theorem 4.1 Let a 0 be coercive on V (Ω) and a 1 be coercive on H 1 0 (H). Let f ∈ A(Ω). Let u 0 be the solution of the mixed problem (2.6). We assume that u 0 admits an analytic extension ū0 ∈ A(Ω ap ). For ε > 0 and any p ≥ 1, let u ap [p, ε] be the solution of problem (2.10) with V ap (D) = V p and assuming exact integration. Then there exist c > 0 and γ > 0 such that for all ε > 0 small enough and all p ≥ 1

u 0 -u ap [p, ε] H 1 (Ω) ≤ c(e -pγ + √ ε). (4.2)
Proof: We use Corollary 3.3. By the coercivity assumption on a 0 , we find that the semi-norm | • | a 0 is equivalent to the H 1 (Ω)-norm. Thus, relying on estimate (3.11), it suffices to find v ap ∈ V p such that

u 0 -v ap H 1 (Ω) ≤ ce -pγ and v ap H 1 (H) ≤ c. ( 4
∈ V p (Ω ap ) such that ū0 -v p H 1 (Ω ap ) ≤ c 0 e -pγ . (4.4)
with a similar control of the elementwise tensor H 1 -norm. Therefore, in particular,

v p H 1 (Γ ap ) ≤ c 1 (4.5)
with a similar control of nodal values and, if n = 3, of H 1 -norm on edges. Then using tensor trace liftings as mentioned above, we find that there exists ṽp ∈ V p (H ap ) such that

ṽp Σ = v p Σ and ṽp H 1 (H ap ) ≤ c 2 . (4.6)
We define v ap by v p on Ω ap and ṽp on H ap and we deduce (4.3) from (4.4)-(4.6).

In the case where the grid is matching with the interface Σ, the estimate (4.2) is improved (see also Remark 3.4 on this question -why such an improvement does not hold in the general case).

Theorem 4.2 Under the assumptions of Theorem 4.1 we assume moreover that Ω ap = Ω. Then there exist c > 0 and γ > 0 such that for all ε > 0 small enough and all p ≥ 1

u 0 -u ap [p, ε] H 1 (Ω) ≤ c(e -pγ + ε). (4.7)
2 Let Î = (-1, 1) be the reference segment. Tensor H 1 -norms are defined as follows in dimension n = 2 and n = 3:

u 2 H 1,1 ( Î2 ) = 1 α 1 =0 1 α 2 =0 ∂ α 1 x 1 ∂ α 2 x 2 u 2 L 2 ( Î2 )
and

u 2 H 1,1,1 ( Î3 ) = 1 α 1 =0 1 α 2 =0 1 α 3 =0 ∂ α 1 x 1 ∂ α 2 x 2 ∂ α 3 x 3 u 2 L 2 ( Î3 )
Proof: Let us choose the degree p. When the grid is matching the interface Σ, the terms u (j) [p] of the expansion of u ap [p, ε] in powers of ε can be described as discrete FEM solutions: Let u + j [p] and u - j [p] be the restrictions of u (j) [p] to Ω = Ω ap and H = H ap , respectively. For j ≥ 1, u + j [p] is the discrete solution of the mixed problem in Ω with Neumann data on Σ coming from u - j-1 [p], and u - j [p] is the discrete solution of the Dirichlet problem in H with Dirichlet data on Σ coming from u + j [p]. The uniformity of continuity constants with respect to p can be deduced.

If Ω ap coincides with D, we may even have exponential convergence with ε = 0: Theorem 4. [START_REF] Abedian | The finite cell method for the J 2 flow theory of plasticity[END_REF] We assume that Ω ap = D (i.e., H ap = ∅). Let a 0 be coercive on V (Ω). Let f ∈ A(Ω) and let u 0 be the solution of the mixed problem (2.6). We assume that u 0 admits an analytic extension ū0 ∈ A(Ω ap ). Let u ap [p] be solution of (here we assume exact integration)

Find u ap ∈ V p (D), ∀v ∈ V p (D), a 0 (u ap , v) = D f (x)v(x) dx, (4.8)
Then there exist c > 0 and γ > 0 such that for all p ≥ 1

u 0 -u ap [p] H 1 (Ω) ≤ c e -pγ . (4.9)
Proof: This is a consequence of Céa Lemma 3.1 if we have proved that u ap [p] does exist. Let us choose the degree p. It suffices to show that the kernel K 0 defined as

K 0 = {v ∈ V p (D) : ∀u ∈ V p (D), a 0 (u, v) = 0} is reduced to {0}. Let v ∈ K 0 . Then a 0 (v, v) = 0. Since v Ω belongs to V (Ω)
, we deduce from the coercivity property of a 0 that v Ω ≡ 0. By assumption any element K of the mesh T has a non-empty intersection with Ω. Since v is a polynomial on K which is zero on K ∩ Ω, it is zero over the whole of K. Hence v ≡ 0, which ends the proof. The important point for the analysis of Theorem 4.1-4.3 to hold is that the boundary of H is analytic, so that the analytic continuation of solutions across Γ is possible. Near corners and edges, solutions with analytic data belong to weighted analytic spaces (also called countably normed spaces), see [START_REF] Babu | Regularity of the solution of elliptic problems with piecewise analytic data. I. Boundary value problems for linear elliptic equation of second order[END_REF][START_REF]Regularity of the solution of elliptic problems with piecewise analytic data. II. The trace spaces and application to the boundary value problems with nonhomogeneous boundary conditions[END_REF] when n = 2, [START_REF] Guo | Regularity of the solutions for elliptic problems on nonsmooth domains in R 3 . I. Countably normed spaces on polyhedral domains[END_REF][START_REF]Regularity of the solutions for elliptic problems on nonsmooth domains in R 3 . II. Regularity in neighbourhoods of edges[END_REF] for preliminary results and [START_REF] Costabel | Analytic regularity for linear elliptic systems in polygons and polyhedra[END_REF][START_REF]Weighted analytic regularity in polyhedra[END_REF] for complete results when n = 3. Nevertheless, we will see in Section 6 an example where these assumptions are not satisfied, i.e., a mere p-extension is used for the Lamé system on a square, and where pre-asymptotic exponential convergence for ε = 0 can be observed in an error range which is relevant for engineering practice. △ Remark 4.5 The polynomial recovery property of finite element approximations is only fulfilled in the limit case ε = 0, as for all finite ε a model error is included and the modified problem will in general not have the polynomial as exact solution. Yet, in all numerical tests using a sufficiently small ε, deviations from a polynomial solution were negligible. △

In Theorems 4.1-4.3, we have assumed exact numerical integration. In the FCM, the numerical integration is based on a subpartition of cells, the sub-cells, which are geometrically refined near the boundary Γ, see Figure 16. The numerical integration is exact, except on a set of very fine sub-cells covering Γ; This latter contribution to the error can be kept exponentially small with respect to p, see Figure 17.

The capabilities of hp quadrature are illustrated by the paper [START_REF] Chernov | Exponential convergence of hp quadrature for integral operators with Gevrey kernels[END_REF].

Finite cell approach

The bilinear form (5.3) as described in the previous subsection is discretized by means of the finite cell method. To this end, the fictitious domain D is subdivided into a mesh T consisting of n c cells with the corresponding nodal coordinates denoted as X c , X c+1 with 1 ≤ c ≤ n c . On each cell hierarchic shape functions N i based on integrated Legendre polynomials [START_REF] Szab | Finite element analysis[END_REF][START_REF] Üster | The p-version of the finite element method for three-dimensional curved thin walled structures[END_REF] are applied to discretize the trial and test functions. The discretization of the bilinear form results in a matrix composed of two parts: the stiffness matrix and the mass matrix. Since in general the cells do not conform with the geometry, the integrand of the cell stiffness matrix

K c ij = X c+1 Xc α dN i dx dN j dx dx = 2 X c+1 -X c 1 -1 α dN i d ξ dN j dξ dξ , i, j = 1, 2, 3, ... (5.8)
and the cell mass matrix

M c ij = X c+1 Xc α M N i N j dx = X c+1 -X c 2 1 -1
α M N i N j dξ , i, j = 1, 2, 3, ... (5.9) might be discontinuous. In (5.8) and (5.9) x, ξ denote the global and local coordinates, which are related to each other by a linear mapping function. Note, that we distinguish between α and α M that are defined by (2.9). Whereas α jumps from 1 in Ω to ε in H, α M jumps from k 2 to 0 (inside the hole, we set the mass matrix to 0). The integration of the cell matrices is carried out by applying a composed Gaussian quadrature. To account for the hole, i.e. the jump of α, α M , the corresponding cell is divided for the purpose of (exact) integration into n sc sub-cells, so that on each sub-cell α and α M are constant. In this way it is possible to perform an exact computation of the stiffness and mass matrix with n G = p + 1 Gaussian points applied on the sub-cells which are introduced just for integration purposes. Considering a mesh with one cell only, the minimum number of sub-cells needed for an exact integration is n sc = 3.

Evaluation of the error

In order to quantify the efficiency and accuracy of the finite cell method we briefly present in this section the definition of the error. Thanks to availability of the exact solution, the error e = uu ap (5.10)

of the finite cell approximation can be evaluated directly. In the following we compute the error in the 5.11) by considering the contribution in the domain Ω only, i.e. ignoring the results of the finite cell method in the hole H. Since the computation of the error in H 1 norm (5.11) involves the integration of nonpolynomials a Gaussian quadrature will not yield exact values. Therefore we apply an composite Gaussian quadrature as described in the previous section in order to reliably determine the error.

H 1 norm e 2 H 1 = -1/4 -1 (e ′ 2 + k 2 e 2 ) dx + 1 1/4 (e ′ 2 + k 2 e 2 ) dx (

Numerical examples for Neumann boundary conditions at the hole

In the following we present several numerical results obtained with the finite cell method discretizing the problem described in Section 5.1. We choose k = 3 and compute the error in terms of Equation (5.11).

Non-matching grid with one cell

First, we choose one finite cell with nodal coordinates X 1 = -1 and X 2 = 1 to discretize the fictitious domain and perform a p-extension with p = 1, 2, 3, ..., 20. In this example, α and likewise α M is set to 0 inside the hole. The integration of the stiffness and mass matrix is carried out exactly. A comparison of the exact solution and the finite cell approximation for p = 20 is given in Figure 2. It can be seen that one cell very accurately represents the exact solution. Note that in the hole the FCM approximation presents a smooth behavior, connecting the two branches of the exact solution. To quantify the efficiency more precisely, the error e 2 H 1 of the FCM approximation with p = 1, 2, 3, ... is plotted in Figure 3 against the polynomial degree. From this it is evident that an exponential convergence can be obtained although the mesh consisting of one cell only does not conform to the geometry. It is also noted that the convergence of the problem with the even solution is faster. Based on the Céa Lemma 3.1 used with the form a 0 , we prove this exponential convergence if we know the existence of a polynomial v p of degree ≤ p such that

exact solution approximation -1 +1 x u(x) exact solution approximation -1 +1 x u(x)
u -v p H 1 (-1,-1 4 ) + u -v p H 1 ( 1 4 ,1) ≤ c e -pγ .
(5.12)

Lemma 5.1 Let λ ∈ (0, 1). Let g be a function defined and analytic on the union of intervals [-1, -λ]∪ [λ, 1]. There exists c > 0 and γ > 0 and for all p ≥ 1 a polynomial v p ∈ P p (-1, 1) such that

g -v p H 1 (-1,-1 4 ) + g -v p H 1 ( 1 4 ,1) ≤ c e -pγ .
(5.13)

Proof: Considering the even and odd parts of g, we may reduce to the case when g is either even or odd.

Even case: By the formula G(t) = g( √ t) we define an analytic function

G ∈ A[ √ λ, 1] such that ∀x ∈ [-1, -λ] ∪ [λ, 1], g(x) = G(x 2 ).
Let p = 2q be a positive (even) integer. By [START_REF] Schwab | p-and hp-finite element methods, theory and applications in solid and fluid mechanics[END_REF]Thm. 3.20], there exists

Φ q ∈ P q ( √ λ, 1) satisfying the estimate G -Φ q H 1 ( √ λ,1) ≤ c e -qγ ′ . (5.14)
We set v p (x) = Φ q (x 2 ) and have proved (5.13) with γ = γ ′ /2.

Odd case: Similarly, by the formula

G(t) = t -1/2 g( √ t) we define an analytic function G ∈ A[ √ λ, 1] such that ∀x ∈ [-1, -λ] ∪ [λ, 1], g(x) = xG(x 2 ).
Let p = 2q + 1 be a positive (odd) integer. Like above, there exists Φ q ∈ P q ( √ λ, 1) satisfying the estimate (5.14). We set v p (x) = xΦ q (x 2 ) and have proved (5.13) as before.

Remark 5.2 The steps in both curves in Figure 3 originate from the improvement of convergence when the parities of the solution and the polynomials coincide, since only one element is used. Both curves indicate an exponential convergence, involving however a larger prefactor in the odd case. This is due to distinct approximation properties of the auxiliary functions G in the even and odd case (G(t) = g( √ t) and G(t) = t -1/2 g( √ t), respectively). △

Matching grid with three cells

Next, we consider the odd problem discretized by a mesh with three cells where the layout is such that the hole is precisely covered by one cell. The nodal coordinates X c of the mesh correspond to {-1, -0.25, 0.25, 1}. In this case, the integration of the cell matrices can be carried out exactly by a standard Gaussian quadrature with n G = p + 1 without the necessity of introducing sub-cells. The aim of this example is to consider the influence of the parameter ε. In Figure 4 the results of the FCM obtained with p = 20 for two different values of ε, i.e. ε = 10 -01 and ε = 10 -14 are presented. Small deviations from the exact solution are observed in the case of ε = 10 -01 . These deviations are due to the fact that a value different from ε = 0 inside the hole corresponds to a modification of the original problem, replacing the hole by a (very) soft material. Therefore, it can not be expected that a p-extension of the FCM converges to the exact solution of the problem when ε = 0 inside the hole. In order to study the influence of ε let us consider the convergence of a p-extension for different values of ε inside the hole, see Figure 5. It can be observed that the error e 2 H 1 converges exponentially down to a certain threshold which depends on the value of ε. The smaller we choose ε the higher is the achievable accuracy. The influence of ε is investigated more systematically in Figure 6, where the error e 2 plotted against ε obtained with three cells with a polynomial degree of p = 20. From this it can be observed that a p-extension on the mesh with three cells, with nodes being aligned to the hole, yields exponential convergence up to the error ε 2 in quadratic energy. This is in coherence with Theorem 4.2. Again, we consider a mesh with three cells but this time non-matching with the hole. The nodal coordinates X c of the mesh correspond to {-1, -0.3, 0.3, 1}. Therefore the nodes of the middle element are slightly outside of the hole, or in other words the hole is located completely inside the middle element. A comparison of the exact solution with the FCM approximation with p = 20 and ε = 10 -14 is presented in Figure 7. The convergence of a p-extension applying the FCM for different values of ε exact solution approximation -1 +1

exact solution approximation -1 +1 x u(x) exact solution approximation -1 +1 x u(x)
x u(x) shows again an exponential convergence up to a certain threshold depending on the chosen value of ε. However, in this example where the grid is not matching with the hole, the convergence is not as fast as in the case of the matching grid. The convergence of the error e 2 error e 2 H 1 of a p-extension applying the FCM for different values of ε is plotted in Figure 11. First of all, we observe a similar behavior as in the previous example, i.e. an exponential convergence can be obtained which is limited by the value of ε. However, increasing the polynomial degree further on can result also in an increase of the error. This effect can be explained by the poor conditioning of the resulting equation system observed by the increase of the number of iterations of the preconditioned conjugate gradient method which is applied to solve the overall equation system. The poor conditioning is due to the fact that one cell is completely inside the hole and therefore almost no stiffness is related to the corresponding degrees of freedom of that element. Increasing the polynomial degree further on deteriorates the situation and therefore round-off errors start to accumulate. Considering the scale of the y-axis of Figure 11 reveals that still very accurate results can be obtained with the FCM. Figure 12 presents, as in the previous examples, the dependency of the error e 2 H 1 on ε. Although the conditioning problem interferes in this investigation, the numerical results are again in good coherence with Theorem 4.1, stating that the error in quadratic energy converges linearly in ε.

Influence of integration error

Next, we consider computations based on a mesh with two equidistant cells, i.e. with nodal coordinates X c corresponding to {-1, 0, 1}. A comparison of the exact solution with the FCM approximation with p = 20 and ε = 10 -14 is given in Figure 13. Note, that an exact integration has been carried out by applying a composite Gaussian integration. Since the convergence of the error e 2 H 1 with respect to exact solution approximation -1 +1

x u(x) a p-extension and the choice of ε is very similar to the results of the non-matching grid with the hole being completely located inside the middle element, see Section 5.4.3, they are not presented here in detail. Summarizing the results, a p-extension on the (non-matching) mesh with two cells yields an exponential convergence up to the error ε in quadratic energy e 2 H 1 . Since in two and three spatial dimensions an exact integration of the stiffness and mass matrix is in general impossible, we investigate also the influence of the quality of the quadrature. The integration of the matrices of the two equidistant cells can be carried out exactly when applying n sc = 4 uniform sub-cells with n G = p + 1 Gaussian points on each sub-cell. Here we choose n sc such that we cannot perform an exact integration in order to investigate the influence of the quadrature. We choose ε = 0 in order to exclude a modelling error and focus only on the influence of the integration. In Figure 14 the convergence of the error e 2 H 1 is plotted as a function of the number of sub-cells n sc in a double logarithmic style. From the figure it is evident that the convergence is up to n -2 sc in quadratic energy. Here the size L sc of the sub-cells is uniform. We notice that the error e H 1 is proportional to L sc . In fact, the integration error is concentrated in the two sub-cells which cross Γ. In practice, a geometrical refinement towards Γ is used, so that the size of the sub-cells crossing Γ is very small. In [START_REF] Abedian | The finite cell method for the J 2 flow theory of plasticity[END_REF][START_REF]Performance of different integration schemes in facing discontinuites in the finite cell method[END_REF] an adaptive integration algorithm based on quadtree/octree schemes is proposed. To this end, the leaves of the quadtree/octree are introduced as the sub-cells on which the composite integration of the underlying cell is performed. The refinement process of the adaptive integration can be controlled by applying different criteria. One possibility is to refine the quadtree/octree by geometrical considerations, i.e. only those sub-cells are refined which are intersected by an interface or boundary. Alternatively the convergence during the adaptive integration of selected quantities such as the area or volume or selected entries of the stiffness matrix can be considered to control the refinement of the sub-cells.

A 2D benchmark of linear elasticity

Finally, we study a two-dimensional benchmark problem of linear elasticity which was defined in [START_REF] Stein | Error-Controlled Adaptive Finite Elements in Solid Mechanics[END_REF] to compare different adaptive finite element strategies. The benchmark problem is a two-dimensional plate under plane strain condition. Due to symmetry it is sufficient to discretize one quarter of the system, see Figure 15. The width and height are b = h = 100 mm and the radius is r = 10 mm. Linear isotropic elasticity with Young's modulus E = 206900 MPa and Poisson's ratio of ν = 0.29 is assumed. The plate is loaded by a traction of p = 450 MPa. The quantities to be computed are given in Table 1 in which also the reference values are listed. The plate is discretized with 2 × 2 quadrilateral cells on 1: Reference values which the tensor product space utilizing hierarchic shape functions [START_REF] Szab | Finite element analysis[END_REF][START_REF] Üster | The p-version of the finite element method for three-dimensional curved thin walled structures[END_REF] is used to discretize the trial and test functions. In Figure 16 the FCM grid as well as the sub-cells that are introduced for integration purposes only are presented. The sub-cells are generated in a fully automatic way by means of a space partitioning scheme based on a quadtree. We set ε = 10 -14 inside the hole independent of the degree p and the quadtree is refined towards the boundary of the circle. Note that in practical computations ε is usually set to 10 -6 , yielding a sufficiently small model error for engineering problems. The larger ε poses, due to a small condition number of the resulting linear equation system, less demands on the solver. The leafs of the quadtree correspond to the sub-cells that are used for the adaptive quadrature [START_REF]Performance of different integration schemes in facing discontinuites in the finite cell method[END_REF]. On each of the sub-cells a Gaussian quadrature is performed to accurately compute the stiffness matrix of the cell that is cut by the circle. As a first result the relative error in energy norm is plotted in Figure 17 demonstrating clearly the exponential convergence of the p-extension. Similar to the one-dimensional test example the optimal type of convergence for the FCM performing a pextension can be obtained also in the case of holes in two-dimensions. In order not to hinder the optimal convergence rate, the integration of the cell cut by the circle has been carried out with a very high accuracy, demonstrating that the exponential convergence can be observed within a p-extension from p = 1 up to p = 20. In practice, however, a finer grid of cells with a moderate polynomial degree of p = 6, ..., 8 would lead to accuracies which are relevant for engineering decisions. The FCM results for

Conclusions

A mathematical analysis of the finite cell method has been presented, proving exponential rate of convergence which has been observed earlier by numerical investigations, only. The sufficient smoothness conditions that we require on the exact solution are similar to those that ensure the exponential convergence for the classical p or hp-version of the finite element method. Furthermore, the dependence of an inherent modelling error on the computational scheme's penalization parameter was proved and numerically confirmed. This high order fictitious domain approach thus not only yields significant advantages over FEM concerning engineering applications by virtually eliminating the necessity to generate a finite element mesh, it also guarantees convergence properties, which are unachievable for low order methods.

This paper concentrates on Neumann boundary conditions at the transition from the physical to the fictitious domain. Numerical experiments using Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] to apply Dirichlet conditions show [START_REF] Schillinger | Small and large deformation analysis with the p-and B-spline versions of the finite cell method[END_REF][START_REF] Ruess | Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method[END_REF] that also in this case exponential convergence rates can be obtained. A mathematical investigation confirming this observation still has to be done in future work.

Figure 1 :

 1 Figure 1: Domain Ω and fictitious domain D The aim is to solve an elliptic equation Lu = f on Ω with mixed boundary conditions, namely Dirichlet conditions on Γ and Neumann conditions on Σ by solving Dirichlet problems on D. So D plays the role of a fictitious domain.

Remark 4 . 4

 44 If corners (and edges in dimension n = 3) are present on the external part Σ of the boundary, a local hp-extension could be used near Σ, leaving unchanged the FCM around the holes.

Figure 2 :

 2 Figure 2: Comparison of exact solution and FCM approximation with one cell and p = 20; odd problem (upper part), even problem (lower part)

  odd problem, p = 1, 2, 3, ... even problem, p = 1, 2, 3, ...

Figure 3 :

 3 Figure 3: p-Extension on mesh with one element, ε = 0, exact integration

Figure 4 :

 4 Figure 4: Comparison of exact solution and FCM approximation with p = 20 and ε = 10 -01 (upper part) and ε = 10 -14 (lower part)

ε 2 H 1

 21 = 10 -01 ε = 10 -02 ε = 10 -04 ε = 10 -06 ε = 10 -08 ε = 10 -10 ε = 10 -12 ε = 10 -14 error e polynomial degree p

Figure 5 :

 5 Figure 5: p-Extension on a mesh with three cells for different values of ε

Figure 6 :

 6 Figure 6: Influence of ε on the error

Figure 7 :

 7 Figure 7: Comparison of exact solution and FCM approximation with p = 20 and ε = 10 -14

H 1 with respect to ε is plotted in

Figure 8 :Figure 9 .Figure 9 :Figure 10 :

 89910 Figure 8: p-Extension on mesh with three elements for different ε values Figure 9. From this it is evident that the error in quadratic energy depends linearly on ε. This numerical

ε 2 H 1

 21 = 10 -01 ε = 10 -02 ε = 10 -04 ε = 10 -06 ε = 10 -08 ε = 10 -10 ε = 10 -12 ε = 10 -14 error e polynomial degree p

Figure 11 :Figure 12 :

 1112 Figure 11: p-Extension on mesh with three elements for different ε values

Figure 13 :

 13 Figure 13: Comparison of exact solution and FCM approximation with p = 20 and ε = 10 -14

Figure 14 :

 14 Figure 14: Convergence with respect to the number of sub-cells n sc . On each sub-cell n G = p + 1 Gaussian points are used.

Figure 15 : 5 U

 155 Figure 15: Square plate with a circular hole taking advantage of symmetry

Figure 16 :

 16 Figure 16: FCM grid with 2 × 2 cells (left) and corresponding sub-cells for integration purposes (right)

2 ×

 2 2 cells, p = 1, 2, 3, ...

Figure 17 : 2 × 2 ×

 1722 Figure17: Convergence of the error in energy norm for a p-extension on 2 × 2 cells the displacement u y at point 4 and the displacement u x at point 5 are depicted in Figure18showing a fast convergence of the p-extension also for point-wise quantities like displacements. More challenging than the results at points 4 and 5 are those quantities which are computed directly at the boundary of the hole. Therefore we study also the convergence of the displacement u x and stress component σ yy at point 2, see Figure19. Again, a fast convergence towards the reference values can be observed even for the stress component σ yy .

Figure 18 : 2 ×

 182 Figure 18: Convergence of the displacement u y at point 4 and u x at point 5 for a p-extension on 2 × 2 cells

2 ×

 2 2 cells, p = 1, 2, 3, ...

Figure 19 :

 19 Figure 19: Convergence of the displacement u x and stress component σ yy at point 2 for a p-extension on 2 × 2 cells

  Finally, we deduce (3.3) by dividing (3.10) by |w ap | a , taking the sup in w ap ∈ V ap , and coming back to(3.5).The second estimate (3.4) is obtained by dividing (3.9) by |w ap | a ap and using (3.2) next.

2

ap RHS of (3.9) .

(3.10)

  .3) Since ū0 ∈ A(Ω ap ), the fundamental approximation result of the p-version[START_REF] Guo | The h-p version of the finite element method for elliptic equations of order 2m[END_REF][START_REF] Babu | The p and h-p versions of the finite element method, basic principles and properties[END_REF][START_REF]Approximation properties of the h-p version of the finite element method[END_REF][START_REF] Schwab | p-and hp-finite element methods, theory and applications in solid and fluid mechanics[END_REF] provides exponential convergence. In dimension n = 1, this is a consequence of the analysis performed in[START_REF] Schwab | p-and hp-finite element methods, theory and applications in solid and fluid mechanics[END_REF] Section 3.3]. For dimension n = 2 (and also n = 1), it is easier to rely on [13, Appendix A], where various interpolants are constructed, providing error estimates of exponential type for analytic functions in tensor H 1 -norms 2 . The procedure to pass from n = 1 to n = 2 relies on tensorial properties of the discrete and continuous spaces. A similar procedure provides corresponding results in dimension n = 3. These interpolants reproduce nodal values. We obtain C 0 conformity using tensor trace liftings from edges (and faces in dimension n = 3), see[START_REF] Costabel | Exponential convergence of hp-FEM for Maxwell equations with weighted regularization in polygonal domains[END_REF] Prop. A.2]. Thus we can find v p

The operator of extension by 0 from V ′ (Ω) into H -1 (D) has to be understood as the dual of the restriction operator from H 1 0 (D) onto V (Ω).

1D test problem with Neumann boundary conditions at the hole 5.1 Problem definition and exact solution

As the simplest possible model for a hole in 1D, we choose the two component domain

Thus, the "hole" H is the interval (-1 4 , 1 4 ). The associated fictitious domain D is D = (-1, 1).

(5.2)

We consider the family of bilinear forms, indexed by the coefficient ε

Here, the coefficient ε inside the hole has been set to zero only for the mass matrix and small for the stiffness matrix. It corresponds to our general setting with

The domain Ω is symmetric with respect to the origin. We investigate two problems with different symmetry properties: the first one is odd, and the second, even.

The odd problem is described by

with the corresponding exact solution

e k(0.5-x) + e kx e -0.5k + e k if x > 0, -e k(0.5+x) + e -kx e -0.5k + e k if x < 0.

(5.5)

The even problem reads

where

e k(0.5-x) + e kx e -0.5k + e k if x > 0 e k(0.5+x) + e -kx e -0.5k + e k if x < 0

(5.7)

denotes the exact solution.