Multimesh finite elements with flexible mesh sizes

We analyze a new framework for expressing finite element methods on arbitrarily many intersecting meshes: multimesh finite element methods. The multimesh finite element method, first presented in Johansson et al. (2019), enables the use of separate meshes to discretize parts of a computational domain that are naturally separate; such as the components of an engine, the domains of a multiphysics problem, or solid bodies interacting under the influence of forces from surrounding fluids or other physical fields. Furthermore, each of these meshes may have its own mesh parameter. In the present paper we study the Poisson equation and show that the proposed formulation is stable without assumptions on the relative sizes of the mesh parameters. In particular, we prove optimal order a priori error estimates as well as optimal order estimates of the condition number. Throughout the analysis, we trace the dependence of the number of intersecting meshes. Numerical examples are included to illustrate the stability of the method.


Introduction
The multimesh finite element method presented in [1] extends the finite element method to arbitrarily many overlapping and intersecting meshes.This is of great value for problems that are naturally formulated on domains composed of parts, such as complex domains composed of simpler parts that may be more easily meshed than their composition.This is of particular importance when the parts are moving, either relative to each other or relative to a fixed background mesh, as part of a time-dependent simulation or optimization problem [2,3].Fig. 1 provides some illustrative examples.Here, as in [1], we consider the Poisson equation with stationary interfaces to simplify the analysis.
Another approach is to use a matching mesh and make use of elements with polytopic shapes.Methods with this capability include the PolyDG method [51,52], hybrid high order methods [53,54], virtual element methods [55,56], and mimetic methods [57,58].
The contributions of this paper is first a generalization of the formulation in [1] for the Poisson equation to allow for meshes of arbitrary mesh sizes.In the formulation, each mesh has its own mesh size and can be placed in a general position.The properties of the mesh arrangement are encoded in terms of the maximum number of overlapping meshes at any point in the domain.Naturally, this number may be much lower than the total number of meshes.The second contribution is a detailed analysis of the method.We carefully trace the dependency of the number of intersecting meshes in the coercivity of the method, in the error estimates and in the analysis of the condition number.The analysis holds for two and three dimension as well as for higher order elements, and extends previous works on cut finite elements for overlapping meshes and interface problems to much more general mesh arrangements and mesh sizes.We restrict ourselves to two dimensions in the numerical examples.See also [59][60][61] where related formulations for the Stokes problem are presented and analyzed.
In the remainder of this paper, we analyze the multimesh finite element method for the Poisson problem for an arbitrary number of intersecting meshes and arbitrarily mesh sizes, and present numerical examples.We will start with reviewing the notation from [1] in Section 2, following a presentation of the multimesh finite element method 3. We then proceed to establish standard results such as consistency and continuity of the method in Section 4. Showing coercivity, interpolation error estimates, a priori error estimates and a condition number estimate require more work, which is why we dedicate individual Sections to these in 5, 6, 7, and 8 correspondingly.We end the paper with numerical results in Section 9, conclusions in Section 10 and acknowledgments in the last section.

Notation
We first review the notation for domains, interfaces, meshes, overlaps, function spaces and norms used to formulate and analyze the multimesh finite element method.For a more detailed exposition, we refer to [1].

Notation for domains
3, be a domain with polygonal boundary (the background domain).
The method can be extended to include situations where the subdomains may intersect the boundary by using weak enforcement of boundary conditions.If cut elements appear at the boundary some stabilization of the formulation is needed, for instance face based least squares control of jumps in derivatives across faces in the vicinity of the boundary [17] or an extension procedure [62].

Notation for interfaces
• Let the interface Γ i be defined by

Notation for meshes
• Let Kh,i be a quasi-uniform [63] premesh on Ωi with mesh parameter h i = max K ∈ Kh,i diam(K ), i = 0, . . ., N (see Fig. 4(a)).(The illustrations (a) and (b) also appear in [1].) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) . ., N be the active meshes (see Fig. 4(b)).
• The multimesh is formed by the active meshes placed in the given ordering (see Fig. 5 . ., N be the active domains.

Notation for overlaps
• Let O i denote the overlap defined by 5(a) for an example.
• For i < j, let be a function indicating which overlaps are non-empty.For ease of notation, we further let δ ii = 1 for i = 0, . . ., N .
be the maximum number of overlaps.Note that N O is bounded by N but is smaller if not all meshes intersect with each other.See Fig. 5(c) for an example.
, the number of meshes below mesh i with non-empty intersection.

Notation for function spaces
• Let W s p (ω) denote the standard Sobolev spaces on ω ⊂ Ω with norm denoted by ∥ • ∥ W s p (ω) and semi-norm . The special case p = 2 is denoted by H s (ω) and the space with p = 2 and zero trace is denoted by H s 0 (ω) (see also e.g.[63,64]).The Euclidean norm on R N is denoted by | • | N .The corresponding inner products are labeled accordingly.The same notation is used for the Lebesgue measure and absolute value.It will be clear from the argument which is used.
• In order to define function spaces on the different meshes in the multimesh configuration we recall the concept of external direct sums of vector spaces, see for instance [65].Let the external direct sum X of the vector spaces X i be denoted by X = ⨁ N i=0 X i , Here X is the vector space with elements x ∈ X which are N + 1 tuples of the form x = (x 0 , . . ., x N ), where x i ∈ X i for i = 0, . . ., N , equipped with component wise addition and scalar multiplication If the vector spaces X i are equipped with inner products (•, •) X i and associated norms ∥ • ∥ X i we let Note that the direct sum is a purely algebraic construction, which is not dependent on the fact that the component spaces X i have specific further properties for instance being function spaces, and since we always have a finite number of component spaces the direct sum is identical to the direct or Cartesian product of spaces, see [65] for further details.• Let H s (Ω h,i ), s ≥ 0, be the standard Sobolev spaces of order s on the domain Ω h,i , for i = 0, . . ., N , and define the so called multimesh Sobolev space as the external direct sum (2.6) • Let V h,i be a continuous piecewise polynomial finite element space on the partition K h,i of Ω h,i , and let the multimesh finite element space V h be the external direct sum of the finite element spaces V h,i , i = 0, . . ., N .To simplify the presentation we assume homogeneous Dirichlet boundary conditions on ∂Ω with strong implementation of the boundary conditions in the finite element space, i.e., v = 0 on ∂Ω for v ∈ V h,0 .Other boundary conditions can be enforced using standard techniques.• To represent a function v ∈ H s (Ω ), s ≥ 0, as a multimesh function let E † : H s (Ω ) ↪→ ⨁ N i=0 H s (Ω h,i ) be the embedding defined by • To interpret the multimesh function v ∈ ⨁ N i=0 L 2 (Ω h,i ) as a function in L 2 (Ω ) we define an embedding E : We note that it follows from the definitions of E and E † that since where we used the fact Ω i ⊂ Ω h,i in the last step.For implementation purposes the following equivalent definition is useful which corresponds to picking the top-most mesh in the case when x belongs to several meshes.
• For s ≥ 0 we define the sum of V h and ⨁ N i=0 H s (Ω h,i ) by where each of the component spaces

Notation for jumps and averages
• To formulate the stabilization form we define a jump operator for functions v ∈ V h on the overlaps

.14)
• To formulate the Nitsche method we define the jump and average operators on the interface segment Γ i j = Γ i ∩ Ω j for i > j (cf.Section 2.2), as follows ) where v i ∈ V h,i and v j ∈ V h, j .The weights κ i and κ j are defined by and we note that Remark 2.2.By the definitions in Sections 2.2 and 2.4 we have Therefore, the two jump operators [•] on Γ i j and • on O ik are compatible on Γ i j in the sense that The indices are swapped since it is most natural to partition Γ i with respect to domains below i, and O i with respect to domains above i.We will use this to show coercivity in Proposition 5.4 and in the interpolation estimate in Proposition 6.1.

Notation for norms
• Let c > 0 and C > 0 be constants.The inequality x ≤ C y is denoted by x ≲ y.The equivalence cx ≤ y ≤ C x is denoted by x ∼ y. • Let ∥ • ∥ s h denote the semi-norm defined by and let ∥ • ∥ h denote the norm Note that for the norm ∥ • ∥ h , the domain of integration extends to each active domain Ω h,i , meaning that each overlap will be counted (at least) twice.• The definition of the energy norm, which we will denote by ||| • ||| h , will be defined after the presentation of the finite element method in (4.2).

Finite element method
As a model problem we consider the Poisson problem where Ω ⊂ R d is a polygonal domain.The multimesh finite element method for (3.1) is to find where for v, w ∈ V h , and we recall that the Dirichlet boundary condition u = 0 on ∂Ω is for simplicity enforced strongly in V h,0 .
Here, β 0 > 0 is the Nitsche (interior) penalty parameter and β 1 > 0 is a stabilization parameter.Note the relation s h (v, v) = β 1 ∥v∥ 2 s h between the stabilization term s h and the norm ∥ • ∥ s h (2.20).

Energy norm, consistency, Galerkin orthogonality and continuity
In the forthcoming analysis we will need to compare the exact solution u ∈ H s (Ω ) and the finite element solution u h ∈ V h .To facilitate the analysis we will see that it is natural to represent the exact solution u ∈ H s (Ω ) as a multimesh function in ⨁ N i=0 H s (Ω h,i ), using the E † operator (2.8), and define the error by For clarity, when there is no risk of misunderstanding we use the simplified notation u = E † u and write 2) The numbering of the terms will be used to alleviate the analysis of the method.In order for A h and the norm ||| • ||| h to be well defined also on E † u we note that the traces of the normal flux appearing in the forms are well defined if u ∈ H 3/2+ϵ (Ω ), ϵ > 0. We then have using Γ i j ⊂ ∂ Ωi and standard trace inequalities noting that Ωi does not depend on h.In view of these observations we define We note that We next establish the consistency, Galerkin orthogonality and continuity of the form A h .
Proof.This result follows by for v = (v 0 , . . ., v N ) ∈ V h , multiplying (3.1a) by v i , integrating by parts on Ω i ⊂ Ω h,i and summing the contributions For the boundary term we note the decomposition Using summation by parts and then swapping the indices i and j, the second term takes the form Now on Γ i j we obtain using the fact that n j = −n i and the definition of the jump (2.15) and average (2.16) with weights summing to one (2.18), Here all traces of the gradient are well defined in view of (4.3).Now observing that u| which combined with the observation that and similarly Proposition 4.2 (Galerkin Orthogonality).The form A h satisfies the Galerkin orthogonality; that is, ) Proof.The result follows directly from Proposition 4.1 and (3.2). □ Proposition 4.3 (Continuity).The form A h is continuous; that is, where V is defined in (4.4).
Proof.The result follows by repeated use of the Cauchy-Schwarz inequality.□

Coercivity
To prove that the form A h is coercive, we will make use of the following lemma.
Lemma 5.1.For all v ∈ V h , we have Proof.Take an element K ∈ K h,i and observe that K = ∪ N j=i K ∩ Ω j .It follows that Here we have made use of the inequality where we have used Ω h,i ∩ Ω j = O i j ⊆ Ω j for i < j.Note that the second sum is empty for i = N .Summing over all domains, we obtain by (2.2) which proves the estimate.□ Remark 5.2.There is a dependence on the maximum number of overlaps N O in Lemma 5.1.In practice, N O is of moderate size and this dependence is not an issue.The interpolation error estimates and condition number estimates (shown below) have a different kind of dependence.
Remark 5.3.Using an inverse bound of the form (see e.g.[63]) one can show that the stabilization term s h (v, w) may alternatively be formulated as Using Lemma 5.1, we may now proceed to prove the coercivity of the bilinear form.
Proposition 5.4 (Coercivity).The form A h is coercive.More precisely, for β 0 and β 1 large enough, we have (5.12) Proof.We first note that (5.13) Now, for l = i or l = j, let denote the set of elements in K h,l which intersect Γ i j .Using an inverse estimate (see [12]), we have where the constant is independent of the position of Γ i j .It follows that (5.16) (5.17) where we have noted that Γ i j (the part of Γ i bordering to Ω j for j < i) is empty if the overlap O ji (the part of Ω h, j intersected by Γ i for j < i) is empty, as indicated by δ ji .See also Remark 2.2.
For l = i, we thus obtain the estimate while for l = j, we obtain the similar estimate = N O ∥∇v∥ 2 h . (5.27) Proceeding with ⋆ using the Cauchy-Schwarz inequality with weight ϵ(h i + h j ), where ϵ is a positive number, and Young's inequality 2ab ≤ a 2 + b 2 we obtain (5.32) In (5.31) we used that (5.34) where we used the definition (2.17) of the weights κ l to obtain In (5.32) we made use of (5.23) and (5.27).By Lemma 5.1, we may now estimate the ∥ • ∥ h norm in terms of the ∥ • ∥ Ω norm and the ∥ • ∥ s h norm to obtain Combining (5.13) and (5.37), we find that and thus, by choosing ϵ small enough and then β 0 and β 1 large enough, The coercivity now follows by noting that term I I I in (4.2) may be controlled by terms I and I I as above in the estimate of ⋆. □ Remark 5.5.By continuity (4.17), coercivity (5.12) and a continuous l h (v), there exists a unique solution to (3.2) by the Lax-Milgram theorem (see e.g.[63]).

Interpolation error estimate
To construct an interpolation operator into V h , we pick a standard interpolation operator into V h,i , where π h,i satisfies the standard interpolation error estimate (see e.g.[63]) Here, N h (K ) denotes the set of elements that share a vertex with K .We then define the interpolation operator which by composition with E † , see (2.8), provides the interpolation operator To prove an interpolation error estimate for π h , recall Remark 2.1 regarding the non-intersecting boundaries and let U δ (Γ i j ) denote the tubular neighborhood of Γ i j defined by where B δ (x) is the ball of radius δ centered at x; see Fig. 6.In addition, let Proposition 6.1 (Interpolation Error Estimate).The interpolation operator π h satisfies the interpolation error estimate where and the norm is defined by (6.9) Proof.We first let η = v −π h v denote the interpolation error and recall the numbering of the terms in the definition of the energy norm ||| • ||| h (4.2).Starting with term I , we have For term I I , we have, since For term I I I , recall the inverse estimate (5.16) and note that K h, j(Γ i j ) ⊆ U δ (Γ i j ) with δ ∼ h j .Thus, For term I V , we first handle the jump term as in I I and then proceed as for I I I , Therefore, there are only two terms in I -I V that need to be estimated.First which follows immediately by (6.2).Second, we make use of the disjoint partition of Γ i and noting that we obtain Due to the maximum norm, this estimate also holds with η i replaced by η j , and the desired estimate holds.□ Remark 6.2.Note that the stronger control v i ∈ W k+1 ∞ (U h (Γ )) enables us to establish the estimate (6.7) with the constant given by (6.8) which only have weak dependence on the configuration of the overlapping mesh configuration encoded in terms of N O i and |Γ i |.

A priori error estimates
We may now prove the following optimal order a priori error estimates.The estimates are supported by the numerical results presented in Fig. 7.For details on these results, we refer to the accompanying paper [1].Fig. 7. Rate of convergence in the L 2 (Ω ) (left) and H 1 0 (Ω ) (right) norms for p = 1 (blue), p = 2 (red), p = 3 (yellow) and p = 4 (purple), where p is the polynomial degree of the finite element approximation.For each p, the convergence rate is shown for N = 1, 2, 4, 6, 16, 32 meshes (six lines) and the errors for N = 0 (the standard single mesh discretization) are marked with × and dashed lines.(These illustrations also appear in [1].) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Theorem 7.1 (A priori Error Estimates).The finite element solution u h of (3.2) satisfies the following a priori error estimates: Proof.The proof of (7.1) follows the standard procedure of splitting the error and using the energy norm interpolation error estimate from Proposition 6.1, To estimate the second term on the right-hand side of (7.3), we use the coercivity (Proposition 5.4), Galerkin orthogonality (Proposition 4.2) and continuity (Proposition 4.3) of A h to obtain It follows that Combining (7.3) and (7.7) with the interpolation error estimate of Proposition 6.1 now yields (7.1).
To prove (7.2), we use a standard duality argument (see e.g.[63]).Let φ be the solution to the dual problem with ψ ∈ L 2 (Ω ).Using elliptic regularity we then have from which it follows that φ ∈ V .Using the fact that A h is symmetric it follows from consistency (Proposition 4.1) that We now take v = e = u − u h and use the Galerkin orthogonality (Proposition 4.2), continuity (Proposition 4.3) and a standard interpolation inequality on each set Ω h,i (note that we cannot use stronger regularity than φ ∈ H 2 (Ω ) since ψ ∈ L 2 (Ω ) and thus the interpolation bound 6.1 is not applicable for φ) to obtain (e, ψ) Ω = A h (e, φ) (7.11) ) 1/2 (7.16) where we used the fact that the maximum number of overlapping meshes is N O and in the last step we have used the standard elliptic regularity estimate (see e.g.[63]).Note that we have continuity (7.13) also for functions in H 3/2+ϵ (Ω ), ϵ > 0, as noted in Proposition 4.3.The desired estimate (7.2) now follows from (7.17) by (7.1) and taking ψ = e.□

Condition number estimate
To prove a bound on the condition number, we first introduce some notation and definitions.Let {ϕ i, j } M i j=1 be the finite element basis of V h,i .We then have the expansion for each part v i of a multimesh function v = (v 0 , . . ., v N ).Collecting all expansion coefficients for the 1 + N parts into a vector v of dimension M = ∑ N i=0 M i , the total stiffness matrix Â for the multimesh system is defined by with condition number To derive an estimate of κ( Â) we make use of the following lemmas.
Lemma 8.1 (Inverse Inequality).It holds that Proof.Recall the definition of the energy norm (4.2).We first note that For term I I , we have by recalling (6.12) Term I I I may be estimated similarly to obtain For term I V , we have by recalling (6.16) The desired estimate now follows using the standard inverse inequality (5.10).□ Lemma 8.2 (Poincaré Inequality).It holds that where Proof.First note that by a Taylor expansion argument and Lemma 5.1, we have To control the first term on the right-hand side in (8.14), let φ ∈ H 2 (Ω ) be the solution to the dual problem where ψ ∈ L 2 (Ω ).Multiplying the dual problem with v ∈ V h and integrating by parts, we obtain using the Cauchy-Schwarz inequality Now we continue with the second factor in (8.19).Using the trace inequality with constant independent of the position of an interface γ (see [12]), we have By the construction of U δ (Γ i j ), see (6.5), we have K h,i (Γ i j ) ⊆ U δ (Γ i j ) with δ ∼ h i .Furthermore, by the Hölder inequality [64] with coefficients r, s such that 1/r + 1/s = 1 we have with p = 2r and 1/s = 1 − 1/r = 1 − 2/ p.To determine p in (8.28), we use the Sobolev embedding W l q (Ω ) ⊆ W k p (Ω ) [64] with k = 1, l = 2 and q = 2.This is motivated by the fact that due to elliptic regularity and ψ ∈ L 2 (Ω ), we have φ ∈ H 2 (Ω ).Since the embedding holds for 1/ p − k/d = 1/q − l/d [64], we obtain p = 2d/(d − 2), where p = ∞ for d = 2. Thus Cf. [64] regarding the last inequality for d = 2, 3. Returning to the second factor in (8.19) we thus have, using (8.21), (8.28) and (8.30) together with a standard duality argument (see e.g.[63]), elliptic regularity, a stability estimate, the Poincaré equality and ignoring higher order terms that The bound on ∑ i (v i , ψ) Ω i in the left-hand side in (8.16) with ψ = v now reads ) 1/2 (8.37) since ∥v∥ Ω is bounded.
To conclude, recall (8.14) and insert (8.37) to obtain the desired estimate Proof.Since K h,i is conforming and quasi-uniform we have the equivalence see e.g.[63].It follows that Recall the definition of the matrix norm The inequality thus reads Dividing by |ŵ| M and using the definition of the matrix norm (8.43) now gives By using (8.49) and (8.57) in the definition of the condition number (8.3), we obtain the desired estimate (8.40).□ The estimate for the condition number is supported by the numerical results presented in Fig. 8.The slope is found to be −1.76.The details on this example is found in [1].

Numerical results
To demonstrate the applicability and robustness of the multimesh finite element formulation, we present here a couple of numerical examples.For additional examples, we refer to the companion paper [1].

Convergence under variable mesh size
For the first example, we construct two multimesh configurations I and I I , each consisting of three parts (overlapping meshes) as show in Fig. 9.We consider a simple Poisson problem with analytical solution u(x, y) = sin(π x) sin(π y). (9.1) The goal is to study the convergence under refinement of the three meshes for each of the two test cases.Starting from initial coarse meshes with equal mesh sizes, we refine each part separately, using 8 different mesh sizes, and compute the L 2 (Ω ) and H 1 0 (Ω ) error norms.A piecewise linear finite element basis is used for both configurations.The refinement procedure is as follows.First we will refine part 0 in 8 steps, then part 1 in 8 steps, and finally part 2 in 8 steps.Then we swap the order and refine part 1 first, followed by parts 0 and 2. We do this for all permutations of the order of the parts; in total there are 3! combinations.This procedure is performed for both I and I I .Configuration I is a nested configuration (but not hierarchical).Configuration I I is generated by placing the second and third parts in a "random" position on top of the background mesh of Ω 0 = [0, 1] 2 .Specifically, we have ) ) as illustrated in Fig. 9.The meshes are refined with mesh sizes 2 −k , k = 3, . . ., 10, i.e., 8 steps.Thus, the mesh size ratio between two parts are in this example at most 2 7  = 128.In Fig. 10 we show the L 2 (Ω ) and H 1 0 (Ω ) errors for the two configurations.As expected, the different curves start and end in the same point.Moreover, we see that during refinement of the first part, errors decrease but flatten.This is due to the fact that the errors from the other two, unrefined, parts dominate.When the second part starts being refined, the errors drop but will again flatten since the errors are dominated by the third and last unrefined part.Refining this part results in a sharp decrease in the error.Due to the effect of dominating errors from different parts, we observe two L-shaped drops for each refinement permutation, for both configurations and for both error quantities.
There is no significant L-shape decrease for the first refined part, but it would be possible to construct a multimesh configuration such that this would be the case.For the example with this analytical solution, the first part would have a dominating error if the area of the part would be dominating.
It is worth noting is that the errors decrease smoothly and the method is stable despite the large differences in mesh size.

Boundary layer resolution
To demonstrate the potential of the multimesh formulation for local adaptation, we consider the boundary value problem For ϵ → 0, the PDE reduces to u = 0 which is compatible with the boundary condition on Γ 0 .As a consequence, the solution for small ϵ is u ≈ 0 away from the boundary Γ 1 and then an exponential transition to u = 1 close to Γ 1 .The width of the boundary layer is ∼ ϵ.The multimesh finite element formulation is identical to (3.2) with the We consider a model problem where the domain Ω is defined by [0, 1] 2 \ ω, where ω is the shape of the standard NACA 6409 standard.We let Γ 0 be the boundary of the unit square and let Γ 1 be the boundary of the airfoil.The solution exhibits a boundary layer of width ϵ on the airfoil boundary.
To discretize the problem, we let Kh,0 be a uniform mesh of the unit square with mesh size H = 2 −(6+k) and let Kh,1 be a boundary-fitted mesh of width w = 0.1 • 2 −k for k = 0, 1, 2, 3, 4. The boundary layer parameter is chosen as ϵ = w/2.The mesh size h ≪ H is chosen to well resolve the boundary layer.Note that we intentionally take w small relative to the boundary layer width so as not to get the entire boundary layer transition on the finer mesh, in order to illustrate better the robustness of the method and the coupling of the solution represented on the background mesh and the boundary-fitted mesh on the interface Γ .If instead we take ϵ = w/10, the solution would transition quickly to u ≈ 0 on the interface Γ .Fig. 11 shows the solution for k = 1, 2, 3, 4, clearly demonstrating the decreasing width of the boundary layer with increasing k.Note the smooth transition of the solution going from the representation on the coarse background mesh to the fine boundary-fitted mesh.In Fig. 12, a 3D view is plotted for both solution components for k = 0, 1, 2, 3, 4. Finally, Fig. 13 shows detailed plots of the solution close to the boundary layer for k = 0 and k = 4.

Conclusions
We have analyzed a general framework for discretization of the Poisson equation posed on a domain defined by an arbitrary number of intersecting meshes with arbitrary mesh sizes.The analysis show that for sufficiently large Nitsche and stabilization parameters, the method is optimal and stable.As expected, there is a dependency on the maximum number of intersecting meshes in the coercivity, error analysis and in the condition number estimate.This   was seen numerically in the accompanying paper [1], and here we are able to quantify this dependence.In addition, the numerical results presented in this paper show that the method is indeed stable when the meshes involved have vastly different mesh sizes.
As mentioned in the introduction, the multimesh method may be advantageous in the case of dynamic domains, since remeshing may be avoided.This is due to the fact that the computational geometry routines automatically identify the elements constituting the active meshes, and this can easily be done every time the domains move.Although so far only studied for two-dimensional problems [2] reports a speed up.The same approach is also applied in [3].
Future work involves extending the implementation to include three-dimensional meshes, which is a challenge due to requirements of efficient and accurate computational geometry routines in the case of arbitrary many intersecting meshes.That the multimesh formulation is valid in the case of two meshes in three dimensions is explored in [60] for the Stokes problem.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1. (Top left) The flow around a propeller may be computed by immersing a mesh of the propeller into a fixed background mesh.(Top right) The geometry of a composite object may be discretized by superimposing meshes of each component.(Bottom)The interaction of a set of solid bodies may be simulated using individual meshes that move and intersect freely relative to each other and a fixed background mesh.(These illustrations also appear in[1].)

Fig. 5 .
Fig. 5. (a) Given three ordered triangles K 0 , K 1 and K 2 , the overlaps are O 01 in green, O 02 in red and O 12 in blue.(b) The multimesh of the domains in Fig. 2(b) consists of the active meshes in Fig. 4(b).(c) Example with N = 3 domains and N O = 2 intersecting meshes.(Theillustrations (a) and (b) also appear in[1].)(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Configuration I (left) and and configuration I I (right) exemplified by refined part 2 for I and refined part 1 for I I .The coarse and fine meshes have mesh sizes 2 −3 and 2 −6 respectively.

Fig. 10 .
Fig.10.Errors during refinement of multimesh configurations I (left column) and I I (right column).Colors and markers indicate which part is being refined.The refinement procedure starts with all parts having a mesh size of 2 −3 resulting in approximately 10 2 degrees of freedom.Each part is then refined individually and sequentially as described in the text, until all parts have a mesh size of 2 −10 , resulting in a total of approximately 10 6 degrees of freedom.

Fig. 12 .
Fig. 12. Solution of the boundary layer problem (9.6) for k = 0, 1, 2, 3, 4. The left column shows the solution represented on the background mesh, the middle column shows the solution on the overlapping boundary-fitted mesh and the right column shows the composite multimesh solution.

Fig. 13 .
Fig. 13.(Top) A 3D view of the matching solutions to the boundary layer problem (9.6) on the background mesh and the overlapping boundary-fitted mesh for k = 0. (Middle) The corresponding 2D view for k = 0. (Bottom) A detailed zoom close to the tip of the airfoil for the finest mesh (k = 4).