CutFEM without cutting the mesh cells: a new way to impose Dirichlet and Neumann boundary conditions on unfitted meshes

TL;DR

This paper introduces a novel CutFEM approach that avoids integrating over cut elements, maintaining optimal convergence rates for Poisson problems with Dirichlet or Neumann boundary conditions on unfitted meshes.

Contribution

It proposes a new method that eliminates the need for integration over cut elements in CutFEM, simplifying implementation while preserving accuracy.

Findings

Achieves optimal convergence rates theoretically.

Numerical results confirm the method's effectiveness.

Simplifies CutFEM implementation without loss of accuracy.

Abstract

We present a method of CutFEM type for the Poisson problem with either Dirichlet or Neumann boundary conditions. The computational mesh is obtained from a background (typically uniform Cartesian) mesh by retaining only the elements intersecting the domain where the problem is posed. The resulting mesh does not thus fit the boundary of the problem domain. Several finite element methods (XFEM, CutFEM) adapted to such meshes have been recently proposed. The originality of the present article consists in avoiding integration over the elements cut by the boundary of the problem domain, while preserving the optimal convergence rates, as confirmed by both the theoretical estimates and the numerical results.

Equations

$$-\Delta u=f\text{ in }\Omega ,u=g\text{ on }\Gamma$$

$$-\Delta u=f\text{ in }\Omega ,\frac{\partial u}{\partial n}=g\text{ on }\Gamma$$

$${\mathcal{T}}_h=\{T\in {\mathcal{T}}_h^{\mathcal{O}}:T\cap \Omega \ne \varnothing \},{\Omega }_h=({\cup }_{T\in {\mathcal{T}}_h}T{)}^{\circ }$$

$$-\Delta u=f\text{ in }{\Omega }_h,u=g\text{ on }\Gamma .$$

$${\int }_{{\Omega }_h}\nabla u\cdot \nabla v-{\int }_{{\Gamma }_h}\frac{\partial u}{\partial n}v+{\int }_{\Gamma }u\frac{\partial v}{\partial n}+\frac{\gamma }{h}{\int }_{\Gamma }uv={\int }_{{\Omega }_h}fv+{\int }_{\Gamma }g\frac{\partial v}{\partial n}+\frac{\gamma }{h}{\int }_{\Gamma }gv,\forall v\in H^1({\Omega }_h)$$

$$-\Delta u=f\text{ in }{\Omega }_h,\frac{\partial u}{\partial n}=g\text{ on }\Gamma .$$

$${\int }_{{\Omega }_h}\nabla u\cdot \nabla v-{\int }_{{\Gamma }_h}\frac{\partial u}{\partial n}v+{\int }_{\Gamma }\frac{\partial u}{\partial n}v={\int }_{{\Omega }_h}fv+{\int }_{\Gamma }gv$$

$${\int }_{{\Omega }_h}\nabla u\cdot \nabla v$$

$$={\int }_{{\Omega }_h}fv+{\int }_{\Gamma }gv,$$

$$V_h=\{v_h\in H^1({\Omega }_h):v_h{|}_T\in {\mathbb{P}}_1(T)\forall T\in {\mathcal{T}}_h\}$$

$$a_h(u_h,v_h)={\int }_{{\Omega }_h}f{v}_h+{\int }_{\Gamma }g\frac{\partial v_h}{\partial n}+\frac{\gamma }{h}{\int }_{\Gamma }g{v}_h\forall v_h\in V_h$$

$$a_h(u,v)={\int }_{{\Omega }_h}\nabla u\cdot \nabla v-{\int }_{{\Gamma }_h}\frac{\partial u}{\partial n}v+{\int }_{\Gamma }u\frac{\partial v}{\partial n}+\frac{\gamma }{h}{\int }_{\Gamma }uv+\sigma h\sum _{E\in {\mathcal{F}}_{\Gamma }}{\int }_E\left[\frac{\partial u}{\partial n}\right]\left[\frac{\partial v}{\partial n}\right]$$

$${\mathcal{F}}_{\Gamma }=\{E\text{ (an internal facet of }{\mathcal{T}}_h)\text{ such that }\exists T\in {\mathcal{T}}_h:T\cap \Gamma \ne \varnothing \text{ and }E\in \partial T\}$$

$${\tilde{V}}_h=\left\{v_h\in V_h:{\int }_{{\Omega }_h}=0\right\}$$

$$Z_h=\{z_h\in H^1({\Omega }_h^{\Gamma }{)}^d:z_h{|}_T\in {\mathbb{P}}_1(T{)}^d\forall T\in {\mathcal{T}}_h^{\Gamma }\}$$

$${\mathcal{T}}_h^{\Gamma }=\{T\in {\mathcal{T}}_h:T\cap \Gamma =\varnothing \}\text{ and }{\Omega }_h^{\Gamma }=({\cup }_{T\in {\mathcal{T}}_h^{\Gamma }}T{)}^{\circ }$$

$$a_h^N(u_h,y_h;v_h,z_h)={\int }_{{\Omega }_h}f{v}_h+{\int }_{\Gamma }g{v}_h+{\gamma }_{div}{\int }_{{\Omega }_h^{\Gamma }}f\mathrm{div}{z}_h\forall (v_h,z_h)\in {\tilde{V}}_h\times Z_h$$

$$a_h^N(u,y;v,z)={\int }_{{\Omega }_h}\nabla u\cdot \nabla v$$

$$+{\gamma }_1{\int }_{{\Omega }_h^{\Gamma }}(y+\nabla u)\cdot (z+\nabla v)+\sigma h{\int }_{{\Gamma }_h^i}\left[\frac{\partial u}{\partial n}\right]\left[\frac{\partial v}{\partial n}\right]$$

$$\Vert v{\Vert }_{0,{\Omega }_h^{\Gamma }}\le C\left(\sqrt{h}\Vert v{\Vert }_{0,\Gamma }+h|v{|}_{1,{\Omega }_h^{\Gamma }}\right)$$

$$\Vert v{\Vert }_{0,{\Omega }_h^{\Gamma }({\xi }_1)\cap {\mathcal{O}}_k}^2$$

$$⩽C{\int }_0^L{\int }_{-a({\xi }_1)}^{b({\xi }_1)}\left(v^2({\xi }_1,0)+{\left({\int }_0^{{\xi }_2}\frac{\partial v}{\partial {\xi }_2}({\xi }_1,t)dt\right)}^2\right)d{\xi }_{2}d{\xi }_1$$

$$⩽C{\int }_0^L\left(h{v}^2({\xi }_1,0)+h^2{\int }_{-a({\xi }_1)}^{b({\xi }_1)}{\left|\frac{\partial v}{\partial {\xi }_2}({\xi }_1,t)\right|}^{2}dt\right)d{\xi }_1$$

$$\le Ch\Vert v{\Vert }_{0,\Gamma \cap {\mathcal{O}}_k}^2+C{h}^2{\Vert \nabla v\Vert }_{0,{\Omega }_h^{\Gamma }\cap {\mathcal{O}}_k}^2$$

$$\Vert v{\Vert }_{0,{\Omega }_h^{\Gamma }}\le C\left(h^{1/2}\Vert v{\Vert }_{0,\Gamma }+h^{3/2}\Vert \nabla v{\Vert }_{0,\Gamma }+h^2\Vert v{\Vert }_{2,{\Omega }_h^{\Gamma }}\right)$$

$$|v{|}_{1,{\Omega }_h^{\Gamma }}\le C\left(\sqrt{h}\Vert \nabla v{\Vert }_{0,\Gamma }+h|v{|}_{2,{\Omega }_h^{\Gamma }}\right)$$

$$\Vert v{\Vert }_{0,{\Gamma }_h}\le C\left(\Vert v{\Vert }_{0,\Gamma }+h\Vert \nabla v{\Vert }_{0,\Gamma }+h^{3/2}\Vert v{\Vert }_{2,{\Omega }_h^{\Gamma }}\right)$$

$$\Vert v{\Vert }_{0,{\Gamma }_h\cap {\mathcal{O}}_k}^2$$

$$⩽C{\int }_0^L\left(v^2({\xi }_1,0)+h^2\frac{\partial v}{\partial {\xi }_2}({\xi }_1,0{)}^2+h^3\left|{\int }_0^{b({\xi }_1)}\frac{{\partial }^{2}v}{\partial {\xi }_2^2}({\xi }_1,t{)}^{2}dt\right|\right)d{\xi }_1$$

$$⩽C(\Vert v{\Vert }_{0,\Gamma \cap {\mathcal{O}}_k}^2+h^2\Vert \nabla v{\Vert }_{0,\Gamma \cap {\mathcal{O}}_k}^2+h^3\Vert v{\Vert }_{2,{\Omega }_h^{\Gamma }\cap {\mathcal{O}}_k}^2)$$

$$\sum _{E\in {\mathcal{F}}_{\Gamma }}\Vert v{\Vert }_{0,E}^2\le C(\Vert v{\Vert }_{0,\Gamma }^2+h|v{|}_{1,{\Omega }_h^{\Gamma }}^2).$$

$$\Vert v{\Vert }_{0,{\Gamma }_h}^2\le C(\Vert v{\Vert }_{0,\Gamma }^2+h|v{|}_{1,{\Omega }_h^{\Gamma }}^2).$$

$$\Vert v{\Vert }_{0,E}^2\le C\left(\frac{1}{h}\Vert v{\Vert }_{0,T}^2+h|v{|}_{1,T}^2\right)$$

$$\sum _{E\in {\mathcal{F}}_{\Gamma }}\Vert v_h{\Vert }_{0,E}^2\le C\left(\frac{1}{h}\Vert v{\Vert }_{0,{\Omega }_h^{\Gamma }}^2+h|v{|}_{1,{\Omega }_h^{\Gamma }}^2\right)$$

Full text

1 Introduction.

In this article, we propose a new approach to the numerical solution of boundary value problems for partial differential equations using finite elements on non-matching meshes, circumventing the need to generate a mesh accurately fitting the physical boundaries or interfaces. Such approaches, classically known as the fictitious domain methods, have a long history dating back to [1] in the case of the Poisson problem with Dirichlet boundary conditions. They were later popularized, by Glowinski and co-workers, cf. [2] for example, and successfully applied in the context of particular flow simulations [3, 4]. The basic idea of the classical fictitious domain method is to embed the physical domain $\Omega$ into a bigger simply shaped domain $\mathcal{O}$, to extend the physically meaningful solution on $\Omega$ by a fictitious solution on $\mathcal{O}\setminus\Omega$ using the same governing equations as on $\Omega$, and to impose the boundary conditions on $\partial\Omega$ by Lagrange multipliers. At the numerical level, this means that one can work with simple meshes on $\mathcal{O}$, but one also needs a mesh on the physical boundary $\partial\Omega$ for the Lagrange multiplier, which should be coarser than the first one in order to satisfy the inf-sup condition [5]. One is thus not completely free from meshing problems. Another unfortunate feature of such methods is their poor accuracy: one cannot expect the convergence order to be better than $1/2$ with the respect to the mesh size. Closely related penalty methods are well suited for both Dirichlet and Neumann boundary conditions [6, 7], may be simpler to implement in practice than the fictitious domain methods with Lagrange multipliers, but share with them the poor convergence properties.

More recently, several optimally convergent finite element methods on non-matching meshes were proposed following the XFEM or CutFEM paradigms. XFEM (extended finite element method) was initially introduced in [8] for applications in the structural mechanics on cracked domains (Neumann boundary conditions on the crack). Its ability to impose Dirichlet boundary conditions was demonstrated in [9, 10] and a properly stabilized version with proved optimal convergence was proposed in [11]. The CutFEM methods [12] were first introduced in a series of papers by Burman and Hansbo [13, 14, 15]. The common feature of all these methods is that the simple background mesh is only used to define the finite element space (only the mesh elements having non empty intersection with $\Omega$ are kept in the computational mesh), but the solution is no longer extended from the physical domain $\Omega$ to a larger fictitious domain. The integrals over $\Omega$ are thus maintained in the finite element formulation. The boundary conditions are imposed either through Lagrange multipliers living on the same mesh as the primary solution [11, 13] or by the Nitsche method [14, 15] stabilized by the ghost penalty [16]. The optimal convergence, i.e. the error estimates of the same order as those for finite element methods on a comparable matching mesh, are established for all the methods above.

As already mentioned, XFEM/CutFEM methods contain the integrals over $\Omega$ in their formulations. This can be cumbersome in practice. Citing [13] “the only remaining difficulty of implementation is the actual integration on the boundary and on parts of elements cut by the boundary. This difficulty however is expected to arise in any optimal order fictitious domain method.” We attempt in the present article to prove that this statement is not entirely true. We propose in fact an optimal order fictitious domain method that does not involve the integrals on $\Omega$ and thus does not require to perform the integration on parts of elements cut by the boundary. Our method can be regarded as a mix between the classical fictitious domain approach and CutFEM. As in CutFEM, we use the the computational mesh constructed by keeping only the elements from the background mesh on the embedding domain $\mathcal{O}$ having non empty intersection with $\Omega$. On the other hand, we extend the solution from $\Omega$ to the computational domain, which is now only slightly larger than $\Omega$. In fact, the extension is done on a narrow band of width of order of the mesh size, contrary to the extension to entire $\mathcal{O}$ as in the classical fictitious domain approach.111The idea of constructing numerically a smooth extension to the whole $\mathcal{O}$ is explored in [17] resulting in an optimally convergent method. The price to pay is the necessity to solve an optimization problem by an iterative process, which can be expensive in practice. This minimizes the effect of choosing a “wrong” extension and enables us, with the help of a proper stabilization mainly borrowed from CutFEM, to preserve the optimal convergence without integration on the cut mesh elements. Unfortunately, the integration on the actual boundary should still be performed, so that some non trivial quadrature techniques are still needed. We believe however that an integration over the surface (resp. the curve) cut by the mesh is simpler to implement than that on the cut mesh elements in 3D (resp. 2D).

Let us now give a first sketch of the methods that we propose in this article. We restrict ourselves to the simple model problem: the Poisson equation with either Dirichlet boundary conditions

[TABLE]

or Neumann ones

[TABLE]

Here $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ is a domain with smooth boundary $\Gamma=\partial\Omega$, $f$ and $g$ are given functions on $\Omega$ and $\Gamma$ respectively (satisfying the compatibility condition $\int_{\Omega}f+\int_{\Gamma}g=0$ in the Neumann case).

We start by embedding $\Omega$ into a simply shaped domain $\mathcal{O}$ and introduce a quasi-uniform mesh $\mathcal{T}_{h}^{\mathcal{O}}$ on $\mathcal{O}$ consisting of triangles/tetrahedrons of maximum diameter $h$ that can be cut by the boundary $\Gamma$ in an arbitrary manner. The computational mesh $\mathcal{T}_{h}$ is obtained from $\mathcal{T}_{h}^{\mathcal{O}}$ by dropping out all the mesh element lying completely outside $\Omega$:

[TABLE]

as illustrated in Fig. 1. $\Omega_{h}$ is thus the domain occupied by the computational mesh, slightly larger than $\Omega$. We denote its boundary by $\Gamma_{h}:=\partial\Omega_{h}$.

We outline first the derivation of our method in the case of Dirichlet boundary conditions, following the preliminary version in [18]. We assume that the right-hand side $f$ is extended from $\Omega$ to $\Omega_{h}$ and imagine (for the moment) that (1) can be solved on the extended domain $\Omega_{h}$ while still imposing the boundary conditions on $\Gamma$:

[TABLE]

We keep here the same notations $u$ and $f$ for the functions on $\Omega_{h}$ as for the originals on $\Omega$. Integration by parts over $\Omega_{h}$ and imposing the boundary conditions weakly on $\Gamma$ as in the antisymmetric variant [19] of the Nitsche method [20] yields

[TABLE]

with $\gamma>0$. Here, $n$ on $\Gamma$ (resp. $\Gamma_{h}$) denotes the unit normal looking outwards from $\Omega$ (resp. $\Omega_{h}$). Our finite element method is then based on the weak formulation (5) adding to it a ghost penalty stabilization to assure the coerciveness of the bilinear form. The method will be fully presented in Subsection 2.1 and analyzed in Subsections 3.2 and 3.3. Unlike the preliminary version in [18], the optimal error estimates are here proved under the natural assumptions on the regularity of the data: $f\in L^{2}(\Omega_{h})$, $g\in H^{3/2}(\Gamma)$ (vs. $f\in H^{1}(\Omega_{h})$ in [18]). We stress again that formulation (5) does not contain integrals on $\Omega$. Otherwise, the resulting method is very close to the antisymmetric Nitsche CutFEM method from [16].

We turn now to the case of Neumann boundary conditions (2). We start again by extending $f$ from $\Omega$ to $\Omega_{h}$ and imagine (for the moment) that (2) can be solved on the extended domain $\Omega_{h}$ while still imposing the boundary conditions on $\Gamma$:

[TABLE]

Integration by parts over $\Omega_{h}$ and imposing the boundary conditions weakly on $\Gamma$ would yield

[TABLE]

Such a formulation does not seem to lead to a reasonable finite element method. It is difficult to imagine a stabilization that would make the bilinear form on the left-hand side above coercive; unlike the Dirichlet case, one cannot rely on the smallness of $u$ near $\Gamma$ to control the non-positive terms in the bilinear form. Note that a Nitsche method for the Neumann boundary conditions is proposed in [21] in the mesh conforming case. However, the terms added there to the natural variational formulation tend to weaken its coercivenes, rather than to enhance it.

Fortunately, we are able to devise an optimally convergent method by introducing a reconstruction $y=-\nabla u$ of the gradient of $u$ on the band of cut mesh elements $\Omega_{h}^{\Gamma}$, cf. Fig. 1 on the right, and by replacing the weak formulation above with the following one

[TABLE]

with $\gamma_{1}>0$. The terms multiplied by $\gamma_{1}$ serve to impose $y=-\nabla u$ on $\Omega_{h}^{\Gamma}$. Our finite element method (13), based on the weak formulation (7) with the addition of grad-div and ghost stabilizations, will be fully presented in Subsection 2.2 and analyzed in Subsections 3.4 and 3.5, under the assumption of some extra regularity on the data: $f\in H^{1}(\Omega_{h})$, $g\in H^{3/2}(\Gamma)$ so that $u\in H^{3}(\Omega)$. It involves the additional vector-valued unknown $y_{h}$, which results in some extra cost in practice, as compared with a simpler method for Dirichlet boundary conditions. Fortunately, this extra cost is negligible as $h\to 0$ since $y_{h}$ lives only on the mesh elements cut by $\Gamma$ which constitute smaller and smaller proportion of the total number of elements as the meshes are refined.

Our finite element methods for both Dirichlet and Neumann cases are presented in more detail in the next section. Their well-posedness and the optimal error estimates are proved in Section 3. We restrict ourselves here to $P_{1}$ continuous finite elements on a triangular/tetrahedral mesh. An extension to higher-order elements $P_{k}$ and a possible adaptation to Robin boundary conditions are sketched in Section 4. We present some numerical experiments in Section 5, restricting ourselves to the $P_{1}$ elements. Some final conclusions are provided in section 6.

2 Presentation of the methods

2.1 Dirichlet boundary conditions

We present first our discretization of Problem (1). Recall that $\mathcal{T}_{h}$ is a quasi-uniform mesh obtained from a larger background mesh retaining only the elements lying inside $\Omega$ or cut by $\Gamma$, cf. (3) and Fig. 1, and $\Omega_{h}$ is the corresoinding domain with boundary $\Gamma_{h}$. We inspire ourselves from the variational formulation (5), introduce the finite element space

[TABLE]

with $\mathbb{P}_{1}$ denoting the set of polynomials of degree $\leq 1$, and introduce the following discrete problem:

Find $u_{h}\in V_{h}$ such that

[TABLE]

where

[TABLE]

and $\gamma,\sigma$ are some positive numbers properly chosen in a manner independent of $h$. The last term in (10) is the ghost penalty [16]. It is crucial to assure the coerciveness of $a_{h}$. The notations here are as follows: $[\cdot]$ stands for the jump over an internal facet of mesh $\mathcal{T}_{h}$ and

[TABLE]

The ghost penalty is thus a properly scaled sum of the jumps of the normal derivatives over all the internal facets of the mesh either cut by $\Gamma$ themselves or owned by a mesh element cut by $\Gamma$.

2.2 Neumann boundary conditions

We turn now to Problem (2). We inspire ourselves with the variational formulation (7) and introduce a subspace of $V_{h}$ (continuous $P_{1}$ finite elements on mesh $\mathcal{T}_{h}$ as defined above) appropriate to the treatment of Neumann boundary conditions

[TABLE]

and an auxiliary finite element space

[TABLE]

Here $\mathcal{T}_{h}^{\Gamma}$ represents the cut elements of the mesh and $\Omega_{h}^{\Gamma}$ the corresponding subdomain of $\Omega_{h}$:

[TABLE]

cf. Fig. 1, right. We shall also denote by $\Gamma_{h}^{i}$ the internal boundary of $\Omega_{h}^{\Gamma}$, i.e. the ensemble of the facets separating $\Omega_{h}^{\Gamma}$ from the mesh elements inside $\Omega$, so that $\partial\Omega_{h}^{\Gamma}=\Gamma_{h}\cup\Gamma_{h}^{i}$.

Our finite element problem is: Find $u_{h}\in\tilde{V}_{h}$, $y_{h}\in Z_{h}$ such that

[TABLE]

where

[TABLE]

with some positive numbers $\gamma_{div}$, $\gamma_{1}$ and $\sigma$ properly chosen in a manner independent of $h$. In addition to the variational formulation (7), we have introduced here

a grad-div stabilization (the terms multiplied by $\gamma_{div}$) in the vein of [22], which is consistent with the governing equations since $y=-\nabla u$ and thus $\operatorname{div}y=-\Delta u=f$ on $\Omega_{h}^{\Gamma}$;

• 2.

a reduced ghost stabilization (the term multiplied by $\sigma h$ in (14). Contrary to the Dirichlet case (10), it is sufficient to penalize the normal derivative jumps only on the facets separating the cut elements from the uncut ones.

Remark 2.1**.**

The introduction of the subspace $\tilde{V}_{h}$ is helpful to ensure the well-posedness of problem 14 (otherwise, if posed on $V_{h}\times Z_{h}$, the functions $u=1$, $y=0$ would be in the kernel of $a^{N}_{h}$ reflecting the fact that the exact solution is not unique without an additional constraint). We have chosen to impose this constraint on $\Omega_{h}$ rather than on $\Omega$ in line with our desire to avoid the integrals on $\Omega$, difficult to implement in practice.

3 The theoretical error analysis.

3.1 Geometrical assumptions and technical lemmas.

The theoretical analysis of the methods presented above will be done under the following assumptions on the boundary $\Gamma$ and on the subdomain $\Omega_{h}^{\Gamma}$ covered by the cut mesh elements, as defined in (12). Both these assumptions are typically satisfied if the mesh is sufficiently refined with respect to $\Gamma$.

Assumption 1**.**

The subdomain $\Omega_{h}^{\Gamma}$ can be covered by open sets $\mathcal{O}_{k}$, $k=1,\ldots,K$ and one can introduce on every $\mathcal{O}_{k}$ local coordinates $\xi_{1},\xi_{2}$ in the 2D case, $d=2$ (resp. $\xi_{1},\xi_{2},\xi_{3}$ if $d=3$) such that

$\Gamma\cap\mathcal{O}_{k}$* is given by $\xi_{d}=0$ and $\Omega\cap\mathcal{O}_{k}$ by $\xi_{d}<0$;*

• 2.

$\Omega_{h}^{\Gamma}\cap\mathcal{O}_{k}$* is given by $-a(\xi_{1},\xi_{2})<\xi_{d}<b(\xi_{1},\xi_{2})$ with some continuous non-negative functions $a$ and $b$;*

• 3.

$b(\xi_{1},\xi_{2})+a(\xi_{1},\xi_{2})\leq C_{1}h$* with some $C_{1}>0$;*

• 4.

all the partial derivatives $\partial\xi_{i}/\partial x_{j}$ and $\partial x_{i}/\partial\xi_{j}$ are bounded by some $C_{2}>0$.

Moreover, each point of $\Omega_{h}^{\Gamma}$ is covered by at most $N_{int}$ sets $\mathcal{O}_{k}$.

Assumption 2**.**

The boundary $\Gamma$ can be covered by element patches $\{\Pi_{k}\}_{k=1,\ldots,N_{\Pi}}$ having the following properties:

Each patch $\Pi_{k}$ is a connected set;

• 2.

Each $\Pi_{k}$ is composed of a mesh element $T_{k}$ lying inside $\Omega$ and some mesh elements cut by $\Gamma$, more precisely $\Pi_{k}=T_{k}\cup\Pi_{k}^{\Gamma}$ where $T_{k}\in\mathcal{T}_{h}$, $T_{k}\subset\bar{\Omega}$, $\Pi_{k}^{\Gamma}\subset\mathcal{T}_{h}^{\Gamma}$, and $\Pi_{k}^{\Gamma}$ contains at most $M$ mesh elements;

• 3.

$\mathcal{T}_{h}^{\Gamma}=\cup_{k=1}^{N_{\Pi}}\Pi_{k}^{\Gamma}$;

• 4.

$\Pi_{k}$* and $\Pi_{l}$ are disjoint if $k\neq l$.*

In what follows, we suppose both assumptions above to hold true and use the notation $C$ for positive constants (which can change from one instance to another) that depend only on $C_{1},C_{2},N_{int},M$ from the assumptions above and on the mesh regularity. We also recall that mesh $\mathcal{T}_{h}$ is supposed quasi-uniform.

Assumption 1 is quite standard and essentially tells us that the boundary $\Gamma$ is smooth and not too wiggly on the length scale $h$, so that one can be sure that the band of cut mesh elements is of width $\sim h$. Assumption 2 is slightly more technical. Let us explain the construction of element patches evoked there, cf. Fig. 2. Each patch $\Pi_{k}$ is assigned to a facet $E_{k}\subset\Gamma_{h}^{i}$ separating the cut elements from the interior ones. To form the patch $\Pi_{k}$, one takes a mesh element $T_{k}$ lying inside $\Omega$ and attached to $E_{k}$. One then picks up several cut elements touching $E_{k}$ to form $\Pi_{k}^{\Gamma}$ and set $\Pi_{k}=T_{k}\cup\Pi_{k}^{\Gamma}$. If, again, the boundary $\Gamma$ is not too wiggly with respect to the mesh $\mathcal{T}_{h}$, one can partition the cut elements between the patches $\Pi_{k}$ so that each patch contains a small number of elements (typically from 2 to 4 in 2D, and slightly more in 3D).

We begin with two technical lemmas. The first one is in the vein of Poincaré inequalities taking into account the assumption that $\Omega_{h}^{\Gamma}$ is a strip of width $\sim h$ around $\Gamma$.

Lemma 3.1**.**

For all $v\in H^{1}(\Omega_{h}^{\Gamma})$,

[TABLE]

Proof.

Consider the 2D case ($d=2$). Recalling Assumption 1, we can pass to the local coordinates $\xi_{1},\xi_{2}$ on every set $\mathcal{O}_{k}$ assuming that $\xi_{1}$ varies between [math] and $L$, and to use the bounds on $\xi_{2}$ and on the mapping $(x_{1},x_{2})\mapsto(\xi_{1},\xi_{2})$ to write

[TABLE]

Summing over all neighborhoods $\mathcal{O}_{k}$ gives (15). The proof in the 3D case is the same up to the change of notations. ∎

Corollary 3.2**.**

For all $v\in H^{2}(\Omega_{h}^{\Gamma})$,

[TABLE]

Proof.

Applying Lemma 3.1 to $\nabla v$ we get

[TABLE]

Substituting this into (15) leads to (16). ∎

Lemma 3.3**.**

For all $v\in H^{2}(\Omega_{h}^{\Gamma})$,

[TABLE]

Proof.

As in the proof of Lemma 3.1, we only consider the 2D case. Using again Assumption 1, we can represent $v|_{\Gamma_{h}}$ inside every $\mathcal{O}_{k}$ by a Taylor expansion with the integral remainder

[TABLE]

Summing this over all $\mathcal{O}_{k}$ gives (18). ∎

Lemma 3.4**.**

For all $v\in H^{1}(\Omega_{h}^{\Gamma})$,

[TABLE]

and

[TABLE]

Proof.

Let $E$ be a mesh facet belonging to a mesh element $T\in\mathcal{T}_{h}$. Recall the well-known trace inequality

[TABLE]

Summing this over all $E\in\mathcal{F}_{\Gamma}$ gives

[TABLE]

leading, in combination with (15), to (19). The proof of (20) is similar: it suffices to take the sum over the facets composing $\Gamma_{h}$. ∎

3.2 Dirichlet boundary conditions: coerciveness of $a_{h}$.

We turn to the study of method (9)–(10) for the problem with Dirichlet boundary conditions. Our first goal is to prove the coerciveness of the bilinear form $a_{h}$ uniformly in $h$. The proof of this result, cf. Lemma 3.6, will be based on the following

Lemma 3.5**.**

For any $\beta>0$ one can choose $0<\alpha<1$ depending only on the mesh regularity and on parameter $M$ from Assumption 2 such that

[TABLE]

for all $v_{h}\in V_{h}$.

Proof.

Choose any $\beta>0$, consider the decomposition of $\Omega_{h}^{\Gamma}$ in element patches $\{\Pi_{k}\}$ as in Assumption 2, cf. also Fig. 2. Introduce

[TABLE]

where the maximum is taken over all the possible configurations of a patch $\Pi_{k}$ allowed by the mesh regularity and over all the piecewise linear functions on $\Pi_{k}$. The subset $\mathcal{F}_{k}\subset\mathcal{F}_{\Gamma}$ gathers the facets internal to $\Pi_{k}$. Note that the quantity under the $\max$ sign in (3.2) is invariant under the scaling transformation $x\mapsto hx$ and is homogeneous with respect to $v_{h}$. Recall also that the patch $\Pi_{k}$ contains at most $M$ elements. Thus, the maximum is indeed attained since it is taken over a bounded set in a finite dimensional space (one can restrict the maximization to patches of diameter of order 1, since the quantity to maximize is invariant under rescaling).

Clearly, $\alpha\leq 1$. Supposing $\alpha=1$ would lead to a contradiction. Indeed, if $\alpha=1$ then we can take $\Pi_{k}$, $v_{h}$ yielding this maximum and suppose without loss of generality $|v_{h}|_{1,\Pi_{k}}=1$. We observe then

[TABLE]

since $|v_{h}|_{1,\Pi_{k}}^{2}=|v_{h}|_{1,T_{k}}^{2}+|v_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}$. This implies $\nabla v_{h}=0$ on $T_{k}$ and $[\frac{\partial v_{h}}{\partial n}]=0$ on all $E\in\mathcal{F}_{k}$, thus $\nabla v_{h}=0$ on $\Pi_{k}$, which contradicts $|v_{h}|_{1,\Pi_{k}}=1$.

This proves $\alpha<1$. We have thus

[TABLE]

for all $v_{h}\in V_{h}$ and all the admissible patches $\Pi_{k}$. Summing this over $\Pi_{k}$, $k=1,\ldots,N_{\Pi}$ yields (21). ∎

Lemma 3.6**.**

Provided $\sigma$ is sufficiently big, there exists an $h$-independent constant $c>0$ such that

[TABLE]

for all $v_{h}\in V_{h}$.

Proof.

Let $v_{h}\in V_{h}$ and $B_{h}$ be the strip between $\Gamma$ and $\Gamma_{h}$, i.e. $B_{h}=\Omega_{h}\setminus\Omega$ and $\partial B_{h}=\Gamma\cup\Gamma_{h}$. Recall that we assume the normal $n$ to look outward from $\Omega_{h}$ (resp. $\Omega$) on $\Gamma_{h}$ (resp. $\Gamma$). The outward looking normal on $\partial B_{h}$ coincides with $n$ on $\Gamma_{h}$ and with $-n$ on $\Gamma$. We have thus

[TABLE]

Substituting this into the definition (10) of $a_{h}$ yields

[TABLE]

Noting that $B_{h}\subset\Omega_{h}^{\Gamma}$ we can use (21) combined with the Young inequality (for any $\varepsilon>0$) and (19) to write

[TABLE]

Taking $\varepsilon$ sufficiently small and $\sigma$ sufficiently big this bounds $a(v_{h},v_{h})$ from below by $c{\left|\kern-0.51663pt\left|\kern-0.51663pt\left|v_{h}\right|\kern-0.51663pt\right|\kern-0.51663pt\right|}_{h}^{2}$ as claimed. ∎

3.3 Dirichlet boundary conditions: the error estimates.

We can now establish an optimal $H^{1}$ error estimate for the method (9)–(10) using the coerciveness of $a_{h}$ provided by Lemma 3.6. An $L^{2}$ error estimate will follow in Theorem 3.8.

Theorem 3.7**.**

Suppose $f\in L^{2}(\Omega_{h})$, $g\in H^{3/2}(\Gamma)$ and let $u\in H^{2}(\Omega)$ be the solution to (1), $u_{h}\in V_{h}$ be the solution to (9)–(10). Provided $\sigma$ is sufficiently big, there exists an $h$-independent constant $C>0$ such that

[TABLE]

Proof.

Under the Theorem’s assumptions, the solution to (1) is indeed in $H^{2}(\Omega)$ and it can be extended to a function $\tilde{u}\in H^{2}(\Omega_{h})$ such that $\tilde{u}=u$ on $\Omega$ and $\|\tilde{u}\|_{2,\Omega_{h}}\leq C(\|f\|_{0,\Omega}+\|g\|_{3/2,\Gamma})$, cf. [23]. Clearly, $\tilde{u}$ satisfies

[TABLE]

with $\tilde{f}:=-\Delta\tilde{u}$. It entails a Galerkin orthogonality relation

[TABLE]

We have then using the standard nodal interpolation $I_{h}:C(\bar{\Omega}_{h})\to V_{h}$ and recalling Lemma 3.6

[TABLE]

with $e_{u}=\tilde{u}-I_{h}\tilde{u}$. Using the definition (10) of $a_{h}$, we can now bound term by term

[TABLE]

We have used here Lemma 3.4 to bound $h\sum_{E\in\mathcal{F}_{\Gamma}}\left\|\left[\frac{\partial v_{h}}{\partial n}\right]\right\|_{0,E}^{2}$ and $\|v_{h}\|_{0,\Gamma_{h}}$. This entails thanks to the usual interpolation estimates

[TABLE]

Moreover, $f=\tilde{f}$ on $\Omega$ so that

[TABLE]

thanks to Lemma 3.1 (recall that $\Omega_{h}\setminus\Omega\subset\Omega_{h}^{\Gamma}$). We conclude

[TABLE]

Combining the estimates above with the triangle inequality proves ${\left|\kern-0.51663pt\left|\kern-0.51663pt\left|u_{h}-\tilde{u}\right|\kern-0.51663pt\right|\kern-0.51663pt\right|}_{h}\leq Ch(\|f\|_{0,\Omega_{h}}+\|g\|_{3/2,\Gamma})$, as claimed. ∎

Remark 3.1**.**

The proof above does not rely on a solution to the non-standard boundary value problem (4) on $\Omega_{h}$. We rather use the well defined solution $u$ to problem (1) and extend it to $\Omega_{h}$.

The following theorem gives an $L^{2}$ estimate for the error. It is sub-optimal, although the numerical experiments reveal the optimal convergence rate $O(h^{2})$, similar to the state of the art in the study of the non-symmetric Nitsche method.

Theorem 3.8**.**

Under the assumptions of Theorem 3.7, there exists an $h$-independent constant $C>0$ such that

[TABLE]

Proof.

Let us introduce $w:\Omega\to\mathbb{R}$ such that

[TABLE]

By elliptic regularity, $\|w\|_{2,\Omega}\leq C\|u-u_{h}\|_{0,\Omega}$. Let $\tilde{w}$ be an extension of $w$ from $\Omega$ to $\Omega_{h}$ preserving the $H^{2}$ norm estimate. Applying Corollary 3.2 to $\tilde{w}$ and Lemma 3.1 to $\nabla\tilde{w}$ yields

[TABLE]

Similarly, applying Lemma 3.3 to $\tilde{w}$ yields

[TABLE]

Taking $w_{h}=I_{h}\tilde{w}$, we can summarize the bounds above together with the interpolation estimates as

[TABLE]

Using Galerkin orthogonality (24) and estimates (27) we arrive at (recall $B_{h}=\Omega_{h}\setminus\Omega$)

[TABLE]

which gives the announced error estimate in $L^{2}(\Omega)$ norm thanks to the error estimate in the triple norm ${\left|\kern-0.51663pt\left|\kern-0.51663pt\left|\cdot\right|\kern-0.51663pt\right|\kern-0.51663pt\right|}$ already established in the proof of Theorem 3.7. We also recall $\|\tilde{f}\|_{0,\Omega_{h}^{\Gamma}}=\|\Delta\tilde{u}\|_{0,\Omega_{h}}\leq C\left(\|f\|_{0,\Omega_{h}}+\|g\|_{3/2,\Gamma}\right)$. ∎

3.4 Neumann boundary conditions: coerciveness of $a^{N}_{h}$.

We turn to the study of method (13)–(14) for the problem with Neumann boundary conditions. Our first goal is to prove the coerciveness of the bilinear form $a^{N}_{h}$ uniformly in $h$. To this end, we note that $a^{N}_{h}$ can be rewritten using the divergence Theorem as

[TABLE]

provided $y$ and $v$ are of regularity $H^{1}$. We have denoted here again $B_{h}=\Omega_{h}\setminus\Omega$. The following lemma will allow us to control $\int_{B_{h}}y\cdot\nabla v$ while the other term $\int_{B_{h}}v\operatorname{div}y$ can be controlled thanks to the grad-div stabilization.

Lemma 3.9**.**

For any $\beta>0$, there exist $0<\alpha<1$ and $\delta>0$ depending only on the mesh regularity and on parameter $M$ in Assumption 2 such that

[TABLE]

for all $v_{h}\in V_{h},z_{h}\in Z_{h}$.

Proof.

The boundary $\Gamma$ can be covered by element patches $\{\Pi_{k}\}_{k=1,\ldots,N_{\Pi}}$ as in Assumption 2. Choose any $\beta>0$ and consider

[TABLE]

where the maximum is taken over all the possible configurations of a patch $\Pi_{k}$ allowed by the mesh regularity and over all the piecewise linear functions $v_{h}$ and $z_{h}$ on $\Pi_{k}$. Note that the quantity under the $\max$ sign in (30) is invariant under the scaling transformation $x\mapsto hx$, $z_{h}\mapsto\frac{1}{h}z_{h}$, $v_{h}\mapsto v_{h}$ and is homogeneous with respect to $v_{h}$, $z_{h}$. Thus, the maximum is indeed attained since it is taken over a bounded set in a finite dimensional space.

Clearly, $\alpha\leqslant 1$. Supposing $\alpha=1$ would lead to a contradiction. Indeed, if $\alpha=1$, we can then take $\Pi_{k}$, $v_{h}$, $z_{h}$ yielding this maximum (in particular, $|v_{h}|_{1,\Pi_{k}}^{2}+\|z_{h}\|_{0,\Pi_{k}^{\Gamma}}^{2}>0$). We observe then

[TABLE]

and consequently (recall $|v_{h}|_{1,\Pi_{k}}^{2}=|v_{h}|_{1,T_{k}}^{2}+|v_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}$)

[TABLE]

This implies $\|z_{h}+\nabla v_{h}\|_{0,\Pi_{k}^{\Gamma}}$=0, i.e. $\nabla v_{h}=-z_{h}$ on $\Pi_{k}^{\Gamma}$. Since $\nabla v_{h}$ is piecewise constant and $z_{h}$ is continuous, it means that $\nabla v_{h}=-z_{h}=$const on $\Pi_{k}^{\Gamma}$. We also have $\left\|\left[\frac{\partial v_{h}}{\partial n}\right]\right\|_{0,\partial T_{k}\cap\partial\Pi_{k}^{\Gamma}}=0$, which implies $\nabla v_{h}=$const on the whole $\Pi_{k}$. Since, by (31), $\nabla v_{h}=0$ on $T_{k}$, we have finally $\nabla v_{h}=0$ on $\Pi_{k}$ and $z_{h}=0$ on $\Pi_{k}^{\Gamma}$, which is in contradiction with $|v_{h}|_{1,\Pi_{k}}^{2}+\|z_{h}\|_{0,\Pi_{k}^{\Gamma}}^{2}>0$.

Thus $\alpha<1$ and

[TABLE]

for all $v_{h},z_{h}$ and all admissible patches $\Pi_{k}$. We now observe

[TABLE]

for any $\varepsilon>0$. We have used here the following estimate valid by Young inequality

[TABLE]

[TABLE]

Taking $\varepsilon$ sufficiently small, redefining $\alpha$ as $\alpha\left(1+\frac{\varepsilon}{2}\right)$ and putting $\delta=\left(\beta+\frac{\alpha}{2}+\frac{\alpha}{2\varepsilon}\right)$ we obtain (29). ∎

Lemma 3.10**.**

Provided $\gamma_{div},\gamma_{1}$ are sufficiently big, there exists an $h$-independent constant $c>0$ such that

[TABLE]

with

[TABLE]

Proof.

Expression (28) for $a^{N}_{h}$ implies for all $v_{h}\in\tilde{V}_{h}$, $z_{h}\in Z_{h}$

[TABLE]

A combination of (29) with the Young inequality (for any $\varepsilon>0$) yields

[TABLE]

We have for any $v_{h}\in\tilde{V}_{h}$

[TABLE]

with a mesh independent constant $C_{p}$. This is the standard Poincaré inequality (recall that $\int_{\Omega_{h}}v_{h}=0$) but one should check that $C_{p}$ does not depend on $\Omega_{h}$. This can be easily done thanks to the Poincaré inequality on $\Omega$. Indeed, denoting $\bar{v}_{h,\Omega}=\frac{1}{|\Omega|}\int_{\Omega}v_{h}$ we have

[TABLE]

and

[TABLE]

The ratio of domain measures $|\Omega_{h}\setminus\Omega|/|\Omega|$ is uniformly bounded under the assumption of $h$ sufficiently small. We have thus, by Poincaré inequality on $\Omega$ (with a constant which is obviously mesh independent) and then by Lemma 3.1

[TABLE]

Finally, invoking the trace inequality $\|v_{h}\|_{0,\partial\Omega}\leqslant C(\|v_{h}\|_{0,\Omega}+|v_{h}|_{1,\Omega})$ and recalling that $h$ is sufficiently small, we obtain (33).

Substituting (33) into (32) entails

[TABLE]

Taking $\varepsilon,\beta$ sufficiently small and $\gamma_{1},\gamma_{div}$ sufficiently big this bounds $a^{N}_{h}(v_{h},z_{h};v_{h},z_{h})$ from below by $c{\left|\kern-0.51663pt\left|\kern-0.51663pt\left|v_{h},z_{h}\right|\kern-0.51663pt\right|\kern-0.51663pt\right|}_{h}^{2}$ as claimed. ∎

3.5 Neumann boundary conditions: error estimates.

We can now establish an optimal $H^{1}$ error estimate for the method (13)–(14) using the coerciveness of $a^{N}_{h}$ provided by Lemma 3.10. An $L^{2}$ error estimate will follow in Theorem 3.12.

Theorem 3.11**.**

Suppose $f\in H^{1}(\Omega_{h})$, $g\in H^{3/2}(\Gamma)$ and let $u\in H^{3}(\Omega)$ be a solution to (2), $(u_{h},y_{h})\in\tilde{V}_{h}\times Z_{h}$ be the solution to (13). Provided $\gamma_{div},\gamma_{1}$ are sufficiently big, there exists an $h$-independent constant $C>0$ such that

[TABLE]

Proof.

Under the Theorem’s assumptions, the solution to (2) is indeed in $H^{3}(\Omega)$ and it can be extended to a function $\tilde{u}\in H^{3}(\Omega_{h})$ such that $\tilde{u}=u$ on $\Omega$ and $\|\tilde{u}\|_{3,\Omega_{h}}\leq C(\|f\|_{1,\Omega}+\|g\|_{3/2,\Gamma})$. Introduce $y=-\nabla\tilde{u}$ on $\Omega_{h}^{\Gamma}$. Clearly, $\tilde{u}$, $y$ satisfy

[TABLE]

with $\tilde{f}:=-\Delta\tilde{u}$. It entails a Galerkin orthogonality relation

[TABLE]

We have then using the standard nodal interpolation $I_{h}:C(\bar{\Omega}_{h})\to\tilde{V}_{h}$ and recalling Lemma 3.10

[TABLE]

with $e_{u}=\tilde{u}-I_{h}\tilde{u},e_{y}=y-I_{h}y$. Recalling (28), we can bound

[TABLE]

By the usual interpolation estimates this entails

[TABLE]

since $y=-\nabla\tilde{u}$ on $\Omega_{h}^{\Gamma}$. Moreover,

[TABLE]

thanks to (33). We recall now that $f=\tilde{f}$ on $\Omega$ so that, thanks to Lemma 3.1,

[TABLE]

and conclude

[TABLE]

Combining the estimates above with the triangle inequality proves ${\left|\kern-0.51663pt\left|\kern-0.51663pt\left|u_{h}-\tilde{u},y_{h}-y\right|\kern-0.51663pt\right|\kern-0.51663pt\right|}_{h}\leq Ch(\|f\|_{1,\Omega_{h}}+\|g\|_{5/2,\Gamma})$. ∎

Remark 3.2**.**

The proof above does not rely on a solution to the non-standard boundary value problem (6) on $\Omega_{h}$. We rather use the well defined solution $u$ to problem (2) and extend it to $\Omega_{h}$. The optimal convergence is then obtained at the expense of a stronger than usual assumption on the right-hand side in (1): we need $u\in H^{3},f\in H^{1}$ where as $u\in H^{2},f\in L^{2}$ suffices for standard $P_{1}$ finite elements on a conforming mesh.

Theorem 3.12**.**

Under the assumptions of Theorem 3.11, there exists an $h$-independent constant $C>0$ such that

[TABLE]

Proof.

Let us introduce $w:\Omega\to\mathbb{R}$ such that

[TABLE]

with $c\in\mathbb{R}$ chosen to ensure that this problem is well posed, i.e. $c=\frac{1}{|\Omega|}\int_{\Omega}(u-u_{h})$. By elliptic regularity, $\|w\|_{2,\Omega}\leq C\|u-u_{h}-c\|_{0,\Omega}$. Let $\tilde{w}$ be an extension of $w$ from $\Omega$ to $\Omega_{h}$ preserving the $H^{2}$ norm estimate and set $w_{h}=I_{h}\tilde{w}$. Integration by parts and interpolation estimates yield

[TABLE]

We now rewrite the bilinear form (28) as

[TABLE]

so that the Galerkin orthogonality relation (35) with $v_{h}=w_{h}-\bar{w}_{h,\Omega_{h}},z_{h}=0$ and $\bar{w}_{h,\Omega_{h}}$ the average of $w_{h}$ over $\Omega_{h}$ becomes

[TABLE]

Recalling the definition of the triple norm from Lemma 3.10, this leads to

[TABLE]

By Lemma 3.1 and interpolation estimates

[TABLE]

and similarly

[TABLE]

Finally,

[TABLE]

Combining all the estimates above we arrive at

[TABLE]

We conclude recalling the bound on the triple norm of the error from the proof of Theorem 3.11, (36), and the regularity estimate $\|\tilde{w}\|_{2,\Omega_{h}}\leq C\|u-u_{h}-c\|_{0,\Omega}$. ∎

4 Extension to $P_{k}$ finite elements and Robin boundary conditions

The methods presented above can be easily extended to $P_{k}$ finite elements giving optimal convergence of order $h^{k}$ in the $H^{1}$-norm, $k\geq 2$. We thus consider in this section the finite element space

[TABLE]

with $\mathbb{P}_{k}(T)$ representing the polynomials of degree $\leq k$, propose some modifications to be introduced to the methods above, first for Dirichlet boundary conditions and then for Neumann-Neumann ones, and outline the convergence proofs in both these cases.

4.1 Dirichlet boundary conditions

The method (9)–(10) is modified as follows:

Find $u_{h}\in V_{h}^{(k)}$ s.t.

[TABLE]

with

[TABLE]

Note that the additional stabilization term with the product $(\Delta u)(\Delta v)$ is strongly consistent since $-\Delta u=f$. Also note that the ghost penalty term is extended to control the normal derivatives of all orders up to $k$, cf. [15].

Revisiting the theoretical analysis of Section 3 reveals the following:

Lemma 3.5 remains valid thanks to the extended ghost penalty. Indeed, the only thing to recheck in its proof is the following implication: if $\nabla v_{h}=0$ on $T_{k}$ and all the norms of the jumps contained in ghost penalty vanish, then $\nabla v_{h}=0$ on $\Pi_{k}^{\Gamma}$. This is true since the extended ghost penalty controls all the derivatives present in our finite element space.

• 2.

The first relation (22) in the proof of Lemma 3.6 now becomes

[TABLE]

since one can no longer assume $\Delta v_{h}=0$ on $T$. The new term with $\Delta v_{h}$ can be then controlled thanks to the additional $\Delta\cdot\Delta$ stabilization term in the bilinear form $a_{h}$ so that Lemma 3.6 remains valid.

• 3.

In Theorem 3.7, we should now suppose $f\in H^{k-1}(\Omega_{h})$, $g\in H^{k+1/2}(\Gamma)$ so that the solution $u$ to (1) is in $H^{k+1}(\Omega)$. The proof of the Theorem can be then followed using appropriate Sobolev spaces and standard interpolation estimates to $P_{k}$ finite elements. The only non-trivial change in the proof concerns the estimate for $\|f-\tilde{f}\|_{0,\Omega_{h}^{\Gamma}}$, which should be changed to

[TABLE]

This can be proved by refining the argument of Lemma 3.1. We recall that $f=\tilde{f}$ on $\Omega$ and both $f$ and $\tilde{f}$are in $H^{k-1}(\Omega)$ so that all the derivatives up to order $k-2$ of $f$ and $\tilde{f}$ coincide on $\Gamma$. The bound (39) then follows as in Lemma 3.1 employing a Taylor expansion with the integral remainder of order $k-1$.

Keeping in mind the modifications above, we can easily establish the convergence of the $P_{k}$ version of our method: given $f\in H^{k-1}(\Omega_{h})$ and $g\in H^{k+1/2}(\Gamma)$, and supposing $\sigma$ sufficiently big, one has

[TABLE]

The $L^{2}$ error estimate of order $h^{k+1/2}$ can be also proved as in Theorem 3.8.

4.2 Neumann boundary conditions

We turn now to Problem (2). The goal is to extend the method (13)–(14) to $P_{k}$ finite elements. We thus recall the space $V_{h}^{(k)}$ as in (38), introduce its subspace $\tilde{V}_{h}^{(k)}$ as in (11) and the auxiliary finite element space

[TABLE]

Our finite element problem is: Find $u_{h}\in\tilde{V}_{h}^{(k)}$, $y_{h}\in Z_{h}^{(k)}$ solving (13) with the bilinear form $a^{N}_{h}$ modified as follows

[TABLE]

We have added here extra penalization terms controlling the jumps of the higher normal derivatives of $u$ on $\Gamma_{h}^{i}$ and an additional ghost stabilization term controlling the jumps of derivatives of $y$ on the cut facets, denoted by $\mathcal{F}^{cut}_{\Gamma}$ ($\mathcal{F}^{cut}_{\Gamma}$ is thus the collection of interior facets of mesh $\mathcal{T}_{h}^{\Gamma}$; $\mathcal{F}^{cut}_{\Gamma}\subset\mathcal{F}_{\Gamma}$).

Keeping in mind the modifications above, we can easily establish the convergence of the $P_{k}$ version of our method: given $f\in H^{k}(\Omega_{h})$ and $g\in H^{k+1/2}(\Gamma)$, ans supposing $\gamma_{div},\gamma_{1}$ sufficiently big, one has

[TABLE]

The $L^{2}$ error estimate of order $h^{k+1/2}$ can be also proved as in Theorem 3.8. Note that we still need to assume some extra regularity: $f\in H^{k}(\Omega_{h})$ contrary to $f\in H^{k-1}(\Omega_{h})$ in the case of Dirichlet boundary conditions, or to what whould be necessary for the optimal convergence of the standard finite element method.

4.3 Robin boundary conditions

We can also consider the problem

[TABLE]

with $\kappa>0$. It is straightforward to adapt the method of the Neumann case (13)–(14) to the present Robin case. To this end, we rewrite the boundary conditions in (40) as $\frac{1}{\kappa}u+\frac{\partial u}{\partial n}=\frac{1}{\kappa}g$. By the same considerations as in the Neumann case we arrive then at the following method for (40) using $P_{k}$ finite elements (without the constraint of zero integral over $\Omega_{h}$): Find $u_{h}\in V_{h}^{(k)}$, $y_{h}\in Z_{h}^{(k)}$ such that

[TABLE]

The proofs above can be easily adapted to show the optimal convergence of this scheme both in $H^{1}(\Omega)$ and $L^{2}(\Omega)$ norms.

Alternatively, one can adapt the method of the Dirichlet case (9)–(10) mimicking the approach of [21]. In the context of unfitted meshes ($\Gamma\not=\Gamma_{h}$), choosing the antisymmetric variant of the Nitsche method and adding the ghost penalty, this would give the scheme222This approach has been suggested by an anonymous reviewer.: Find $u_{h}\in V_{h}^{(k)}$ such that

[TABLE]

The advantage of (42) over (41) is in the absence of the additional variable $y_{h}$. However, the analysis of this scheme is yet to be done and is out of the scope of the present paper.

5 Numerical experiments.

We shall illustrate our methods (9)–(10) and (13)–(14) by numerical experiments in a 2D domain $\Omega$ defined by a level-set function $\varphi$:

[TABLE]

where $(r,\theta)$ are the polar coordinates $r=\sqrt{x^{2}+y^{2}}$, $\theta=\arctan\frac{y}{x}$, and $R=0.47$ (this example is taken from [11]). To construct the computational mesh, we embed $\Omega$ into the square $\mathcal{O}=(-0.5,0.5)^{2}$, introduce a regular $N\times N$ criss-cross mesh on $\mathcal{O}$, and drop the triangles outside $\Omega$ to produce $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\Gamma}$, as illustrated at Fig. 1 in the case $N=16$. In some of our experiments, the domain $\Omega$ will be rotated by an angle $\theta_{0}$ counter-clockwise around the origin. This is achived by redefining $\theta$ as $\theta=\arctan\frac{y}{x}-\theta_{0}$.

All the computations are done in FreeFEM [24] taking advantage of its level-set capabilities (the key word levelset in numerical integration commands int1d and int2d to deal with the integrals over $\Gamma$ and over $\Omega$, cf. Subsection 5.2). Note that the numerical integration on $\Gamma$ introduces an additional error which is not covered by our theoretical analysis: $\Gamma$ is in fact approximated by a sequence of straight segments with the endpoints obtained by approximate intersections of $\Gamma$ with the edges of $\mathcal{T}_{h}$. This error is of order $h^{2}$ and thus can be assumed negligible in the case of $P_{1}$ finite elements (which is the only case studied numerically in this paper). We presume that a subtler approximation should be employed when dealing with higher order finite elements, as in [25].

5.1 Dirichlet boundary conditions

We start by solving numerically the Dirichlet problem (1) in domain (43) with zero right-hand side $f=0$ and a non-homogeneous boundary condition $g$ set so that the exact solution is given by

[TABLE]

We employ the method (9)–(10) taking the following parameter values: $\gamma=1$, $\sigma=0.01$. Fig. 3 represents the solution obtained on a $16\times 16$ mesh. We observe that the numerical method captures well the exact solution, the perturbation being essentially concentrated in the narrow fictitious domain $B_{h}=\Omega_{h}\setminus\Omega$. Fig. 4 reports the evolution of the error under the mesh refinement (always using the regular criss-cross meshes as mentioned above) and confirms the optimal convergence order of the method in both $H^{1}$ and $L^{2}$ norms.

In order to explore the robustness of the method with respect to the placement of the physical domain on the computational mesh, we now redo the calculations above, rotating $\Omega$ by a series of angles $\theta_{0}$ ranging from [math] to $\frac{2\pi}{7}$ as described in the preamble of this Section. For each rotation angle $\theta_{0}$, the boundary $\Gamma$ cuts the triangles of the background mesh in a different manner, creating sometimes the “dangerous” situations when certain mesh triangles of $\mathcal{T}_{h}$ have only a tiny portion inside the physical domain $\Omega$. Fig. 5 presents the errors in $L^{2}(\Omega)$ and $H^{1}(\Omega)$ norms as functions of $\theta_{0}$ at different discretization levels. We observe that the errors do not vary much from one position to another, especially when measured in the $H^{1}(\Omega)$ norm. Moreover, this variability decreases as the meshes are refined.

Finally, we explore the influence of the parameters $\gamma$ and $\sigma$ in (10) on the precision of the method. The error norms on a fixed $32\times 32$ mesh for various choices of $\gamma$ and $\sigma$ are presented at Fig. 6. We observe that the method is not very sensitive to the parameters, especially when the error is measured in the $H^{1}$ norm. The accuracy does not seem to deteriorate catastrophically even in the limit $\gamma,\sigma\to 0$. This is somewhat surprising in view of the theoretical analysis of Section 3 and may indicate that a subtler theory could reveal the optimal convergence properties of the method (9)–(10) without stabilization. In our numerical experiments above, we have preferred however to remain on a safer side and have chosen a rather large value of the Nitsche parameter $\gamma=1$ while keeping the ghost stabilization parameter $\sigma=10^{-2}$ small. Note that larger values of $\gamma$ seem to make the method more sensible and unpredictive with respect to the choice of $\sigma$. On the other hand, the error is clearly monotonically increasing with $\sigma$ in the regime $\gamma\leq 1$.

5.2 Comparisons with CutFEM

As already mentioned in Introduction, our method (9)–(10) is very close to CutFEM methods, the essential difference being that we avoid the integration over $\Omega$, i.e. the numerical integration over the cut mesh elements. Although such an integration is in principle difficult to implement (in particular, it is more difficult than the integration over the segments of the boundary $\Gamma$), it is already available in FreeFEM, so that we can easily compare numerically the performance of CutFEM with that of our method.

We have considered the following variants of CutFEM with $V_{h}$ denoting everywhere the $P_{1}$ finite elements space on the mesh $\mathcal{T}_{h}$, as in (8).

a)

A version with Lagrange multipliers approximated by $P_{0}$ finite elements on the cut triangles $\mathcal{T}_{h}^{\Gamma}$ [13]:

Find $u_{h}\in V_{h}$, $\lambda_{h}\in W_{h}:=\{\mu_{h}\in L^{2}(\Omega_{h}^{\Gamma}):\mu_{h}|_{T}\in\mathbb{P}_{0}(T)\ \forall T\in\mathcal{T}_{h}^{\Gamma}\}$ s.t.

[TABLE]

b)

A version based on the symmetric Nitsche method [14]:

Find $u_{h}\in V_{h}$ s.t.

[TABLE]

c)

A version based on the antisymmetric Nitsche method (similar to the above, only the sign in front of the third term changes) [16]:

Find $u_{h}\in V_{h}$ s.t.

[TABLE]

The results are presented at Fig. 7. We consider there the same setup as in Fig. 5, i.e. domain $\Omega$ given by (43) rotated around the origin at a series of angles $\theta_{0}$, the exact solution given by (44). The stabilization parameters $\gamma$ and $\sigma$ are given in the captions (we have taken the same parameters for the antisymmetric version as for our method in Fig. 5, while a larger value for $\gamma$ was necessary for the symmetric version). Comparing the results in Figs. 5 and 7, we observe that all the methods have overall almost the same accuracy in the $H^{1}(\Omega)$ norm, but CutFEM produces better results (gaining by a factor around 2) when the error is measured in the $L^{2}(\Omega)$ norm. However, the performance of CutFEM (with the exception of the antisymmetric Nitsche variant) can drastically degrade at certain position of the domain $\Omega$ with respect to the mesh (the spikes of the error on the graphs as functions of the rotation angle $\theta_{0}$). This is certainly an implementation issue, which can be conceivably explained by inaccuracies in the numerical integration over the cut triangles and/or by an incorrect determination of such triangles due to the round-off errors. We do not attempt here to investigate this issue further, and merely note that the absence of integration over the cut triangles in our method (9)–(10) permits us to avoid some delicate implementation issues.

5.3 Neumann boundary conditions

We now turn to the Poisson-Neumann problem (2). Our test case is similar to that used before: the domain is given by (43) and the exact solution by (44). The right-hand side in (2) is thus $f=0$ and the boundary condition $g$ is set up as $g=n\cdot\nabla u$ with the normal $n$ defined via the levelset function $\varphi$ in (43) as $n=\frac{\nabla\varphi}{|\nabla\varphi|}$. We employ the method (13)–(14) taking the following parameter values: $\gamma_{div}=1$, $\gamma_{1}=10$, $\sigma=0.01$. The constraint $\int_{\Omega}u_{h}=0$ is enforced with the help of a Lagrange multiplier, thus increasing the size of the system matrix by 1. A convergence study under the mesh refinement is reported at Fig. 8. It confirms the optimal convergence order of the method in both $H^{1}$ and $L^{2}$ norms.

In order to explore the robustness of the method with respect to the placement of the physical domain on the computational mesh, we now redo the calculations above, rotating $\Omega$ by a series of angles $\theta_{0}$ ranging from [math] to $\frac{2\pi}{7}$, same as in the Dirichlet case. Fig. 9 presents the errors in $L^{2}(\Omega)$ and $H^{1}(\Omega)$ norms as functions of $\theta_{0}$ at different discretization levels. Overall, the errors are of the same order as in the Dirichlet case, but they are now much more sensitive to positioning the domain with respect to the mesh. This is especially true for the $L^{2}$ error whose variability does not fade out when the meshes are refined.

The origin of this unfortunate phenomenon can be partially attributed to possible bad conditioning of the matrix of method (13)–(14). Indeed, the scaling of the grad-div stabilization term there is inconsistent with other terms. Loosely speaking, assuming that $y$ scales as $\frac{1}{h}u$, the scaling of the term with $\gamma_{div}$ with respect to the mesh size is $1/h^{2}$ in 2D, while the scaling of all the other terms is 1. We have tried to recover the correct scaling of this term in numerical experiments reported in Fig. 10. The coefficient $\gamma_{div}$ is taken there proportional to $h^{2}$ rather than a mesh-independent constant. This scheme is thus not covered by the preceding theory but it gives essentially the same results in practice as the scheme with constant $\gamma_{div}$. In particular, the lack of convergence in the $L^{2}$ norm on certain geometrical configurations is still observed on the most refined meshes.

It is also interesting to compare the results produced by the scheme (13)–(14) with those of CutFEM: Find $u_{h}\in\tilde{V}_{h}$ s.t.

[TABLE]

The results are reported in Fig. 11 taking the stabilization parameter $\sigma=0.01$. The error curves in the $H^{1}(\Omega)$ norm are practically the same as those for (13)–(14). We also observe once again a big variability of the $L^{2}(\Omega)$ with respect to the geometry. However, these variations seem to fade out under the mesh refinement, contrary to the scheme (13)–(14).

Extensive numerical experiments on the influence of the stabilization parameters $\gamma_{div}$, $\gamma_{1}$ and $\sigma$ in (14) have been also conducted. Similar to the Dirichlet case, the accuracy does not seem to deteriorate catastrophically in the limit $\gamma_{div},\sigma\to 0$ and the method is not too much sensitive to the parameters in a wide range. We choose not to go into further details since the behaviour with respect to 3 parameters is difficult to summarize efficiently in figures or tables.

5.4 A test case in 3D

In a slight variation to the preceding numerical experiments and to the theoretical analysis of the paper, we present here the results for the diffusion equation in a 3D domain:

[TABLE]

with a non-constant diffusion coefficient $D$. We take $\Omega\subset\mathbb{R}^{3}$ as the ball of radius 1 centred at the origin. Denoting $r$ the distance to the origin, we choose the diffusion coefficient and the exact solution as

[TABLE]

and adjust the right-hand side $f$ accordingly. The method (9)–(10) is modified as follows: Find $u_{h}\in V_{h}$ (the $P_{1}$ continuous finite elements on mesh $\mathcal{T}_{h}$) such that

[TABLE]

The numerical experiments have been done in FreeFEM using the HYPRE linear solvers via PETSc. The computational meshes $\mathcal{T}_{h}$ were obtained from the background uniform tetrahedral meshes on a sufficiently large cube by getting rid of the tetrahedra outside $\Omega$, cf. some example in Fig. 12. The following stabilization parameters were used: $\gamma=1$, $\sigma=0$. A mesh convergence study is reported in Fig. 13 and confirms the optimal convergence rates.

6 Conclusions

We have presented an optimally convergent method of the fictitious domain/XFEM/CutFEM type avoiding the numerical integration on cut mesh elements. The numerical experiments confirm the optimal convergence and the robustness of the method (with a possible slight deficiency on the level of $L^{2}$ errors in the case of Neumann boundary conditions). We have restricted ourselves to simple model problem (Poisson equation, easily generalisable to the diffusion equation) in this first publication, but the methods should be applicable in more realistic settings as, for example, elasticity problem on a cracked domain (the method would then be similar to XFEM of [8] but avoiding the integration on the mesh elements cut by the crack). Further research could be aimed, apart from extensions to more complex problems, at a finer understanding of the influence of stabilization parameters.

Acknowledgements

I am grateful to Grégoire Allaire for stimulating discussions and suggestions concerning the implementation. The implementation in FreFEM++ would not be possible without advices and help generously provided by Frédéric Hecht and Pierre Jolivet. I have also greatly appreciated the insightful and thought provoking comments of one of the anonymous reviewers.