CutFEM: Discretizing geometry and partial differential equations

We discuss recent advances on robust unfitted finite element methods on cut meshes. These methods are designed to facilitate computations on complex geometries obtained, for example, from computer‐aided design or image data from applied sciences. Both the treatment of boundaries and interfaces and the discretization of PDEs on surfaces are discussed and illustrated numerically. Copyright © 2014 John Wiley & Sons, Ltd.


INTRODUCTION
Research in numerical methods for solving problems in science and engineering has mainly focused on techniques for approximating models described by partial differential equations (PDEs), while the important coupling to the geometrical description of the domain has been largely overlooked. In recent years, the need for a unified approach has been recognized, and this area is today receiving rapidly increasing interest. This interest is motivated by the demand for efficient and robust techniques for solving problems involving complex and evolving geometries. The use of geometric descriptions in the computational method, that are closely linked to the geometric data acquisition, can dramatically reduce the computational cost of preprocessing, or transformation, of acquired geometry descriptions into representations suitable for the computational method at hand.
For instance, the simulation of blood flow dynamics in vessel geometries requires a series of highly non-trivial steps to generate a high quality, full 3D finite element mesh from biomedical image data [1]. Similar challenging and computationally costly preprocessing steps are required to transform geological image data into conforming domain discretizations, which respect complex structures such as faults and large scale networks of fractures; see, for instance, [2].
An example of a successful paradigm for the integration of geometry and computation is given by the isogeometric analysis pioneered by Hughes and co-workers [3]. Here, the merging of computation and geometry is obtained by adopting the functions used for geometry representation as basis for the computational method leading to new approaches in the discretization of PDEs.
The idea behind CutFEM is to make the discretization as independent as possible of the geometric description and minimize the complexity of mesh generation, while retaining the accuracy and 473 robustness of a standard finite element method. In particular, we will show later how recent stabilization techniques can be applied to make both the accuracy of the approximation and the system condition number independent of the mesh/boundary intersection and physical parameters. Thanks to this robustness of the discretization, powerful linear algebra techniques developed for finite element methods can be made to bear on the solution of the linear systems obtained by the CutFEM discretization.
In the CutFEM approach, the boundary of a given domain is represented on a background grid, for instance, using a level set function. The background grid is then also used to represent the approximate solution of the governing PDEs. CutFEM builds on a general finite element formulation for the approximation of PDEs, in the bulk and on surfaces, that can handle elements of complex shape and where boundary and interface conditions are built into the discrete formulation. This way, CutFEM can ease the burden of mesh generation by requiring only a low-quality and even non-conform surface mesh representation of the computational geometry. The integration of the geometry in the discrete formulation leads to a method that can be applied equally well to CAD generated geometries and to geometries obtained from biomedical or geological image data.
In this paper, we give some examples of how CutFEM combined with Nitsche's method [4] is implemented for a range of problems with increasing complexity. The use of Nitsche's method for unfitted interface problems and fictitious domain methods has been developed in [5-13, 28, 29, 48]. Other related approaches based on Lagrange multipliers or discontinuous Galerkin methods have been suggested in the following works [14][15][16][17][18][19][20].
The paper is organized as follows. In Section 2, we give an overview of some archetypal problems with associated CutFEM discretizations. Then we will discuss the crucial question of robustness in Section 3. The representation of the geometry using level sets is discussed in Section 4 and implementation issues are reviewed in Section 5. A series of numerical illustrations for different coupling problems on non-trivial geometries are presented in Section 6.

A Poisson model problem
Let be a bounded domain in two or three space dimensions, with an interface dividing into two non-overlapping subdomains 1 and 2 , so that WD 1 [ 2 [ , with the interface defined by D 1 \ 2 . For simplicity, we assume that the subdomains are polyhedral (or polygonal in R 2 ) and that is polygonal (or a broken line). For any sufficiently regular function u in 1 [ 2 , we define the jump of u on by u WD u 1 j u 2 j , where u i D uj 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 . For definiteness, we define n as the outward pointing unit normal to 1 on . We shall consider a problem with piecewise constant data˛so that˛i D˛j i , and a reaction coefficient c.x/ > 0 that may be discontinuous across . Our first model problem is the following variant of Poisson's equation: where we used the flux vectors q i WD ˛i ru i .
In a standard finite element method, the jump in the normal derivative resulting from the continuity of the flux q WD ˛ru, when˛1 ¤˛2, can be taken into account by letting coincide with mesh lines. In [5,7], another approach was taken in that (1) was solved approximately using piecewise polynomial finite elements on a family of conforming shape regular triangulations T h of that were independent of the location of the interface . The approximation was then allowed to be discontinuous inside elements, which intersected the interface.
To recall this method, we will use the following notation for mesh related quantities. For any element K, let K i D K \ i denote the part of K in i . By G h WD ¹K 2 T h W K \ ¤ ;º, we here denote the set of elements that are intersected by the interface. For an element K 2 G h , let K WD \ K be the part of in K. We shall seek a discrete solution U D .
As may intersect two edges of a triangle arbitrarily, the size of the parts K i is not fully characterized by the meshsize parameters. Thus, to guarantee stability of this method using elements with internal discontinuities, further conditions on the combinations of numerical fluxes (co-normal derivatives) must be imposed by choosing appropriate mesh and geometry dependent weights Ä.
The approach suggested in [5,7] was to choose the numerical fluxes by introducing weights where meas.K/ denotes the size (area or volume) of K, and to use a weighted mean where @ n WD n r, as the flux on (note that Ä 1 C Ä 2 D 1). Thus, for an intersected element, a mean numerical flux that takes the different meshsizes into account was computed. However, for problems in which there is a very large difference between the parameters˛i , this approach is not robust. To improve robustness, the weighted average must also take into account the parameters: in (3), we must change the definitions of the Ä i to cf. [21,22]. The method is defined by the variational problem of finding U 2 V h such that where For stability, one must take sufficiently large. Two choices of have been proposed in the literature for these weights. In [21], it was suggested to take 475 and in [22], the choice was advocated. Above O is a mesh and parameter independent constant and With these definitions, the discrete problem (5) is consistent in the sense that, for u solving (1), In Section 3, we will discuss the consequences of the different choices of the penalty parameter and an alternative route to robustness inspired by fictitious domain methods. A FE basis for V h is easily obtained from a standard FE basis on the mesh by the introduction of new basis functions for the elements that intersect . Thus, we replace each standard basis function living on an element that intersects the interface by two new basis functions, namely its restrictions to 1 and 2 , respectively. The collection of basis functions with support in i is then clearly a basis for V h i , and hence, we obtain a basis for V h by the identification D . j 1 ; j 2 /. If the interface coincides exactly with an element edge, no new basis functions are introduced on these elements, but the approximating functions may still be discontinuous over such an edge. As a consequence, there are six non-zero basis functions on each element that properly intersects . Perhaps this process is most easily seen as creating two new separate meshes with doubling of the elements crossed by the interface (Figure 1).
We now want to show that the approximation property of V h is optimal in the following mesh dependent norm: Let I h be the standard nodal interpolation operator and define We then have the following result. Let I h be an interpolation operator defined as in (10).
In the proof of this result, we need to estimate the interpolation error at the interface. To that end, the following trace inequality is necessary: under reasonable mesh assumptions [5][6][7], there exists a constant C T , depending on but independent of the mesh, such that The crucial fact is that the constant in this inequality is independent of the location of the interface relative to the mesh. Optimal interpolation estimates follow, as does optimal convergence of the method irrespective of the location of the interface relative to the mesh. The key elements of the analysis are the robust coercivity of a h . ; / with respect to the norm h , the consistency property (8), and the approximability (11). For details, see [5].

The case of Dirichlet boundaries
We now consider the boundary value problem where D D [ N denotes the boundary of the domain , discretized on a mesh T h that contains but is not fitted to the domain boundary. Denoting the mesh domain T , we consider the following finite element space: By G h WD ¹K 2 T h W K \ ¤ ;º, we denote the set of elements that are intersected by the interface. For an element K 2 G h , let K WD \K be the part of in K. We also introduce the set of element faces F G associated with G h , defined as follows: for each face F 2 F G , there exist two simplices K and K 0 such that F D K \K 0 and at least one of the two is a member of G h . This means in particular that the boundary faces of the mesh T h are excluded from F G . On a face F such that F D K \ K 0 , define the jump of the gradient of v 2 W h by @ n F v D n F rv j K n F rv j K 0 , where n F denotes the unit normal to F , the orientation is chosen arbitrarily. The problem that arises when applying Nitsche's method in this framework is that the inverse inequality corresponding to (12) cannot be robust when the right hand side must be controlled by norms over the physical domain alone, as the cut elements can be arbitrarily small. We thus need to further stabilize the problem. The finite element discretization takes the form: find U 2 W h , such that where and Here, D ; N , and 1 are positive penalty parameters; cf. [12]. The rationale for the stabilization term (16) is that it extends the coercivity from the physical domain to the mesh domain T . Details are given in Section 3.

Other interface conditions of interest
The interface conditions may vary depending on the application and may pose more discretization challenges than the aforementioned model problems. In this section, we mention a few alternatives that can still be handled in the same framework of Nitsche's method.
Heat release at the interface leads to u D 0 on ; q n D g on ; where g is a heat source term. This source leads to a modified right-hand side in (8) so that see [5]. Spring-type boundary conditions common in solid mechanics can be modeled by u D k q n on ; q n D 0 on : Here, k denotes the compliance of distributed springs on the interface. These conditions can be modeled by Nitsche's approach as follows. Let s h D 1=.h= C k/ and modify the bilinear form to The analysis of this method follows the same lines as the one for the original method; cf. [7,23]. Transport also on the interface can be modeled by where r is the tangential derivative on , r WD r n @ n : These conditions model, for example, porous media flow in a medium with a crack represented by ; cf. [24]. Here,˛ represents a porosity along the crack. The bilinear form must now be modified to take the differential equation on into account: formally replacing g in (17) by r .˛ r u/, we see that a consistent bilinear form is given by where ¹aº WD k 2 a 1 C k 1 a 2 , with k 1 and k 2 positive weights. This method requires additional stabilization in general; cf. the numerical example in Section 6.4. Alternative surface transport conditions are given by seeking u W ! R and u W ! R, where denotes the boundary of , such thať u ˇ u C˛@ n u D 0 on ; @ n u r .˛ r u / D g on : Here, u is a concentration on the surface, which is independent of the concentration inside . Applications are found, for example, in cell membrane transport; cf. [25]. Here, we no longer have a distinct side condition and can dispense with Nitsche's method. However, with cut elements, we now need a way to define the discrete approximation of u , which is different from the previous case, where the trace spaces on the cut elements were used to compute r U . An obvious idea is to keep on using the trace spaces on a higher dimensional mesh as suggested by Olshanskii et al. [26]. Thus, we let W h denote the space of continuous piecewise linear polynomials defined on G h and seek U 2 W h and U 2 W h such that and .
or the symmetrized versioň The discretization of the equation on the interface can now be stabilized by adding related to the stabilization in the Dirichlet case but with another scaling because of dimensionality. The bulk variable can be stabilized as in the Dirichlet case.

ENHANCING ROBUSTNESS: GHOST PENALTY
Cutting the mesh can result in boundary elements with very small intersections with the physical domain, or for PDEs on embedded surfaces, bulk elements with very small intersection with the surface domain. This may lead to a poorly conditioned system matrix or failure of stability of the discrete scheme.
Situations that are particularly sensitive are the imposition of Dirichlet boundary conditions [12,27], or domain decomposition on unfitted meshes, where an inf-sup condition has to be satisfied, such as for incompressible elasticity [9,[28][29][30][31]. In these cases, one cannot choose the weights in (3) in a robust way. There are also situations where the weights are already prescribed by other concerns. This is the case in fluid-structure interaction [32,33], where the elastodynamic system has no dissipation by which one can absorb the contribution of the boundary stress term and therefore only the fluid stresses are considered on the boundary. If independent adaptive mesh refinement is performed in the two subdomains of (5), this also imposes a certain choice of weights to ensure robustness, both with respect to large contrast in the physical parameter and in the mesh parameter [34].
For cases such as this, a useful trick [12,27] is to add a penalty term in the interface zone that extends the coercivity to the whole mesh domain, that is, in the O.h/ zone of the mesh domain of each subdomain that does not intersect the associated physical domain. This penalty term must be carefully designed to add sufficient stability, while remaining weakly consistent for smooth solutions. To illustrate this idea, we consider the fictitious domain method for the Poisson problem (13) with˛D 1. In the next section, we demonstrate how these arguments can be extended to the coupled problem (1) and how the two formulations are related.
We observe that by taking v D U in the bilinear form a.U; v/, we have the coercivity krU k 2 L 2 . / 6 a.U; U /: However, to obtain coercivity of the form a h .U; v/ using this stability and the boundary penalty term, the penalty parameter will depend on the cut, as by the Cauchy-Schwarz inequality By the definition of the norm and because rU is piecewise constant, we have, using (12), where we recall the notation K D K \ , K D K \ , m K WD meas d .K / and m K WD Unfortunately, this makes strongly dependent on the cut, because for m K D O.h K /, the volume measure m K can be arbitrarily small. When the penalty parameter becomes strongly dependent on the mesh/boundary intersection, one may encounter problems with both conditioning and accuracy. A solution to this problem is to add the ghost penalty term of (16), denoted by j. ; /, to the form a h . ; / as in the formulation (14). The role of this term is to extend the coercivity from the physical domain to the mesh domain T WD [ G h . Indeed, one may show the following inequality: where c G > 0 is bounded away from zero independent of the mesh/boundary intersection for positive ghost penalty stabilization parameter 1 . The following weak consistency property can also be shown to hold: where the constant C is independent of the mesh/boundary intersection. Coercivity now follows from (24) and (25) as follows: Here, C T is the constant of the trace inequality (12) and c G is the coercivity constant of the stability estimate (25). We conclude by choosing > 2C 2 T c 1 G , where the lower bound is independent of the mesh/boundary intersection, but not of the penalty parameter 1 in j. ; /. Error estimates now follow in a similar fashion as for the standard Nitsche's method, using (27) and (26). Indeed, provided the solution is smooth enough, there holds where the hidden constant is independent of the mesh/boundary intersection. One may also show that the conditioning of the system matrix is bounded independently of the mesh/boundary intersection. For further details, see [12].

Example: perfect conductor, the limit of infinite diffusion
As an example of how the improved robustness works, we will consider the problem (1), with 2 completely enclosed by 1 such that @ @ 1 and WD @ 2 ( Figure 2). It is well known that in the limit˛2 ! 1, the solution u 2 becomes constant in 2 and can therefore be exactly represented by one degree of freedom. We will give two formulations of this problem, one similar to (5), and the other using the reduced model in a framework similar to the fictitious domain approach (14) using only one degree of freedom to represent u 2 . In both formulations, we will use the ghost penalty method so that the weights may be depending only on the diffusivities. This allows us to show that the solution of the domain decomposition approach converges to that of the fictitious domain approach in the limit as˛2 ! 1. We will then compare this to what is obtained if the weights are chosen as in equation (4), with the penalty defined in equation (6) or (7).
First, consider the formulation [35], and U 2 V h such that with a h .U; v/ and L.v/ as defined in (5), but for simplicity with f 2 Á 0 and using the weights and the penalty parameter This choice of weights was first considered in the context of discontinuous Galerkin methods in [36] and then for Nitsche's method in [37]. The ghost penalty terms are designed similarly to the formulation (14). However, here, j i .
It is straightforward to prove coercivity of this formulation using the arguments of the previous section, indeed, Here, we used the inequality and the stability (25) in both subdomains. If we formally let˛2 ! 1 and replace U 2 and v 2 by constant functions, we find the following formulation: find Q where D 4˛1C 2 T c 1 G . It follows by inspection that this formulation is equivalent to (14), with the only difference that in this case the Dirichlet value set on is an unknown constant value. The coercivity and hence the discrete wellposedness can be shown here exactly using the arguments of (27) together with a triangle inequality and a trace inequality to get control of Q U 2 . First note that by a triangular inequality, a trace inequality on the whole domain 1 , and a Poincaré inequality, we have Therefore, by the same arguments as above, there exists a constant c > 0 such that Consequently, both (28) and (32) are coercive, and with some additional effort, one may prove that the approximations indeed have optimal convergence in the H 1 and L 2 -norm provided the exact solution is sufficiently smooth when restricted to each subdomain.
Here, we will instead consider the robustness of the formulation (28). A natural question to ask, as (32) was obtained by taking the formal limit˛2 ! 1 in (28), is if the solutions of the two formulations also coincide in the limit. We will show now that this is the case. Let U D .U 1 ; U 2 / denote the solution of (28) and Q U D . Q U 1 ; Q U 2 / denote the solution of (32). If we set E D Q U U and note that E 2 V h , we may use the coercivity proven in (31): Using the formulation (28) and once again that E 2 V h , we see that Now, we apply the formulation (32) on the form Using (33) and c˛1 6 6 C˛1 when˛1 6˛2, it follows that and we conclude that E ! 0 as˛2 ! 1 meaning that the asymptotic limit of the solution of (28) coincides with that of (32). Another important observation is that for the discretization using ghost penalty and weights depending only on the diffusion, preconditioning the system matrix using diagonal scaling with 1 ;˛2 leads to a system whose condition number is independent of both the mesh/boundary intersection and the contrast in the diffusion (for details, see [35]).
Note that the use of the weights (4) and the penalty parameters (6) or (7) do not allow a similar robust limit formulation. This can be seen in the following way. Consider the penalty parameter for the choice (6) proposed in [21]. Then the limit satisfies which is robust in the diffusion parameter but degenerates into the choice of weights for the fictitious domain method that depends on the element cuts. Hence, for the weights (6), preconditioning using diagonal scaling will make the system robust for large contrast, but unfortunate cuts will still result in a poorly conditioned system matrix. This disadvantage motivated the introduction of the ghost penalty operator in the previous section. The choice (7) on the other hand results in the limit behavior, This implies that U D 0 across the boundary, which may not always be compatible with optimal accuracy. Observe also that the penalty parameter is present in the system matrices corresponding to the bulk discretization in both subdomains. As becomes big, with increasing˛2, it will destroy the conditioning of the matrix corresponding to subdomain 1 .

A numerical illustration
Robustness issues may also appear in the limit of vanishing diffusion in a subdomain. We consider a configuration similar to that of Figure 2, with the square domain WD OE 1; 1 2 . Let 2 be a disc with radius 0.75 centered in the origin and 1 WD n 2 . In this configuration, we solve (1), with uj @ D 0,˛1 D 1, f D 1, cj 1 D 0, cj 2 D 1 and with˛2 2 ¹10 4 ; 10 6 ; 10 9 º. First, we apply the formulation (5) with weights given by (3). The elevations of the solutions are presented in Figure 3. We then solve the same problem using the method (28) with (29)-(30) and give the corresponding elevations in Figure 4. Observe the relatively strong, but local, spurious overshoots that are present close to the layer in    the parameters (29)-(30) eliminates the oscillations by relaxing the continuity constraint in the limit. Indeed, the sharp layer is represented as a discontinuity when it is under resolved. Now, we consider the case where cj 1 D cj 2 D 0,˛2 D 1 and choose˛1 D 10 20 , to illustrate the effect in the limit˛1 ! 1. This corresponds to a situation similar to that explored in Example 3.1, but with the diffusivity going to infinity in the outer domain. In Figure 5, we compare the results of the method (5) with the weights (3) and of (28) with (29)- (30). We see that for this large contrast, the method using the weights (3) exhibits some spurious oscillations close to the boundary, probably owing to an incompatibility between the constraint U D 0 and the weakly imposed condition on the gradient. In other words, the finite element space defined on the cut mesh does not allow for a H 1 -conforming interpolant that also can represent the jump in the gradient. This effect is reminiscent of locking but only present in the vicinity of the interface. The method (28) on the other hand converges to the fictitious domain solution of (14) in the limit for which optimal error estimates exist.

Ghost penalty for surface partial differential equations
Let us consider the Laplace-Beltrami problem: find u W ! R such that where D r r is the Laplace-Beltrami operator. The finite element method takes the form: and we recall that W h is the space of continuous piecewise linears on G h , the union of elements intersecting .
In this case, the ghost penalty term provides additional stability to the discrete problem, which may be arbitrary ill conditioned. The proof of the optimal scaling estimate of the condition number, that is, basically relies on the Poincaré estimate and the inverse estimate Note that, in these inequalities, the ghost penalty term plays a crucial role. Using the Poincaré and inverse estimate together with the coercivity of the bilinear form and the standard scaled equivalence, between the Euclidian norm on the nodal values O v 2 R N of v 2 W h and the L 2 -norm on G h , we may prove (37); see [38] for details.
Furthermore, we note that the ghost penalty term scales in the same way as the bilinear form .r v; r w/ and that we have the interpolation error estimate Using the Galerkin orthogonality, we obtain an optimal error estimate in the energy norm and by duality an L 2 -error estimate, In order to deal with the effect of approximating the boundary , a more elaborate analysis is needed, and we refer to [38] for further details. For a numerical example of the solution of the Laplace-Beltrami equation on a non-trivial surface, see Section 6.3.

DISCRETIZING GEOMETRY IN CutFEM
The starting point for the discretization of the geometry in CutFEM is to immerse an arbitrary geometric description in a background mesh. This mesh is typically chosen structured, to facilitate the handling of data structures, communication, and hierarchic mesh adaption. The discretization space and the variational formulation are then adapted to the geometry so that suitable boundary 486 E. BURMAN ET AL. and interface conditions are imposed weakly as described in Section 2. This approach leads to some challenging implementation issues that we will describe later. It is advantageous to choose one particular geometry description that is versatile and simple for the construction of the discretization. Other geometrical descriptions can then be either included using modules that translate different formats to the one chosen to interface with the code or provided with their own CutFEM modules. In our work, later, we focus on the level-set method [39] for its use both as a description of stationary boundaries and of interfaces evolved by computation. Here, the location of the boundary is given by the zero level set of a function W R d ! R. More precisely, a given domain R d can be decomposed into an inner part 1 , an outer part c 1 , and their common interface by requiring that 8 x 2 8 < : where we define c 1 as the (open) complement of 1 in . To give a few non-trivial examples of analytically given level-set based surface descriptions, we introduce In the following, we choose R D 1:2, r D 0:3.

Swiss Cheese Block
The corresponding surfaces are depicted in Figure 6. Using a level-set description, complex domains can easily be constructed by translating Boolean set operations and geometric transformations into simple manipulations of level-set representations. For instance, given two level set functions 1 and 2 representing the domains 1 and 2 , respectively, the level-set function representing the result of a standard Boolean set operation can be constructed according to the following table:    Figure 7 illustrates how complex surfaces can be generated by combining these operations.

IMPLEMENTATIONAL ASPECTS OF A LEVEL-SET BASED CutFEM METHOD
We now briefly review some of the data structures and algorithms required to efficiently compute a discrete representation of surfaces implicitly given by a level-set and how to evaluate the variational formulation on this discrete geometry. An important step here is the introduction of a subtriangulation of cut elements to facilitate integration, which we will describe in more detail below. This subtriangulation is used only for integration so that the resulting subtriangles are not constrained by the conformity requirements of the finite element space. In addition, their aspect ratio does not impact the approximation properties of the discretization, because they are never used in the analysis.

E. BURMAN ET AL.
In general, the mesh subdivision can be characterized as follows. Referring to the notation from Section 2, we recall that for a tessellation T h of , the (unfitted) surface leads to natural subdivi- Note that, to apply the finite element method to a pure surface problem, we only need to compute a tessellation of with respect to G h . For a fictitious domain problem, we require a subtesselation for the inner part of the cut elements G h;1 D ¹K \ 1 W K 2 G h º and for an interface problem, we need to subtriangulate both the inner and outer parts, that is, also G h;2 D ¹K \ 2 W K 2 G h º. However, as the following section will explain, all these sub-triangulations can be provided simultaneously by using the well-known marching tetrahedra algorithm [39,41].
In addition to this subtesselation algorithm, routines for the computation of integrals over cut cell parts have to be provided. Indeed, for each intersected element K 2 G h , volume integrals have to be evaluated over the parts of K that are covered by the two physical subdomains 1 , K 1 D K \ 1 and 2 , K 2 D K \ 2 . To construct the matrix contributions that impose interface and boundary conditions, surface integrals over interface segments, K D \ K, that lie within the element K 2 G h , also have to be computed.
In the following, we will first detail our implementation of the classical tetrahedra algorithm [39,41] and then give some further details on how to evaluate integrals over cell parts.

Interface approximation and sub-triangulation of cut elements
In the first step of the subtesselation algorithm, the values of the level-set function in the element nodes are used to distinguish between elements that are fully contained in 1 ( < 0 for all nodes) or 2 ( > 0 for all nodes) and the elements that are intersected by the interface ( > 0 and < 0 in some nodes) (Figures 8 and 9). For elements that are intersected by the interface, we perform a subtriangulation of the element allowing us to apply standard quadrature rules to the subtriangles for the integration over the parts K 1 , K 2 , and K . Here, we only consider linear intersections of the zero level-set with the elements  (3) Build a subtriangulation of each cut cell part, that is, K 1 , K 2 and K . The subtriangulation of K i , i D 1; 2, and the sub-triangulation for K are then added to the subtesselation for the inner part of the cut elements G h;1 , the outer part of the cut elements G h;2 , or the tessellation of , respectively (Figure 11), together with a map between the cell parts and the corresponding parent cell. In two space dimensions, the cut cell parts, K 1 and K 2 , are either triangular or quadrilateral ( Figure 10). The triangular part can be added directly to the subtriangulation of the corresponding domain. The quadrilateral part is subdivided into two triangles first and then added to the corresponding subtriangulation. The interface segment K is represented by straight line segment connecting both intersection points. These interface line segments are stored in a separate mesh object with topological dimension 1. The resulting sub-triangulations for 1 and 2 for a circular domain 1 are shown in Figure 11. In three space dimensions, there is a much wider range of cut cases to consider. We can distinguish 14 cases, among which eight have a triangular interface plane cut of the zero level set with the tetrahedron and six with a quadrilateral interface plane cut ( Figure 12). The triangular plane cut can be added directly to the subtriangulation of the two-dimensional   zero level set interface representation, and the quadrilateral planes are subdivided into two triangular parts. For the volume subtriangulation, we decompose the tetrahedron into eight subtetrahedra ( Figure 13) and add the subtetrahedra to the corresponding subtesselations of 1 and 2 depending on the cell flag, that is, depending on whether they lie in 1 or 2 .

Integration over subtriangulation
For the evaluation of integrals over the subtriangulation, we require two mappings 1 w ı p . Here, the linear affine mapping, p , between the reference element and the cell part, transforms a quadrature rule defined on the reference element in terms of quadrature points i and weights w i into a quadrature rule on the subtriangle cell part (Figure 14). The mapping w between the reference element and the whole 'parent' cell is used to map the quadrature points defined on the physical subtriangle to their location in the reference element of the whole parent cell.
More precisely, the quadrature rule over the cell part is given in terms of quadrature points in physical space x p i D p . i / and quadrature weights w p i D w i jdet.J p /j. Here, J p is the Jacobian of the mapping p . These points and weights can then be used to perform the numerical quadrature over the given cell part, for example, for a mass matrix, by Figure 14. Schematics of integration over the subtriangulation involving two coordinate transformations p between the reference element and the sub-triangle cell part and w between the reference element and the whole cell.
Here, A T OEj OEk denotes the matrix entry for the cell integral over the degrees of freedom j and k.
Note here that we have to evaluate basis functions and/or their derivatives at arbitrary points on the reference element i D 1 w .x i / that result from the backward transformation of the quadrature rule over the cell part onto the reference element of the whole 'parent' element.

Software for CutFEM-type method
The technology required for the automated assembly of general cut finite element based variational forms over fictitious domains and embedded surfaces has been implemented as part of the software library libCutFEM. libCutFEM is open source and freely available from http://www.cutfem.org. This software library builds on the finite element library DOLFIN, which is part of the FEniCS project [42] 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 high-level descriptions of finite element variational problems expressed in near-mathematical notation [43].
libCutFEM makes use of this automatization technology in the following way. We specify the variational problem in near-mathematical notation in the domain specific Uniform Form Language [43] (for an example, see Figure 17). Then the code generation components of FEniCS (UFL+FFC+UFC [42]) generate code containing information about the cell integrals of the given forms and information about the elements such as the degrees of freedom maps. In particular, libCut-FEM uses an extension developed in earlier works of the components FFC [44,45] and UFC [46]. In these works, FFC and UFC were extended [33,47,48] to provide generated code for the cell integration over cut cell parts (see right box in Figure 14). These extensions provide the fundamental infrastructure for the evaluation of cut finite element based variational forms defined on fictitious domains and embedded surfaces. Consequently, given a high-level description of the variational formulation, low-level C++ code can be automatically generated for the evaluation of the cut element and surface integrals, in addition to the evaluation of integrals over the standard (non-cut) mesh entities. The generated code takes as input appropriate quadrature points and weights for each cut entity, which are provided by the libCutFEM library. The resulting cell integral matrix contributions are then assembled to the global matrix by routines provided in libCutFEM.
In addition to these quadrature and assembly routines, libCutFEM provides functionality for computing topological and geometric descriptions of the embedded surface and cut elements. Currently, libCutFEM mainly supports the level set based cutting algorithm as described in detail in the previous section. However, alternative geometrical algorithms can be integrated easily thanks to the modular structure of the library and FEniCS, which decouples the variational formulation, the quadrature and assembly routines, and the geometrical and mesh description. In summary, in libCutFEM, one may specify variational forms defined over finite element spaces on fictitious domains in high-level UFL notation [43], define the background mesh T h , and give a description of the surface , and then invoke the functionality provided by the libCutFEM library to automatically assemble the corresponding stiffness matrix. In particular, the numerical experiments presented in Sections 6.1-6.3, corresponding to the variational formulation defined by (50), (58) and (60), have been carried out using this technology. We also present the UFL scripts that were used to obtain the numerical results illustrating the simplicity by which the end user can access the CutFEM paradigm (Figures 17, 18, and 21). Other software projects including CutFEM style technology include GetFem++ [49], LifeV [50], and DUNE-UDG [51].

Stabilized Nitsche fictitious domain method for the Poisson problem
Consider the following Poisson problem on the domain : where denotes the boundary of the domain . Following [12], we consider the finite element space: ( 49) and seek U 2 W h for the stabilized Nitsche finite element formulation where and and Here, the outward unit normal n is determined by the level set through the formula n D r kr k 0 ; where is the level set function. In particular, we seek n h 2 OEW h d such that We consider the solution of Equation (50) in two complex domains: the popcorn (45) and the olympic rings, which are constructed by taking the minimum over five tori specified in Equation (44). We set f D 1, D 10:0, and 1 D 0:1.   The solution u in the popcorn geometry is depicted in Figure 15. Figure 15a shows the fictitious domain mesh of the popcorn geometry colored by the value of U . Figure 15b displays the subtesselation of the popcorn surface, and Figure 15c shows a contour plot of U in a slice through the middle of the popcorn. Solving the Poisson Equation (50) in the olympic ring geometry gives the result depicted in Figure 16, where Figure 16a shows the surface subtesselation, and Figure 16b shows the solution U in a slice through the middle of the olympic rings. The UFL input files to specify Equations (50 and (56) in the libCutFEM/FEniCS framework are displayed in Figures 17 and 18, respectively.

Stabilized Nitsche fictitious domain method for the Stokes' problem
Following [29,30], we now employ the ghost penalty technique to solve Stokes' problem   As explained in Section 3, we extend the coercivity of the bilinear form a h . ; / from the physical to fictitious domain using the ghost penalty As the mixed spaces V h P 0;dc h and V h P 1 h are known to violate the Babuška-Brezzi conditions, we stabilize our scheme by the following symmetric pressure penalties: where F i denotes all the interior faces in the mesh domain T . Note that the integral contributions in (57) are always evaluated on the entire face F , even for faces that are intersected by the boundary . Again, such a ghost penalty ensures the robustness and optimal convergence properties of the stabilized scheme irrespective of the location of the boundary; see [29,30] for detailed proofs. Combining these forms, the stabilized Nitsche fictitious domain method for Stokes' problem reads: where A h .U ; P I v; q/ D a h .U ; v/ C b h .v; P / C b h .U ; q/; J h .U ; P I v; q/ D j u .U ; v/ j p .P; q/ and L h .v; q/ D .f; v/ C h 1 v @ n v C qn; g : We present two numerical examples. First, we compute the Stokes flow in the Swiss cheese domain while the domain is rotating along the´-axis at constant angle velocity. As volume force, we take f D .0; 0; 1/. Figure 19 shows both the pressure approximation and the computed velocity. In the second example, we consider the flow through a blood vessel bifurcation with an aneurysm developed at the bifurcation; see Figure 20. Here, the vessel boundary is given as a surface mesh generated from biomedical image data, demonstrating the great potential of CutFEM technologies also to applications where only non-smooth boundary descriptions can be obtained.

Laplace-Beltrami problem on a surface using a bulk-mesh
We consider the Laplace-Beltrami problem for a surface : find u W ! R solving u D f on ; where denotes the popcorn surface described by the level-set function (45). Then, the cut finite element method for the Laplace-Beltrami problem is to seek U 2 W h such that  where the bilinear form A. ; / is defined by The right-hand side f is set to f .x; y;´/ D x C y C´, satisfying the compatibility condition .f; 1/ D 0. For the bulk mesh, the mesh size is approximately h 0:04. Figure 22 depicts the solution U on the popcorn surface (45). Figure 21 shows the UFL input file used in libCutFEM to specify the Laplace-Beltrami problem.

Porous media flow in a domain with cracks
In this section, we give some preliminary numerical results obtained when solving (19) using the stabilized bilinear form   The stabilization term j . ; / is necessary in the case when˛ =˛becomes large, as we are then basically solving a perturbed Laplace-Beltrami problem.
We consider the following model problem in R 2 : a domain .0; 1/ .0; 1/ with Dirichlet boundary conditions u D 1 at x D 0 and u D 0 at x D 1, zero Neumann conditions elsewhere. Choosing D 1 in and an elliptically shaped crack shown in Figure 23 and setting D 10, we obtain the following results when setting˛ D 0 on the lower arc of the crack and˛ 2 ¹1; 2; 4º in the upper arc. A close-up of the computed velocities˛rU close to the crack, using the computational mesh of Figure 24, is shown in Figure 25, together with a corresponding zoom of the mesh in Figure 26, and isoline plots of the pressures U are given in Figure 27. We note that even a small increase in permeability increases the velocity quite noticeably in the crack. Further developments for system of cracks are under way.

CONCLUDING REMARKS
In this paper, we have given an exposition of recent results on CutFEM combined with Nitsche's method and ghost penalty stabilization. The main theoretical ideas have been discussed, but emphasis has been put on implementation issues in the setting of the FEniCS software project.
We have endeavored to show the versatility of the CutFEM method as a method for fictitious domain computations, overlapping mesh methods, multiphysics coupling between a bulk domain and its surfaces or embedded interfaces, and model coupling problems. In particular, we have pointed out the possibility of posing and solving PDEs on interfaces/surfaces as well as in the bulk.
In conclusion, the results reported here (and in the references) show that CutFEM holds great promise as a versatile and powerful mesh-free method posed on a mesh. Challenges currently being investigated include optimization techniques using surface sensitivities, and moving interfaces, for example, for free surface problems.