Efficient implementation of finite element methods on nonmatching and overlapping meshes in 3D
Abstract
In recent years, a number of finite element methods have been formulated for the solution of partial differential equations on complex geometries based on nonmatching or overlapping meshes. Examples of such methods include the fictitious domain method, the extended finite element method, and Nitsche’s method. In all of these methods, integrals must be computed over cut cells or subsimplices which is challenging to implement, especially in three space dimensions. In this note, we address the main challenges of such an implementation and demonstrate good performance of a fully general code for automatic detection of mesh intersections and integration over cut cells and subsimplices. As a canonical example of an overlapping mesh method, we consider Nitsche’s method which we apply to Poisson’s equation and a linear elastic problem.
verlapping mesh, nonmatching mesh, Nitsche method, discontinuous Galerkin method, immersed interface, XFEM, algorithm, implementation, computational geometry
1 Introduction
A fundamental problem in computational science is to develop methods for the solution of partial differential equations on domains containing one or several objects that may have complex or timedependent geometry. One approach to attacking this problem is to allow overlapping meshes where a mesh representing an object is allowed to overlap a background mesh representing the surroundings of the object; see for instance Chesshire and Henshaw [22], Yu [60], Zhang et~al. [61], Mayer et~al. [50] for various applications.
In solid mechanics, overlapping meshes may be used to represent materials consisting of elastic objects inserted into a surrounding elastic material of another type [37], and in fluid–structure interaction, an overlapping mesh may be used to represent an elastic body immersed in a fluid represented by a fixed background mesh [60, 7, 49]. Another common application [22, 26, 25] is found in mesh generation where a complicated geometry such as, for example, a pipe junction, may be decomposed into simpler parts and one unstructured tetrahedral mesh is created for each part. These components may then be stored, reused and recombined in applications by using an overlapping mesh technique.
Overlapping mesh techniques are of particular interest in simulations that involve moving objects. For such problems, overlapping mesh techniques are an attractive alternative to ALE techniques. The main advantage is that by using an overlapping mesh technique, one avoids deformation of the mesh that may lead to deterioration of the mesh quality and ultimately force remeshing. This is of particular importance in the simulation of fluid–structure interaction where the topology of the fluid domain may change due to deformation of the solid.
The main focus of this work is on the general algorithms and efficient implementation that is required to handle complex problems posed on overlapping meshes in three space dimensions. Many of the presented algorithms and the tools developed are of interest and use for the implementation of various overlapping mesh techniques. To make the discussion concrete, we here focus mainly on Nitsche’s method. In Hansbo et~al. [33], a consistent finite element method for overlapping meshes based on Nitsche’s method was introduced and analyzed. The basic idea is to construct a finite element space by taking the direct sum of the space of continuous piecewise polynomial functions on the overlapping mesh and the restriction of the space of continuous piecewise polynomial functions to the complement of the overlapping mesh, and then impose the interface conditions using Nitsche’s method. It was shown that this approach leads to a stable method of optimal order for arbitrary degree polynomial approximation.
The main challenge in the implementation is to compute the intersection between the overlapping and the background mesh. The result is a set of cut cells which may be arbitrarily complex polyhedra. These arise as the result of subtracting from the tetrahedra of the background mesh a set of tetrahedra from the overlapping mesh. By adopting algorithms and search structures from the field of computational geometry, we show how these issues can be handled in an efficient manner. Furthermore, one needs to compute integrals on the resulting polyhedra. This can be carried out efficiently based on an application of the divergence theorem in combination with potentials to represent an integral on the threedimensional polyhedron as a sum of onedimensional integrals on its edges.
The presented algorithms and implementation are relevant for several other types of related methods, including the extended finite element method (XFEM) [58, 31], nonfitted sharp interface methods [32, 14], and meshtying techniques [25].
1.1 Major contributions of this paper
Our work consists of several contributions. We identify the major techniques used in the implementation of overlapping mesh methods and related methods. We further identify useful data structures and algorithms from the field of computational geometry. As part of our work, existing computational geometry libraries such as CGAL [1] and GTS [2] have been wrapped into the general purpose finite element library DOLFIN [43] and into an extension library on top of it, thereby making these algorithms and data structures more easily accessible to the finite element community. Based on our implementation, we demonstrate for the first time a highly efficient implementation of Nitsche’s method on overlapping meshes for several problems posed in three spatial dimensions, thus opening the possibility of employing Nitschebased overlapping mesh methods for challenging 3D problems such as fluid–structure interaction or domainbridging problems.
1.2 Outline of this paper
In Section 2, we review Nitsche’s overlapping mesh method for a model problem and present in Section 3 the techniques and algorithms we have developed for efficient implementation of Nitsche’s method in three space dimensions. The corresponding implementation and data structures are described in Section 4. Finally, we present in Section 5 numerical examples that demonstrate the convergence of the numerical solution as well as the scaling of the work required to compute mesh intersections relative to standard finite element assembly on matching meshes.
2 Nitsche’s method on overlapping meshes
We here review Nitsche’s method for a simple model problem posed on two overlapping meshes.
2.1 Model problem
Let be a domain in , with boundary , consisting of two (open and bounded) subdomains and separated by the interface . We consider the following elliptic model problem: find such that
(1)  
(2)  
(3)  
(4)  
(5) 
We here use the notation for the restriction of a function to the subdomain for , is the unit normal to directed from into , and denotes the jump in a function over the interface . In addition, the boundary is divided into two subdomains and where Dirichlet and Neumann boundary conditions are applied, respectively.
2.2 Finite element formulation
We consider a situation where a background mesh is given for and another mesh is given for the overlapping domain (see Figure 1). Both meshes are assumed to consist of shaperegular tetrahedra . The mesh covering is constructed by
(6) 
Note that the consists of both regular and cut elements since the mesh is not a conform tetrahedralization of the subdomain in general. For a cut element , the degree of freedoms are defined to be the same as for the original uncut element . This is illustrated in Figure 2.
We let denote the space of continuous piecewise fixedorder polynomials on that vanish on for . The space is constructed as the restriction to of functions in ; that is, . We may then define a finite element space for the whole domain by
(7) 
Nitsche’ method proposed by Hansbo et~al. [33] then takes the form: find such that
(8) 
where
(9)  
(10) 
Here, is a positive penalty parameter and the average is chosen to be the oneside derivative according to Hansbo et~al. [33]. But as pointed out by the authors, any convex combination of the normal derivatives leads to a consistent formulation.
2.3 Summary of theoretical results
When the penalty parameter in (9) is chosen large enough, the form is coercive on the discrete space and one can derive optimal order a priori error estimates in both the energy norm and the norm for polynomials of arbitrary degree . The estimates take the form
(11)  
(12) 
See Hansbo et~al. [33] for details. Here, denotes the standard Sobolev norm of order on . A posteriori error estimates and adaptive algorithms are also presented in Hansbo et~al. [33].
As observed by Areias and Belytschko [6], Nitsche’s method for an arbitrary cutting interface as in Hansbo and Hansbo [32, 34] can be reinterpreted as a particular instance of an extended finite element method. In contrast, the presented formulation on overlapping meshes lacks such a reinterpretation since two unrelated meshes are involved and therefore the discontinuity across the interface cannot be modeled by the enrichment of degrees of freedom within a single cell. Nevertheless, both methods share some features, especially with regard to their implementation.
3 Techniques and algorithms
The main challenges in the implementation of Nitsche’s method on overlapping meshes arise from the geometric computations which are necessary to assemble the discrete system associated with (8)–(10). Naturally, these geometric computations are more involved in 3D than in 2D. Similar challenges are encountered in related methods such as the extended finite element method [31], but also for the simulation of contact mechanics [59]. Different solutions have been proposed, see for example Sukumar et~al. [58] in the case of extended finite elements methods, and Yang and Laursen [59] for contact problems. Here, we take another approach which efficiently solves the issues in the case of the Nitsche overlapping mesh method. In the following, we will first discuss the challenges and their remedies in general terms, and then return to the specific details of our implementation in Section 4.
3.1 Assembly on overlapping meshes
The main implementation challenges arise from the fact that the interface can cross the overlapped mesh in an arbitrary manner, which has two consequences. First, the definition of the finite element space (7) involves the restriction of the function space , defined on the full background mesh of , to the domain . The restriction results in nonstandard element geometries along the interface as it can be observed from the definition 6 of . Secondly, the weak imposition of the interface conditions (2) and (3) by the interface integrals in (9) involves finite element spaces defined on two unrelated meshes. Both make the assembly challenging to implement, compared to a standard finite element method.
To better understand what kinds of challenges the Nitsche overlapping mesh method adds, we briefly recall the general theme of the assembly routine as it is realized nowadays in many finite element frameworks; see, e.g., Logg [41], Bastian et~al. [12]. Nitsche’s method is closely related to the classical discontinuous Galerkin (DG) methods, which makes it natural to depart from the assembly in the DG case. A detailed description of finite element assembly for discontinuous Galerkin methods can be found in Ølgaard et~al. [55]. We assume a variational problem of the following form: find such that
(13) 
where and are bilinear and linear forms, respectively, and and are the discrete trial and test spaces, respectively. For simplicity, we here make the assumption . The solution of the variational problem (13) may then be computed by solving the linear system
(14) 
for , with stiffness matrix and load vector . Here, denotes a basis for . During assembly of the tensors and , one usually iterates over all cells and computes the contributions from each cell as the cell tensors and where is the local finite element basis on (the shape functions). The element tensors are then scattered into the global tensors and by adding the entries according to a given localtoglobal mapping . One may similarly add contributions from each facet of the set of boundary facets (the set of exterior facets) and, for the implementation of a DG method, from each facet of the set of interior facets . Algorithm 1 summarizes the parts of the standard assembly algorithm relevant for our discussion.
for each : 
for each 
for each 
for each : 
for each : 
for each : 
for each : 
for each : 
for each : 
We now consider how the standard assembly algorithm must be modified in the case of Nitsche’s method on overlapping meshes. We first note that the tessellation of the background domain may be decomposed into three disjoint subsets:
(15) 
where , , and denote the sets of not, completely and partially overlapped cells relative to , respectively.
Integrals over the cells of can be assembled using a standard assembly algorithm. Furthermore, integrals over the cells of need not be assembled (the corresponding contributions will be assembled over ). However, assembly must be carried out over , the partially overlapped cells of . For the model problem (8)–(10), this requires the evaluation of integrals of the type
(16) 
on cut elements (polyhedra) where . Examining (8)–(10), we further note that we must assemble the terms
(17) 
This poses an additional challenge, since the integrands involve products of trial and test functions defined on different meshes. Furthermore, the interface consists of a subset of the boundary facets of , but each such facet may intersect several cells of . We therefore partition each facet on into a set of polygons such that each polygon intersects exactly one cell of the overlapping mesh and one cell of the background mesh . This is illustrated in Figure 3. Assembly may then be carried out by summing the contributions from each polygon .
In summary, we identify the following main challenges in the implementation of Nitsche’s method on overlapping meshes:

collision detection: to determine which cells are involved in the intersection between the two meshes;

mesh intersection: to compute the cut cells of the background mesh represented by the polyhedra and the intersection interface represented by the polygons ;

integration on complex polyhedra: to compute integrals on the polyhedra (and the polygons ).
In the following, we discuss these challenges in some detail and introduce concepts, data structures and algorithms to handle them in an efficient manner. We emphasize that the proposed solutions are not limited to the implementation of Nitschetype methods, but may also be used for the implementation of related overlapping mesh methods.
3.2 Collision detection
Topological relations between entities of a single mesh can be described by concepts of connectivity or mesh incidence as presented in Logg [42] and Bastian et~al. [11]. To represent the topological relation between two overlapping (colliding) meshes, we enrich the notation of Logg [42] by the concepts of collision relations and collision maps. The collision relation between and is defined by
(18) 
The collision relation lists pairs of indices of all intersecting cells of the background mesh with boundary facets of the overlapping mesh . Furthermore, each index of an intersected entity is mapped to the set of indices of intersecting entities via a pair of collision maps:
(19)  
(20) 
We note that the collision maps and can be computed from the collision relation .
A naive approach to computing the collision relation between two meshes and would be to intersect each cell of with each cell of , resulting in an complexity, which is not feasible for large meshes. However, efficient algorithms and data structures which reduce the complexity dramatically have been developed in the fields of computer science, computational geometry, and computer graphics. The task of determining whether two objects collide arises naturally in the rendering of a computer game and the objects in question are often represented by meshes. In this paper, we limit the description to those techniques used in our work and refer to the books Ericson [27], Schneider and Eberly [57], and AkenineMöller et al. [3] for a broader overview. To efficiently find all pairs of the collision relation , two important concepts are used: (i) fast intersection tests for pairs of simple geometric objects and (ii) spatial data structures for large sets of objects to accelerate collision queries.
3.2.1 Fast intersection detection
Fast intersection detection for simple geometric entities such as triangles and tetrahedra is based on socalled geometric predicates. Geometric predicates are tests which determine whether two geometric entities do intersect (collide) without computing the actual intersection. See Figure 4 for an example. This saves work since the actual intersection does not need to be computed when the geometric predicate is false.
3.2.2 Spatial data structures
Spatial data structures provide a hierarchical ordering of geometric objects, which allows quick traversal of large sets of geometric objects as part of a collision test. There exist two principal approaches based on either space subdivisions or bounding volume hierarchies, which can both be realized as treelike data structures. As the name indicates, subdivision approaches rely on some sort of geometric subdivision of the entire space embedding the structures of interest, in our case a finite element mesh. Typical examples are binary space partitions (BSP) trees, quadtrees, and octrees. Since the embedding space is subdivided, the leaf of a tree will usually not represent a single mesh entity and therefore further selection procedures are required to determine the actual intersections.
In contrast, a bounding volume hierarchy is a tree which is built from bounding volumes; that is, simple geometric shapes containing the objects to be tested. Examples of such data structures are axis aligned bounding boxes (AABB, see Figures 4 and 5), oriented bounding boxes (OBB) and socalled kDOPs, discrete orientation polytopes described by hyperplanes. The purpose of using bounding volumes is twofold: (i) testing for collision of bounding volumes is cheaper than testing for collision of the bounded objects, and (ii) the simple geometry of the bounding volumes means that they can be stored efficiently in a hierarchical manner.
3.2.3 Building the collision map
Bounding volume trees accelerate asymmetric collision queries between a tree embedding a large set of objects and a single simple object like a tetrahedron. The concept becomes even more powerful when intersection tests between two large sets of geometric primitives are desired, as in our case of two meshes. Then a hierarchical traversal of both trees greatly improves the complexity of the naive approach. Algorithm 2 describes how to employ a hierarchical traversal to compute the collision map.
compute_collisions(, ): 
if : 
return 
else if is_leaf() is_leaf(): 
if : 
else if descend_a(, ): 
for children(): 
compute_collisions(, ) 
else: 
for children(): 
compute_collisions(, ) 
descend_a(, ): 
return is_leaf() is_leaf() 
After the completion of the traversal according to Algorithm 2, which identifies all cells of that intersect the boundary of , it remains to check whether the remaining cells (which do not intersect the boundary of ) are either completely overlapped or not overlapped by . To check this, it is enough to take a single point contained in each cell and check whether it is in . To avoid building a third AABB tree for the entire mesh and to take advantage of the already built tree for , one can use a method known as rayshooting; see AkenineMöller et al. [3]. The idea is simple: counting how often a ray starting at the point in question intersects the surface tells whether it is inside (odd number of intersections) or outside (even number of intersections). To efficiently find all ray–surface intersections, one may reuse the same AABB tree for that was used in Algorithm 2.
To summarize, the collision map and classification of cells according to the splitting (15) can be found efficiently by utilizing AABB trees, fast traversal of pairs of AABB hierarchies, and ray shooting techniques.
3.3 Mesh intersection
The next step is to compute the cut cells (polyhedra) and the interface decomposition ; see Figure 3. This computation may be phrased in terms of socalled boolean operations which are widely used in CAD systems to build complex geometries by performing boolean operations between primitives from a finite set of geometries. Algorithm 3 summarizes the use of boolean operations to compute and . The resulting geometric objects are depicted in Figure 6.
for each intersected by : 
for each : 
for each : 
for each : 
The boolean operations are completely delegated to the computational geometry library GTS [2] which uses a socalled earclipping algorithm to compute the surface tesselation of the intersection objects. The implemented ear clipping algorithms (first introduced in Fournier and Montuno [28]) are known to have complexity, where is the number of vertices of the 3Dpolygon. This number in turn () depends on the mesh quality parameters such as the smallest and widest element angles and will in practice be bounded by a constant.
3.4 Integration on complex polyhedra
The last required step in implementing the overlapping mesh method is to compute the integrals (16) and (17) on the cut cells and on the interface decomposition , respectively. A widely used technique in the implementation of the extended finite element method [52, 24, 58, 15] is to decompose a cell into subcells which are aligned with the interface. The interface is either approximated by a plane within the cell [29, 30] or completely recovered [24, 58], even in the case of higher order interfaces [49]. Inside each tetrahedron of the subtetrahedralization, one may then use a standard integration scheme. However, subtetrahedralization of an arbitrary polyhedron is in general a quite challenging problem [56], and its existence cannot even be guaranteed without adding additional vertices [21], in contrast to the twodimensional case.
Methods for integration over arbitrary polygonal domains without the use of a subtriangulation have recently been presented in Natarajan et~al. [53] and incorporated into the extended finite element method with discontinuous enrichments in Natarajan et~al. [54]. Using Schwarz–Christoffel conformal mappings as the fundamental tool, the technique is strongly bounded to two space dimensions and hard to generalize to three space dimensions.
We here propose an alternative approach, which is based on a boundary representation of the integrals. Despite the fact that this technique has been known for a long time, see for instance [40], it seems to be largely unknown within the finite element community. Our implementation is based on the efficient realization described in Mirtich [51]. This technique may be easily generalized to the computation of the general moment integral over a polyhedron , defined by
(21) 
where denotes a multiindex of length .
The integral is computed in three steps. The first step is to interpret the integrand as the divergence of a polynomial vector field and to rewrite (21) as surface integral:
(22) 
Here, denotes the th unit vector.
Secondly, the plane equation for each facet allows the construction of a projection map , where is some positive, orientationpreserving permutation of . Using the parametrization , one may rewrite the integrals of type in (22) as integrals in the plane:
(23) 
Here, denotes the component of the normal vector .
The third step consists of using a parametrization of and Green’s theorem in the plane to rewrite (23) as a sum of line integrals. Together, the three steps reduce the evaluation of the moment integral (21) to the evaluation of a set of onedimensional integrals with polynomial integrands. We note that the algorithm for the calculation of each moment integral has only a complexity.
The moment integrals can be used in several ways to finally compute the integrals (16) and (17). First, one may compute the volume and the barycenter as the zeroth and first moments, respectively, which immediately provides a quadrature rule for exact integration of polynomials of degree one. Unfortunately, higher order moments may not be directly reinterpreted as quadrature rules. Furthermore, as the construction of quadrature rules for higher degree polynomials can be quite involved [23], it can be advantageous to avoid runtime generation of quadrature rules. We therefore propose to use the moment integrals directly by first interpolating the integrand onto a monomial basis and then summing the contributions:
(24) 
Illconditioning may be avoided for higherorder expansions by replacing the monomial basis by Bernstein polynomials. In the present work, linear Lagrange elements have been used throughout and so simple barycenter quadrature suffices.
4 Implementation and data structures
We now discuss some of the data structures and classes which reflect the abstract concepts and algorithms described in Section 3. For some of these algorithms, we rely on existing implementations as part of the computational geometry libraries CGAL [1] and GTS [2], while other algorithms have been realized in the finite element library DOLFIN [43, 46] as part of this work. Specialized algorithms for the Nitsche overlapping mesh method have been implemented as part of the extension library DOLFINOLM which is built on top of DOLFIN. The code is free/opensource, licensed under the LGPLv3, and available at http://launchpad.net/dolfinolm.
4.1 The finite element library DOLFIN
Our implementation of Nitsche’s method on overlapping meshes is based on the finite element library DOLFIN which is part of the FEniCS project [44, 41] for automated scientific computing. The main feature of FEniCS is the automated treatment of finite element variational problems, based on automated generation of highly efficient C++ code from abstract highlevel descriptions of finite element variational problems expressed in nearmathematical notation [39, 4]. This is combined with builtin tools for working with efficient representations of computational meshes [42] and wrappers for highperformance linear algebra libraries like PETSc [10, 9, 8] and Trilinos [36].
For this work, we have integrated the computational geometry libraries CGAL [1] and GTS [2] with DOLFIN. While CGAL and GTS provide a large part of the functionality needed to compute intersections between tetrahedra, the integration scheme was realized through an adapted version of Mirtich’s code [51]. Furthermore, the assembly routine of DOLFIN was extended to handle integration over cut cells and meshes. Currently, the code for computing local interface integrals has to be implemented manually by the user, but we plan to extend FEniCS, in particular DOLFIN and the form compiler FFC [39, 45], to provide a full automation of Nitsche’s method where a user only needs to supply the abstract variational problem (8).
4.2 The class AABBTree
The search data structure AABBTree has been added to the DOLFIN library. The implementation is based on the computational geometry library CGAL [1]. Basic search queries such as finding one or all cells intersecting a given entity or distance computation are exposed via an IntersectionOperator class. DOLFINOLM complements that functionality with providing a GTSbased AABB tree [2] to allow traversal of two bounding box trees as described in Algorithm 2.
4.3 The class OverlappingMeshes
The class OverlappingMeshes, provided as part of DOLFINOLM, is a key component in our realization of the overlapping mesh method. It mainly computes additional topological and geometric information to describe the overlap of the two meshes and and provides access to this information, in particular the collision relation and the collision maps and described in Section 3.2.
A MeshFunction, as introduced in Logg [42], describes the splitting (15) of the mesh by assigning different integer values to cells in the not overlapped part , the completely overlapped part , and the partially overlapped part , respectively. In the same way, the boundary facets of the overlapping mesh can be marked if the overlapping domain is not completely contained in background domain .
4.4 Mesh iterators for overlapping meshes
Mesh iterators have been advocated by Berti [16, 17] and Logg [42] and used among others by Bastian et~al. [12] and Botsch et~al. [18] as an important abstraction concept in mesh implementations and FEM frameworks [43, 12]. The iterator concept allows to access and iterate over mesh entities such as vertices, edges, facets and cells without knowing the details of the underlying mesh implementation. We have followed the same ideas in the case of overlapping meshes. The two pairs of classes CutCell, CutCellIterator and CutFacetPart, CutFacetPartIterator provide an interface to the intersected cells in and the interface facet partition , respectively. The FacetPart class and the corresponding iterator class mimic the original interface in DOLFIN by giving access to the two adjacent cells in the overlapping and the overlapped meshes, respectively. The interface assembly in Figure 7 presents a important use case.
Similar iterator concepts have been used in Bastian et~al. [13] where a general infrastructure to couple grids interfaced by the DUNE grid framework is presented. The CutFacetPartIterator corresponds to the RemoteIntersectionIterator described in Bastian et~al. [13] specialized to our case of Nitsche’s method on overlapping meshes. Moreover, the CutCellIterator introduced here can be interpreted as an instance of DomainIntersectionIterator in Bastian et~al. [13].
4.5 The class Quadrature and QuadratureRuleCache
The DOLFIN class Quadrature is a lightweight base class which only computes and stores quadrature data such as quadrature points, weights and order for a given polyhedron at runtime. Since the integration order depends on the underlying finite element scheme, the actual computation of the points and weights is meant to be implemented in subclass constructors. The class BarycenterQuadrature is such an instance of a subclass which computes a quadrature rule of order for a given polyhedron based on the algorithm outlined in Section 3.4. For each intersected cell or facet part, a QuadratureRule object can be stored in a QuadratureRuleCache instance to save geometry and quadrature rule recomputation if several integrations have to be performed on the same intersected entity, as is the case for the computation of the stiffness matrix and load vector in the Nitsche overlapping mesh method.
4.6 Forms and assembly
The DOLFIN Form form class represents the mathematical concept of a finite element variational form. This class has been extended, as part of DOLFINOLM, to reflect the domain decomposition character of the overlapping mesh method. A socalled NitscheForm class holds the description of the variational problem on each part and of the domain. The coupling between the two forms is accomplished through a member function tabulate_interface_tensor, which computes the local interface tensor corresponding to (17). This is in addition to the standard tabulate_tensor functions defined in the UFC code generation interface [5] for cell integrals, exterior facet integrals, and interior facet integrals (cf. Algorithm 1). In addition, the NitscheForm class gives access to the two overlapping meshes as well as a quadrature cache to avoid recomputation of quadrature rules.
Degrees of freedom of the cells of that are entirely covered by the overlapping mesh are inactive; that is, those degrees of freedom do not determine the solution and can be assigned an arbitrary value. For practical reasons, inactive degrees of freedom are included in the linear system but are set to zero by inserting ’one’ on the diagonal of the corresponding rows and ’zero’ in the righthand side vector. This is automatically handled by the DOLFIN function Matrix::ident_zeros.
5 Numerical examples
To demonstrate the efficiency of our realization of the overlapping mesh method, we study two test examples. The first example is the Poisson equation. Here, an analytical solution for a suitable source function allows to verify the implementation by means of convergence studies. In addition, timings for the computation of both a standard finite element method and its corresponding Nitsche approximation are presented. A further breakdown of the timings in the Nitsche case gives a clear picture of the additional costs associated with the geometryrelated computations and their effect on the overall computation time. From that perspective, comparing Nitsche with a simple piecewise linear, continuous finite element method represents the most challenging test case, since then the extra work required for geometryrelated computations contributes the most to the overall assembly time.
As a second example, a linear elastic equation is solved with discontinuous material parameters at the interface between two overlapping meshes.
5.1 Poisson equation
We consider the elliptic model problem (1)–(5) on the domain . The overlapping domain is a translation and rotation of the cube according to Figure 9. The source function is given by
(25) 
Homogeneous Dirichlet boundary conditions are used on the entire boundary. The exact solution is given by
(26) 
The penalty parameter in (9) is set to . A code extract from the implementation of the solver based on DOLFINOLM is shown in Figure 8.
The Nitsche approximation was computed on a sequence of meshes with decreasing mesh size for the tessellations and , starting at and stopping at . For the integration over cut cells and interface facets, barycenter quadrature was employed. The resulting linear systems were solved by the preconditioned conjugate gradient method as implemented in PETSc in combination with the algebraic multigrid solver from Hypre used as a preconditioner. For both the standard finite element method and the overlapping mesh method, a mesh independent number of CG iterations was observed (3–4 and 911, respectively). All computations were carried out on a Macbook Pro equipped with a GHz Intel Core i7 processor and GB of RAM (1066 MHz DDR3). The benchmarks were repeated 10 times and averaged to obtain the reported results.
Convergence
As the theoretical results recalled in (11)–(12) predict, an optimal convergence rate in both in the  and norm is observed; see Figure 11. Figure 10 clearly illustrates the smooth transition from the solution part defined on the overlapped mesh to the part defined on the overlapping mesh .
Benchmarks
The comparison in Figure 12 shows that the Nitsche overlapping mesh method is only twice as expensive as the standard finite element method for the same mesh size and approximately the same number of degrees of freedom. A breakdown of the computing time shows that both the assembly and the solution of the linear system become twice as expensive when the overlapping mesh method is used. The latter can be attributed to the observed higher iteration numbers. The CPU time for the linear solve displays a kink in the slope at ca 3 and 4.5 million cells for the standard and Nitsche FEM, respectively, which may be attributed to the problem size increasing beyond a hardwarespecific threshold (cache or memory size) that induces overhead.
A further breakdown of the assembly time is shown at the bottom of Figure 12. We note that while the number of cells scales like on a regular grid, the number of facets scales like . Consequently, all purely interface related operations with linear complexity should scale like , where is the number of cells. This is also the case for the interface decomposition, computation of the cut cells, and integration over the interface and the cut cells. On the other hand, the initial collision detection between the overlapped and the overlapping mesh is somewhat more expensive as it involves the searching of treelike data structures (adding a logarithmic factor to the complexity). However, this extra cost is negligible compared to the total computing time. In conclusion, the timings depicted in Figure 12 indicate that the cost of assembly for the Nitsche overlapping mesh method is comparable to standard finite element assembly, and its relative efficiency increases with increasing mesh size.
A note on illconditioned stiffness matrices
As analyzed in Burman [19] and illustrated by numerical experiments in [47], the stiffness matrix stemming from Nitsche’s method on overlapping meshes and related schemes can be illconditioned if the intersection gives very small elements, compared to the original element size. Different approaches have been investigated to cure the schemes from illconditioned systems, either by choosing proper weights for the interface [38], replacing basis functions on small elements by extensions of basis functions on larger neighboring elements, or through the use of a socalled ghostpenalty [19, 20]. For the current work, the only precaution was to skip all cut cells with a relative measure of size smaller that . In our numerical experiments, we observed a meshindependent iteration number which indicates that illconditioning was not present. However, proper handling of small cells by introducing ghost penalties has been studied recently for the Stokes problem in Massing et~al. [48].
5.2 Linear elasticity
As a second example, we consider a linear elastic body occupying a domain in consisting of two subdomains and with possibly different material parameters. The displacement and the stresses are assumed to be continuous across the interface between the subdomains. The corresponding linear elasticity problem then takes the form: find the displacement such that
(27)  
(28)  
(29)  
(30)  
(31) 
Here, the stress tensor is related to the displacement vector by Hooke’s law
(32) 
where and are the Lamé parameters in for and is the strain tensor.
Nitsche’s method for the Poisson equation (1) proposed by Hansbo et~al. [33] can be adapted to the case of the linear elastic problem (27)–(31) and takes the following form: find such that
(33) 
where
(34)  
(35)  
(36) 
with a positive penalty parameter. Assuming that the elastic material is not nearly incompressible ( remains bounded), we may extend the analysis in Hansbo et~al. [33] and prove optimal order a priori error estimates. See Hansbo and Larson [35] and Becker et~al. [14] for details on discontinuous Galerkin methods for elasticity problems.
Test configuration and numerical results
We consider the linear elasticity problem (27)–(31) with which is overlapped by a propellerlike domain where ; see Figure 3. The righthand side is zero and the boundary conditions on are defined by
(37) 
where represents a combination of pure tangential, rotational force and a normal pressure:
(38) 
The Lamé parameters are given by in for with , , .
The numerical results are shown in Figures 13 and 14. The results indicate a “smooth” transition of the solution from the overlapped mesh to the overlapping mesh.
6 Conclusions and outlook
We have demonstrated that overlapping mesh methods, in particular the Nitsche overlapping mesh method, may be implemented efficiently in three space dimensions through the use of tree search data structures and tools from computational geometry. Numerical tests show that optimal order convergence is obtained, that the overhead of the overlapping mesh method compared to a standard finite element method is small (roughly factor two), and that the overhead is decreasing as the size of the mesh is increased.
In the future, we plan to extend our implementation and the techniques studied in this work to handle fluid–structure interaction problems as well as contact problems. Furthermore, we plan to fully automate the implementation of Nitsche formulations on overlapping meshes by adding code generation capabilities for interface terms to FEniCS.
Acknowledgments
We thank Harish Narayanan for valuable discussions regarding the elasticity test case. This work is supported by an Outstanding Young Investigator grant from the Research Council of Norway, NFR 180450. This work is also supported by a Center of Excellence grant from the Research Council of Norway to the Center for Biomedical Computing at Simula Research Laboratory.
References
 [1] Cgal, Computational Geometry Algorithms Library, software package. URL http://www.cgal.org.
 [2] gts, GNU Triangulated Surface Library, software package. URL http://gts.sourceforge.net/.
 AkenineMöller et al. [2008] T. AkenineMöller, E. Haines, and N. Hoffman. Realtime rendering. A.K. Peters Ltd, 2008.
 Alnæs [2012] Martin S. Alnæs. UFL: a Finite Element Form Language, chapter 17. Springer, 2012.
 Alnæs et al. [2009] Martin S. Alnæs, Anders Logg, KentAndre Mardal, Ola Skavhaug, and Hans Petter Langtangen. Unified framework for finite element assembly. Int. J. Comput. Sci. Eng., 4(4):231–244, 2009.
 Areias and Belytschko [2006] P. Areias and T. Belytschko. A comment on the article “A finite element method for simulation of strong and weak discontinuities in solid mechanics” by A. Hansbo and P. Hansbo. Commentary. Comput. Methods Appl. Mech. Engrg., 195(912):1275–1276, 2006.
 Baiges and Codina [2009] J. Baiges and R. Codina. The fixedmesh ALE approach applied to solid mechanics and fluidstructure interaction problems. International Journal for Numerical Methods in Engineering, 81:1529–1557, 2009.
 Balay et al. [1997] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkhäuser Press, 1997.
 Balay et al. [2010] Satish Balay, Jed Brown, Kris Buschelman, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang. PETSc users manual. Technical Report ANL95/11  Revision 3.1, Argonne National Laboratory, 2010.
 Balay et al. [2011] Satish Balay, Jed Brown, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang. PETSc Web page, 2011. URL {}{}}{http://www.mcs.anl.gov/petsc}{cmtt}.
 Bastian et~al. [2004] P.~Bastian, M.~Droske, C.~Engwer, R.~Klöfkorn, T.~Neubauer, M.~Ohlberger, and M.~Rumpf. Towards a unified framework for scientific computing. In Proc. of the 15th International Conference on Domain Decomposition Method, 2004.
 Bastian et~al. [2010a] P.~Bastian, F.~Heimann, and S.~Marnach. Generic implementation of finite element methods in the Distributed and Unified Numerics Environment (DUNE). Kybernetika, 46(2):294315, 2010a.
 Bastian et~al. [2010b] Peter Bastian, Gerrit Buse, and Oliver Sander. Infrastructure for the coupling of Dune grids. Numerical Mathematics and Advanced Applications 2009, pages 107114, 2010b.
 Becker et~al. [2009] Roland Becker, Erik Burman, and Peter Hansbo. A Nitsche extended finite element method for incompressible elasticity with discontinuous modulus of elasticity. Comput. Methods Appl. Mech. Engrg., 198(4144):33523360, 2009.
 Belytschko et~al. [2001] Ted Belytschko, Nicolas Moës, S.~Usui, and C.~Parimi. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering, 50(4):9931013, 2001.
 Berti [2002] Guntram Berti. Generic programming for mesh algorithms: Towards universally usable geometric components. In Proceedings of the Fifth World Congress on Computational Mechanics,Vienna University of Technology July, 2002.
 Berti [2006] Guntram Berti. GrAL  The Grid Algorithms Library. Future Generation Computer Systems, 22, 2006.
 Botsch et~al. [2002] M.~Botsch, S.~Steinberg, S.~Bischoff, and L.~Kobbelt. OpenMesh  a generic and efficient polygon mesh data structure, 2002.
 Burman [2010] E.~Burman. Ghost penalty. Comptes Rendus Mathematique, 348(2122):12171220, 2010.
 Burman and Hansbo [2012] E.~Burman and P.~Hansbo. Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. Appl. Numer. Math., 62(4), 2012.
 Chazelle and Palios [1990] Bernard Chazelle and Leonidas Palios. Triangulating a nonconvex polytope. Discrete & Computational Geometry, 5(1):505526, 1990.
 Chesshire and Henshaw [1990] G.~Chesshire and W.~D. Henshaw. Composite overlapping meshes for the solution of partial differential equations. Journal of Computational Physics, 90(1):164, 1990.
 Cools [1997] Ronald Cools. Constructing cubature formulae: the science behind the art. Acta Numerica, 6(1):154, January 1997.
 Daux et~al. [2000] Christophe Daux, Nicolas Moï¿½s, John Dolbow, Natarajan Sukumar, and Ted Belytschko. Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering, 48(12):17411760, 2000.
 Day and Bochev [2008] D.~Day and P.~Bochev. Analysis and computation of a leastsquares method for consistent mesh tying. J. Comput. Appl. Math., 218(1):2133, August 2008.
 Dhia and Guillaume Rateau [2005] Hashmi~Ben Dhia and Guillaume Rateau. The Arlequin method as a flexible engineering design tool. International Journal for Numerical Methods in Engineering, 62(11):14421462, 2005.
 Ericson [2005] C.~Ericson. Realtime collision detection. Morgan Kaufmann, 2005.
 Fournier and Montuno [1984] A.~Fournier and D.~Y. Montuno. Triangulating Simple Polygons and Equivalent Problems. ACM Transactions on Graphics, 3(2):153174, 1984.
 Fries [2008] T.~P. Fries. A corrected XFEM approximation without problems in blending elements. International Journal for Numerical Methods in Engineering, 75(5):503532, 2008.
 Fries [2009] T.~P. Fries. The intrinsic XFEM for twofluid flows. International Journal for Numerical Methods in Fluids, 60(4):437471, 2009.
 Fries and Belytschko [2010] T.~P. Fries and T.~Belytschko. The extended/generalized finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering, 84(April):253304, 2010.
 Hansbo and Hansbo [2002] A.~Hansbo and P.~Hansbo. An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. Comput. Methods Appl. Mech. Engrg., 191(4748):55375552, 2002.
 Hansbo et~al. [2003] A.~Hansbo, P.~Hansbo, and Mats~G. Larson. A Finite Element Method on Composite Grids based on Nitsche's Method. ESAIM, Math. Model. Num. Anal., 37(3):495514, 2003.
 Hansbo and Hansbo [2004] Anita Hansbo and Peter Hansbo. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput. Methods Appl. Mech. Engrg., 193(3335):35233540, 2004.
 Hansbo and Larson [2002] P.~Hansbo and M.~G. Larson. Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche's method. Comput. Methods Appl. Mech. Engrg., 191(1718):18951908, 2002.
 Heroux et~al. [2005] Michael~A. Heroux, Roscoe~A. Bartlett, Vicki~E. Howle, Robert~J. Hoekstra, Jonathan~J. Hu, Tamara~G. Kolda, Richard~B. Lehoucq, Kevin~R. Long, Roger~P. Pawlowski, Eric~T. Phipps, Andrew~G. Salinger, Heidi~K. Thornquist, Ray~S. Tuminaro, James~M. Willenbring, Alan Williams, and Kendall~S. Stanley. An overview of the trilinos project. ACM Trans. Math. Softw., 31(3):397423, 2005.
 Jirásek and Belytschko [2002] M.~Jirásek and T.~Belytschko. Computational resolution of strong discontinuities. In Fifth world congress on computational mechanics, pages 712, 2002.
 Johansson and Larson [2012] August Johansson and Mats~G. Larson. A high order discontinuous Galerkin Nitsche method for elliptic problems with fictitious boundary. submitted to Numerische Mathematik, 2012.
 Kirby and Logg [2006] Robert~C. Kirby and Anders Logg. A Compiler for Variational Forms. ACM Trans. Math. Softw., 32(3):417444, 2006.
 Lee and Requicha [1982] Y.~T. Lee and A.~A.~G. Requicha. Algorithms for computing the volume and other integral properties of solids. I. known methods and open issues. Communications of the ACM, 25(9):635641, 1982.
 Logg [2007] Anders Logg. Automating the finite element method. Arch. Comput. Methods Eng., 14(2):93138, 2007.
 Logg [2009] Anders Logg. Efficient representation of computational meshes. Int. J. of Computat. Sci. Eng., 4(4):283295, 2009.
 Logg and Wells [2010] Anders Logg and Garth~N. Wells. DOLFIN: Automated finite element computing. ACM Trans. Math. Softw., 37(2), 2010.
 Logg et~al. [2012a] Anders Logg, KentAndre Mardal, Garth~N. Wells, et~al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012a.
 Logg et~al. [2012b] Anders Logg, Kristian~B. Ølgaard, Marie~E. Rognes, and Garth~N. Wells. FFC: the FEniCS Form Compiler, chapter~11. Springer, 2012b.
 Logg et~al. [2012c] Anders Logg, Garth~N. Wells, and Johan Hake. DOLFIN: a C++/Python finite element library, chapter~10. Springer, 2012c.
 Massing [2012] A.~Massing. Analysis and implementation of finite element methods on overlapping and fictitious domains. PhD thesis, University of Oslo, 2012.
 Massing et~al. [2012] A.~Massing, Mats~G. Larson, A.~Logg, and Marie~E. Rognes. A stabilized Nitsche overlapping mesh method for the Stokes problem. submitted, 2012.
 Mayer et~al. [2009] U.~M. Mayer, A.~Gerstenberger, and W.~A. Wall. Interface handling for threedimensional higherorder XFEMcomputations in fluidstructure interaction. International Journal for Numerical Methods in Engineering, 79(7):846869, 2009.
 Mayer et~al. [2010] U.~M. Mayer, A.~Popp, A.~Gerstenberger, and W.~A. Wall. 3D fluidstructurecontact interaction based on a combined XFEM FSI and dual mortar contact approach. Comput. Mech., 46(1):5367, 2010.
 Mirtich [1996] Brian Mirtich. Fast and accurate computation of polyhedral mass properties. journal of graphics tools, 1(2):3150, 1996.
 Moës et~al. [1999] N.~Moës, J.~Dolbow, and T.~Belytschko. A finite element method for crack growth without remeshing. Int. J. Numer. Meth. Engng, 46:131150, 1999.
 Natarajan et~al. [2009] S.~Natarajan, S.~Bordas, and D.~R. Mahapatra. Numerical integration over arbitrary polygonal domains based on SchwarzChristoffel conformal mapping. International Journal for Numerical Methods in Engineering, 80(1):103134, 2009.
 Natarajan et~al. [2010] Sundararajan Natarajan, D.~Roy Mahapatra, and Stéphane P.~A. Bordas. Integrating strong and weak discontinuities without integration subcells and example applications in an XFEM/GFEM framework. International Journal for Numerical Methods in Engineering, pages 269294, July 2010.
 Ølgaard et~al. [2008] Kristian~B. Ølgaard, Anders Logg, and Garth~N. Wells. Automated code generation for discontinuous galerkin methods. SIAM Journal on Scientific Computing, 31(2):849864, 2008.
 Ruppert and Seidel [1989] J.~Ruppert and R.~Seidel. On the difficulty of tetrahedralizing 3dimensional nonconvex polyhedra. In Proceedings of the fifth annual symposium on Computational geometry  SCG '89, pages 380392, New York, New York, USA, June 1989. ACM Press.
 Schneider and Eberly [2003] Philip~J. Schneider and David~H. Eberly. Geometric tools for computer graphics. Morgan Kaufmann, 2003.
 Sukumar et~al. [2000] N.~Sukumar, N.~Moës, B.~Moran, and T.~Belytschko. Extended finite element method for threedimensional crack modelling. International Journal for Numerical Methods in Engineering, 48(11):15491570, 2000.
 Yang and Laursen [2008] B.~Yang and T.~A. Laursen. A contact searching algorithm including bounding volume trees applied to finite sliding mortar formulations. Comput. Mech., 41(2):189205, 2008.
 Yu [2005] Z.~Yu. A DLM/FD method for fluid/flexiblebody interactions. Journal of Computational Physics, 207(1):127, 2005.
 Zhang et~al. [2004] L.~Zhang, A.~Gerstenberger, X.~Wang, and W.~K. Liu. Immersed finite element method. Comput. Methods Appl. Mech. Engrg., 193(2122):20512067, 2004.