A Stabilized CutFEM for the Darcy Surface Problem

A Stabilized Cut Finite Element Method for the Darcy Problem on Surfaces

Peter Hansbo Department of Mechanical Engineering, Jönköping University, SE-55111 Jönköping, Sweden. peter.hansbo@ju.se Mats G. Larson  and  André Massing Department of Mathematics and Mathematical Statistics, Umeå University, SE-90187 Umeå, Sweden mats.larson@umu.se andre.massing@umu.se
July 4, 2019

We develop a cut finite element method for the Darcy problem on surfaces. The cut finite element method is based on embedding the surface in a three dimensional finite element mesh and using finite element spaces defined on the three dimensional mesh as trial and test functions. Since we consider a partial differential equation on a surface, the resulting discrete weak problem might be severely ill conditioned. We propose a full gradient and a normal gradient based stabilization computed on the background mesh to render the proposed formulation stable and well conditioned irrespective of the surface positioning within the mesh. Our formulation extends and simplifies the Masud-Hughes stabilized primal mixed formulation of the Darcy surface problem proposed in [28] on fitted triangulated surfaces. The tangential condition on the velocity and the pressure gradient is enforced only weakly, avoiding the need for any tangential projection. The presented numerical analysis accounts for different polynomial orders for the velocity, pressure, and geometry approximation which are corroborated by numerical experiments. In particular, we demonstrate both theoretically and through numerical results that the normal gradient stabilized variant results in a high order scheme.

Key words and phrases:
Surface PDE, Darcy Problem, cut finite element method, stabilization, condition number, a priori error estimates
2010 Mathematics Subject Classification:
Primary 65N30; Secondary 65N85, 58J05.

1. Introduction

1.1. Background and Earlier Work

In recent years, there has been a rapid development of cut finite element methods, also called trace finite element methods, for the numerical solution of partial differential equations (PDEs) on complicated or evolving surfaces embedded into . The main idea is to use the restriction of finite element basis functions defined on a -dimensional background mesh to a discrete, piecewise smooth surface representation which is allowed to cut through the mesh in an arbitrary fashion. The active background mesh then consists of all elements which are cut by the discrete surface, and the finite element space restricted to the active mesh is used to discretize the surface PDE. This approach was first proposed in [32] for the Laplace-Beltrami on a closed surface, see also [3] and the references therein for an overview of cut finite element techniques.

Depending on the positioning of the discrete surface within the background mesh, the resulting system matrix might be severely ill conditioned and either preconditioning [31] or stabilization [4] is necessary to obtain a well conditioned linear system. The stabilization introduced and analyzed in [4] is based on so called face stabilization or ghost penalty, which provides control over the jump in the normal gradient across interior faces in the active mesh. In particular, it was shown that the condition number scaled in an optimal way, independent of how the surface cut the background mesh. Thanks its versatility, the face based stabilization can naturally be combined with discontinuous cut finite element methods as demonstrated in [7]. To reduce the matrix stencil and ease the implementation, a particular simple low order, full gradient stabilization using continuous piecewise linears was developed and analyzed in [8] for the Laplace-Beltrami operator. Then a unifying abstract framework for analysis of cut finite element methods on embedded manifolds of arbitrary codimension was developed in [6] and, in particular, the normal gradient stabilization term was introduced and analyzed. Further developments include convection problems [5, 33], coupled bulk-surface problems [9, 25] and higher order versions of trace fem for the Laplace-Beltrami operator [35, 23]. Moreover, extensions to time-dependent, parabolic-type problems on evolving domains were proposed in [27, 34].

So far, with their many applications to fluid dynamics, material science and biology, e.g., [21, 24, 30, 14, 15, 2], mainly scalar-valued, second order elliptic and parabolic type equations have been considered in the context of cut finite element methods for surface PDEs. Only very recently, vector-valued surface PDEs in combination with unfitted finite element technologies have been considered, for instance in the numerical discretization of surface-bulk problems modeling flow dynamics in fractured porous media [20, 10, 1, 19]. But while these contributions employed cut finite element type methods to discretize the bulk equations, only fitted (mixed and stabilized) finite elements methods on triangulated surfaces have been developed for vector surface equation such as the Darcy surface problem, see for instance [18, 28]. The present contribution is the first where a cut finite element method for a system of partial differential equations on a surfaces involving tangent vector fields of partial differential equations is developed.

1.2. New Contributions

We develop a stabilized cut finite element method for the numerical solution of the Darcy problem on a surface. The proposed variational formulation follows the approach in [28] for the Darcy problem on triangulated surfaces which is based on the stabilized primal mixed formulation by Masud and Hughes [29]. Note that standard mixed type elements are typically not available on discrete cut surfaces. Combining the ideas from [28] with the stabilized full gradient formulations of the Laplace-Beltrami problem from [8, 6], the tangent condition on both the velocity and the pressure gradient is enforced only weakly. When employing finite element function from the full -dimensional background mesh, such a weak enforcement of the tangential condition is convenient and rather natural.

To render the proposed formulation stable and well conditioned irrespective of the relative surface position in the background mesh, we consider two stabilization forms: the full gradient stabilization introduced in [8] which is convenient for low order elements, and the normal gradient stabilization introduced in [6] which also works for higher order elements. Through these stabilizations, we gain control of the variation of the solution orthogonal to the surface, which in combination with the Masud-Hughes variational formulation results in a coercive formulation of the Darcy surface problem. In practice, the exact surface is approximated leading to a geometric error which we take into account in the error analysis. We show stability of the method and establish optimal order a priori error estimates. The presented numerical analysis also accounts for different polynomial orders for the velocity, pressure, and geometry approximation.

1.3. Outline

The paper is organized as follows: In Section 2 we present the Darcy problem on a surface together with its Masud-Hughes weak formulation, followed by the formulation of the cut finite element method in Section 3. In Section 4 we collect a number of preliminary theoretical results, which will be needed in Section 5 where the main a priori error estimates in the energy and norm are established. Finally, in Section 6 we present numerical results illustrating the theoretical findings of this work.

2. The Darcy Problem on a Surface

2.1. The Continuous Surface

In what follows, denotes a smooth compact hypersurface without boundary which is embedded in and equipped with a normal field and signed distance function . Defining the tubular neighborhood of by , the closest point projection is the uniquely defined mapping given by


which maps to the unique point such that for some , see [22]. The closest point projection allows the extension of a function defined on to its tubular neighborhood using the pull back


In particular, we can smoothly extend the normal field to the tubular neighborhood . On the other hand, for any subset such that is bijective, a function on can be lifted to by the push forward satisfying


Finally, for any function space defined on , we denote the space consisting of extended functions by and correspondingly, the notation refers to the lift of a function space  defined on .

2.2. The Surface Darcy Problem

To formulate the Darcy problem on a surface, we first recall the definitions of the surface gradient and divergence. For a function the tangential gradient of can be expressed as


where is the standard gradient and denotes the orthogonal projection of onto the tangent plane of at given by


where is the identity matrix. Since is constant in the normal direction, we have the identity


Next, the surface divergence of a vector field is defined by


With these definitions, the Darcy problem on the surface is to seek the tangential velocity vector field and the pressure such that


Here, is a given function such that and is a given tangential vector field.

2.3. The Masud-Hughes Weak Formulation

We follow [28] and base our finite element method on an extension of the Masud-Hughes weak formulation, originally proposed in [29] for planar domains, to the surface Darcy problem. Using Green’s formula


valid for tangential vector fields , a direct application of the original Masud-Hughes formulation is built upon the fact that the Darcy problem (2.8) solves the weak problem


for test functions and with


As in [28] we enforce the tangent condition on the velocity weakly by using full velocity fields in formulation (2.10) and not only their tangential projection. But in contrast to [28] we also employ the identity (2.6) to replace the pressure related tangent gradients with the full gradient in order to simplify the implementation further. Earlier, such full gradient based formulation have been developed for the Poisson problem on the surface, see [35, 8]. With as the velocity space, as the pressure space and as the total space, the resulting Masud-Hughes weak formulation of the Darcy surface problem (2.8) is to seek satisfying


for all , where


Expanding the forms, the bilinear form and linear form can be rewritten as


and consequently, the bilinear form consists of a symmetric positive definite part and a skew symmetric part. For a more detailed discussion on various weak formulation of the surface Darcy problem, we refer to [28]. Finally, note that since is smooth and is the solution to the elliptic problem , the following elliptic regularity estimate holds


3. Cut Finite Element Methods for the Surface Darcy Problem

3.1. The Discrete Surface and Active Background Mesh

For , we assume that the discrete surface approximation is represented by a piecewise smooth surface consisting of smooth dimensional surface parts associated with a piecewise smooth normal field . On we can then define the discrete tangential projection as the pointwise orthogonal projection on the -dimensional tangent space defined for each and . We assume that:

  • and that the closest point mapping is a bijection for .

  • The following estimates hold


Let be a quasi-uniform mesh, with mesh parameter , which consists of shape regular closed simplices and covers some open and bounded domain containing the embedding neighborhood ). For the background mesh we define the active (background) mesh


see Figure 3.1 for an illustration. Finally, for the domain covered by we introduce the notation

Figure 3.1. Set-up of the continuous and discrete domains. (Left) Continuous surface enclosed by a tubular neighborhood . (Right) Discrete manifold embedded into a background mesh from which the active (background) mesh is extracted.

3.2. Stabilized Cut Finite Element Methods

On the active mesh we introduce the discrete space of continuous piecewise polynomials of order ,


Occasionally, if the order is not of particular importance, we simply drop the superscript . Next, the discrete velocity, pressure and total approximations spaces are defined by respectively


where we explicitly permit different approximation orders and for the velocity and pressure space. As in the continuous case, denotes the average operator on defined by . Now the stabilized, full gradient cut finite element method for the surface Darcy problem (2.8) is to seek such that for all ,




For the stabilization form , two realizations will be proposed in this work. First, we consider a full gradient based stabilization originally introduced for Laplace-Beltrami surface problem in [8],

Second, to devise a higher order approximation scheme, the normal gradient stabilization

first proposed and analyzed in [6] and then later also considered in [23], will be employed. In the remaining work, we will simply write and without superscripts as long as no distinction between the stabilization forms is needed.

Remark 3.1.

For the normal gradient stabilization, any -scaling of the form with gives an eglible stabilization operator, as pointed out in [6]. The condition guarantees that the stabilization result 4.1 for a properly scaled norm holds, the condition on the other hand assures that the condition number of the discrete linear system scales with the mesh size similar to the triangulated surface case. We refer to [6] for further details.

4. Preliminaries

To prepare the forthcoming analysis of the proposed cut finite element method in Section 5, we here collect and state a number of useful definitions, approximation results and estimates. In particular, we introduce suitable continuous and discrete norms, review the construction of a proper interpolation operator and recall the fundamental geometric estimates needed to quantify the quadrature errors introduced by the discretization of .

4.1. Discrete Norms and Poincaré Inequalities

The symmetric parts of the bilinear forms and give raise to the following natural continuous and discrete “energy”-type norms


where denotes the semi-norm induced by the stabilization form ,


To show that actually defines a proper norm, we recall the following result from [6].

Lemma 4.1.

For , the following estimate holds


for with small enough.

Remark 4.2.

Simple counter examples show that the sole expression defines only a semi-norm on . For instance, let be defined as the level set of a smooth, scalar function such that on . Take , and let where is the Lagrange interpolant of . Then on but clearly . Defining gives but in this particular case.

Next, we state a simple surface-based discrete Poincaré estimate. For a proof we refer to [4].

Lemma 4.3.

Let , then it holds


for with chosen small enough.

Finally, the previous two lemma can be combined to obtain the following discrete Poincaré inequality for the discrete “energy” norm .

Theorem 4.4.

For , the following estimate holds


for with small enough.

Figure 4.1. control mechanisms for the full gradient and normal gradient stabilization. (Left) While element has only a small intersection with , there are several neighbor elements in (purple) which share the node and have a “fat” intersection with . The appearance of the full gradient in stabilization allows to integrate along arbitrary directions and thus gives raise to the control of through Lemma 4.1. (Right) The fat intersection property for the discrete “normal” tube guarantees that still a significant portion of can be reached when integrating along normal-like paths which start from and which reside completely inside .

4.2. Trace Estimates and Inverse Inequalities

First, we recall the following trace inequality for


while for the intersection the corresponding inequality


holds whenever is small enough, see [26] for a proof. In the following, we will also need some well-known inverse estimates for :


and the following “cut versions” when


which are an immediate consequence of similar inverse estimates presented in [26].

4.3. Geometric Estimates

We now summarize some standard geometric identities and estimates which typically are used in the numerical analysis of surface PDE discretizations when passing from the discrete surface to the continuous one and vice versa. For a detailed derivation, we refer to [12, 32, 11, 13, 7]. Starting with the Hessian of the signed distance function


the derivative of the closest point projection and of an extended function is given by


The self-adjointness of , , and , and the fact that and then leads to the identities


where the invertible linear mapping


maps the tangential space of at to the tangential space of at . Setting and using the identity , we immediately get that


for any elementwise differentiable function on lifted to . We recall from [22, Lemma 14.7] that for , the Hessian admits a representation


where are the principal curvatures with corresponding principal curvature vectors . Thus


for small enough. In the course of the a priori analysis in Section 5, we will need to estimate various operator compositions involving , the continuous and discrete tangential and normal projection operators. More precisely, using the definition , the following bounds will be employed at several occasions.

Lemma 4.5.

All these estimate have been proved earlier, see [12, 13, 8] and we only include a short proof for the reader’s convenience. We start with the bounds summarized in (4.20). An easy calculation shows that from which the desired bound follows by observing that and thus Next, observe that


Turning to (4.21), the first two bounds follow directly from (4.16) and  (4.19) together with the assumption . Finally, unwinding the definition of , we find that which together with the previously derived estimate for gives the stated operator bound. ∎

The previous lemma allows us to quantify the error introduced by using the full gradient in (3.8) instead of . To do so we decompose the full gradient as with . We then have

Corollary 4.6.

For and it holds


Since , the first estimate follows directly from the identity = from (4.14), while the second estimate is a immediate consequence of (4.20). ∎

Next, for a subset , we have the change of variables formula


with denoting the absolute value of the determinant of . The determinant satisfies the following estimates.

Lemma 4.7.

It holds


Combining the various estimates for the norm and the determinant of shows that for


Next, we observe that thanks to the coarea-formula (cf. Evans and Gariepy [17])

the extension operator defines a bounded operator satisfying the stability estimate


for , where the hidden constant depends only on the curvature of .

4.4. Interpolation Operator

Next, we recall from [16] that for , the Clément interpolant satisfies the local interpolation estimates


where consists of all elements sharing a vertex with . Now with the help of the extension operator , an interpolation operator can be constructed by setting , where we took the liberty of using the same symbol. The resulting interpolation operator satisfies the following error estimate.

Lemma 4.8.

For and , the interpolant defined by satisfies the interpolation estimate


Choosing , it follows directly from combining the trace inequality (4.7), the interpolation estimate (4.31), and the stability estimate (4.30) that the first two terms in the definition of satisfies the desired estimate. Since it is enough to focus on the full gradient stabilization for the remaining part. With the same chain of estimates we find that


which concludes the proof. ∎

5. A Priori Error Estimates

We now state and prove the main a priori error estimates for the stabilized cut finite element method (3.6). The proofs rest upon a Strang-type lemma splitting the total error into an interpolation error, a consistency error arising from the additional stabilization term and finally, a geometric error caused by the discretization of the surface. We start with establishing suitable estimates for the consistency and quadrature error before we present the final a priori error estimates at the end of this section.

5.1. Estimates for the Quadrature and Consistency Error

The purpose of the next lemma is two-fold. First, it shows that the full gradient stabilization will not affect the expected convergence order when low-order elements are used. Second, it demonstrates that only the normal gradient stabilization is suitable for high order discretizations where the geometric approximation order needs to satisfy .

Lemma 5.1.

Let . Then it holds


A simple application of stability estimate (4.30) with shows that for ,


Turning to , the pressure part of the normal gradient stabilization can be estimated by


and similarly, for . ∎

Lemma 5.2.

Let be the solution to weak problem (2.10) and assume that . Then


Furthermore, if and , we have the improved estimate


We start with the term . Unwinding the definition of the linear forms and , we get


For the first term, a change of variables together with estimate (4.27) for the determinant yields


where in the last step, we used the norm equivalences (4.29) and the discrete Poincaré inequality (4.4) to pass to . To estimate , we split into its tangential and normal part


Note that for the tangential field , the identities


hold and thus using once more and the fact that allows us to rewrite as