A hierarchical NXFEM for fictitious domain simulations

We suggest a fictitious domain method, based on the Nitsche XFEM method of (Comput. Meth. Appl. Mech. Engrg 2002; 191:5537–5552), that employs a band of elements adjacent to the boundary. In contrast, the classical fictitious domain method uses Lagrange multipliers on a line (surface) where the boundary condition is to be enforced. The idea can be seen as an extension of the Chimera method of (ESAIM: Math. Model Numer. Anal. 2003; 37:495–514), but with a hierarchical representation of the discontinuous solution field. The hierarchical formulation is better suited for moving fictitious boundaries since the stiffness matrix of the underlying structured mesh can be retained during the computations.


INTRODUCTION
Fictitious domain methods were introduced in order to be able to use Cartesian meshes also for solving problems on domains with complex boundaries.The idea is to enforce Dirichlet boundary conditions on a given curve (surface) that is discretized independently of the mesh, cf.Glowinski et al. [5].On this curve, the boundary condition is enforced, typically by use of Lagrange multipliers.The system of equations can then be set up on a Cartesian mesh and the degrees of freedom falling outside of the boundary are discarded.The problem with this approach is that the derivatives of the finite element solution normal to the curve cannot accommodate the jump necessary to achieve optimal order convergence, cf.[4].Another problem is how to choose the discretization of the curve relative to the elements it crosses in order for the problem to be well posed.Guidelines for this are given in [4] but they are by necessity rather vague.
In this paper we introduce an alternative method based on the use of Nitsche's method in the vein of Hansbo et al. [7] (building on [2] and [6]), where overlapping meshes were considered.We shall also employ overlapping meshes in the form of (see Fig. 1) 1. the (structured) mesh on which the problem is set up and 2. a narrow band of elements that overlays the first mesh.
This allows for the direct use of the method proposed in [7], where the elements on the underlying mesh were cut by the overlying and the solution pasted together by use of Nitsche's method [9].The outer boundary of the band can then be used as the Dirichlet boundary.We remark at this point that another strength of the approach in [7] is that any boundary condition can be applied at the outer boundary of the band.This is not so straighforward in a classical fictitious domain method.
However, in [7] this was achieved by modifying the elements of the underlying mesh, which does not allow for the system matrix on the underlying mesh to remain unchanged.If the boundary were to move continuously, the entries in the system matrix would also have to be changed continuously.It is desirable to have a fixed matrix for the underlying problem and to see the imposition of the boundary condition as a set of additional degrees of freedom (as in the original fictitious domain approach where the additional degrees of freedom are the Lagrange multipliers).As has been noted by Areias and Belytschko [1], the method of [7] can be interpreted as an XFEM method (which instead adds degrees of freedom in a hierarchical fashion) by a reordering of the degrees of freedom.In this paper, we seize on the hierarchical concept to construct a fictitious domain method which is formally identical to that of [7] (thus benefitting from the optimal order error analysis therein) while still keeping the underlying system matrix unchanged by the location of the boundary.Our method shares some characteristics with the fat boundary method proposed by Maury [8], which also makes use of an auxiliary band-like domain.However, in [8], the original Poisson problem is recast on the continuous level and requires a fixed point iteration scheme to solve, like in classical domain decomposition methods.In contrast, our method is defined on the discrete level and can be solved monolithically.

The continuous problem
As a model problem, we consider Poisson's equation where Ω is a domain in R 2 with polygonal boundary ∂Ω and f is a given forcing term.We embed Ω in a larger rectangular domain Ω so that Ω is completely contained in the interior on Ω.Finally, we introduce a third domain Ω 1 consisting of a band (the width of which may be mesh dependent) whose outer boundary coincides with ∂Ω and whose inner boundary forms a line, denoted Γ, in the interior of Ω.The remainder of Ω is denoted by Clearly, the extension to three dimensions is straightforward.
We now rephrase the problem (1) as an interface problem.For u in Ω 1 ∪ Ω 2 we define the jump of u on Γ by [u] := u 1 | Γ − u 2 | Γ , where u i = u| Ωi is the restriction of u to Ω i .Conversely, for u i defined in Ω i we identify the pair {u 1 , u 2 } with the function u which equals u i on Ω i .Then we may formulate (1) as: Here n is the outward pointing unit normal to Ω 1 and For a bounded open connected domain D we shall use standard Sobolev spaces H r (D) with norm

Denote by T h
1 the triangulation of Ω 1 and by T h 2 the triangulation of Ω.We shall make a discretization on the whole of Ω even though the solution has no physical significance outside ∂Ω.This is in line with classical fictitious domain methods and means that the stiffness matrix assembled from T h 2 remains fixed even if the domain should change, as it must do in many applications.
We will use the following notation for mesh related quantities.Let h K be the diameter of an element K ∈ T h i and h = max K∈T h i ,i=1,2 h K .To distinguish elements from the two meshes, we will sometimes use indexed element notation K i ∈ T h i for clarity.Furthermore, we introduce G h as the set of elements in and the corresponding mesh-dependent domain We shall also need the mesh dependent boundary ∂Ω Gh , which consists of the edges on elements in G h that form the boundary of Ω Gh .This boundary is also split in the part exterior to Γ, ∂Ω ext Gh and interior to Γ, ∂Ω int Gh .

The nodes on Γ of the elements in T h
1 , together with the points of intersection between elements in T h 2 and Γ, define a partition of Γ, Γ = ∪ j∈Jh Γ j .Note that each part Γ j belongs to two elements, one from each mesh.We denote these elements by K j 1 and K j 2 , respectively.A local meshsize on Γ is defined by For any element K ∈ T h i , let P K = K ∩ Ω i denote the part of K in Ω i .We make the following assumptions regarding the meshes: A1: The triangulations are non-degenerate, i.e., where h K is the diameter of K and ρ K is the diameter of the largest ball contained in K. A2: The meshes have locally compatible meshsize over Γ.More precisely, let K j 1 ∈ T h 1 and K j 2 ∈ T h 2 be the elements which contain a specific part Γ j of Γ.We assume that Here and below, C and c denote generic constants.We shall seek a discrete solution where Note that functions in V h are, in general, discontinuous across Γ; the discontinuity is represented by the hierarchical space In Figure 3 we illustrate the concept: a discontinuous function on a one-dimensional element occupying the domain Ω Gh = (0, 1) (solid line) can be written as the sum of a continuous function (dashed line) from V h 2 and piecewise continuous functions which are zero in the nodes of the element (dash-dotted line) from V h 3 and V h 4 .Note that even though functions in V h 3 and V h 4 are defined on all of Ω Gh , we shall only use those parts that live on the respective side of Γ, corresponding to the situation in Fig. 3

The Finite Element Method
The method is defined by the variational problem: find U ∈ V h such that where, if we denote , with f extended, e.g., by zero outside Ω, and where h is the local meshsize (3).Here, the jump [U] is interpreted as U 1 −(U 2 +U 3 ).The continuity conditions of u and ∂ n u at Γ are satisfied weakly by means of a variant of Nitsche's method [9] for consistent weak enforcement of Dirichlet boundary conditions.
To ensure stability, the parameter λ has to be taken sufficiently large, cf.[7].
To analyze the method and show its equivalence with that of [7] we introduce a second bilinear form and right-hand side It is straightforward to show that the method analyzed in [7] would then read: find where , since we may ignore the solution outside Ω we set φ 2 | Ω\(Ω2∪ΩG h ) = 0. Note that the integrals here are taken only over the domain Ω, unlike the proposed method.However, since the solution is completely decoupled by the cut, it does not matter what we do outside of the domain (e.g., how we extend f and how the boundary conditions on ∂ Ω are specified).In the next section we prove that indeed U * = U| Ω and hence the analysis of [7] carries over to our formulation.

A PRIORI ERROR ESTIMATES
Consider the following mesh dependent norms: Let U denote the solution of (4) and U * the solution of (5).Then U| Ω = U * .PROOF.The proof proceeds in three steps 1. Show the existence of a subspace Ṽh ⊂ V h such that a h (U, φ) = a h * (U, φ) and l( φ) = l * ( φ), ∀ φ ∈ Ṽh .

Show that ∃ξ
3. Apply coercivity and Galerkin orthogonality to the discrete error using the results of [7].

Firstly let
Since the integrals on Γ are the same in a h (•, •) and a h * (•, •) we only need to prove the equivalence of the volume integrals.By using the decomposition Ω (6) where B Γ (U, φ) denotes the integrals over Γ.It then follows that for all φ ∈ Ṽh The equality l( φ) = l * ( φ) is shown in a similar fashion.Secondly observe that Ṽh only imposes constraints on components of the function outside Ω, since the constraint on φ 2 in Ω int Gh is compensated for by the freedom of φ 3 .Hence the existence of the sought ξ U and ξ U * .
Finally we recall the following coercivity result from [7], for some c > 0 and for all v ∈ V h there holds In particular this holds for v = U − U * and hence, since We have the following consistency relation.4) is consistent in the sense that, for u solving (2) there holds

Proposition 3.2. The discrete problem (
The proof is given in [7]. We can now directly take advantage of the theory developed in [7] which shows that we have optimal error estimates for any polynomial degree of the underlying finite element method: Theorem 3.3.For U solving (4) and u solving (2), the following a priori error estimates hold We refer to [7] for the proof.Here we shall only point out one of the crucial points in the analysis.Accuracy of the method is expressed by the orthogonality relation ( 9), but to show convergence we also need stability of the discrete problem as expressed by (8).In order to show that our system matrix is positive definite, we rely on the following inverse inequality, see, e.g., Warburton and Hesthaven [11].
Lemma 3.4.For φ ∈ V h , the following inverse inequality holds: The size of the constant C I can be found by solving a small local eigenvalue problem; explicit bounds are discussed in [11].The constant C I determines the size of λ, we are obliged to take λ > C 2 I p 2 in order to ensure coercivity.Here is the point of the band: if we just cut the mesh with ∂Ω and apply Nitsche's method on the cut elements, sliver elements would be generated that would require extremely large C I in Lemma 3.4, leading to severe ill conditioning of the discrete system.We illustrate the problem in the case p = 1.Thus, consider using applying Nithsche's method directly on an element cut by ∂Ω.Then C I can be found, for a given element K as the largest eigenvalue λ max in the eigenproblem of finding where P 1 K denotes the space of polynomials of degree 1 on K. Then ∇v is constant on each element and thus we have where meas(•) denotes the length, area, or volume of the object in question, and and it follows that and thus we must choose The situation is illustrated in Fig. 4 where ∂Ω is represented by the dashed line.We note that as ∂Ω moves to the left in Fig. 4, the area meas(K ∩ Ω) will approach zero while the length meas(K ∩ ∂Ω) remains bounded from below, which means that λ must grow without bound.This situation is remedied by inserting the band of elements between the cut and the boundary.In Burman and Hansbo [3], where no band was used, this conditioning problem was instead eliminated by use of additional stabilizing terms, limiting the analysis to linear polynomials.In the present formulation, we avoid the use of additional terms since the mesh on the band Ω 1 is always shape regular.

IMPLEMENTATION ASPECTS
The difference between the suggested approach and that of [7] is in the new hierarchical formulation.Conceptually, this means that in the method of [7], there are only two finite element spaces: the one on the band and the one inside the band consisting of elements up to the interface Γ.This gives two sets of unknowns, u 1 and u 2 say, and the system of equations becomes where S 1 is the stiffness matrix from the discretization on the band, S 2 from the cut mesh, B represents the coupling terms, and f 1 , f 2 represents load terms.The degrees of freedom in u 2 outside the cut can then be discarded already at the outset.One problem with this approach is that all matrices B, S 1 , and S 2 change is we want to move the interface.This could be in a time dependent problem, or if we want to use the scheme for the purpose of shape optimization.In the present approach we instead obtain the system where S 2 now denotes the stiffness matrix on the underlying mesh which does not change.Since this matrix is by far the largest of the involved matrices, this means that we only have to recompute small matrices (corresponding in a sense to the multiplier matrices of the original fictitious domain method).
In our implementation, we have used Gaussian quadrature on the interface using the band as the master mesh.We have also used the boundary on the mesh on Ω 1 to perform integration of jump and consistency terms, this boundary will not precisely match the cut mesh if the boundary is curved since we use linear cuts in the elements.This does not, however, affect the convergence rates in our example below.

NUMERICAL EXAMPLE
We consider a problem posed on a disc of radius r 0 = 0.95.With r the length of the radius vector, we use f = r to obtain the exact solution u = u = (r 3 0 − r 3 )/9.The stabilization parameter was set to γ = 10.In Fig. 5 we show the obtained convergence rates in L 2 (Ω)− and H 1 (Ω)−norms, which are optimal.An elevation of the solution is given in Fig. 6, and an elevation of the solution on the whole of Ω is shown in Fig. 7.Note that we have extended f = r to hold on the whole of Ω and imposed zero boundary conditions on ∂ Ω.This is of no consequence since the solution is decoupled at Γ.

CONCLUDING REMARKS
In this contribution we have shown that the NXFEM method is well suited for fictitious domain type simulations.It has optimal convergence for arbitrary polynomial order and does not require Lagrange multipliers to enforce Dirichlet boundary conditions.Indeed, since the boundary conditions are prescribed on a regular mesh, we can handle all types of boundary conditions in the usual way.
Applications for the proposed method include shape optimization, where the boundary of the domain has to be moved in order to calculate sensitivities, and for computations involving objects moving across a background mesh in general.

Figure 3 .
Figure 3.A discontinuous trial function and its split into a continuous and a discontinuous part.

Figure 4 .
Figure 4.An element cut by the boundary.

Figure 5 .
Figure 5. Convergence in broken H 1 and L 2 norms.

Figure 6 .
Figure 6.Elevation of the computed solution on Ω.

Figure 7 .
Figure 7. Elevation of the solution on Ω.