A Stabilized Cut Finite Element Method for the Darcy Problem on Surfaces
Abstract.
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 MasudHughes 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 estimates2010 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 LaplaceBeltrami 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 LaplaceBeltrami 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 bulksurface problems [9, 25] and higher order versions of trace fem for the LaplaceBeltrami operator [35, 23]. Moreover, extensions to timedependent, parabolictype 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 scalarvalued, second order elliptic and parabolic type equations have been considered in the context of cut finite element methods for surface PDEs. Only very recently, vectorvalued surface PDEs in combination with unfitted finite element technologies have been considered, for instance in the numerical discretization of surfacebulk 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 LaplaceBeltrami 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 MasudHughes 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 MasudHughes 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
(2.1) 
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
(2.2) 
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
(2.3) 
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
(2.4) 
where is the standard gradient and denotes the orthogonal projection of onto the tangent plane of at given by
(2.5) 
where is the identity matrix. Since is constant in the normal direction, we have the identity
(2.6) 
Next, the surface divergence of a vector field is defined by
(2.7) 
With these definitions, the Darcy problem on the surface is to seek the tangential velocity vector field and the pressure such that
(2.8a)  
(2.8b) 
Here, is a given function such that and is a given tangential vector field.
2.3. The MasudHughes Weak Formulation
We follow [28] and base our finite element method on an extension of the MasudHughes weak formulation, originally proposed in [29] for planar domains, to the surface Darcy problem. Using Green’s formula
(2.9) 
valid for tangential vector fields , a direct application of the original MasudHughes formulation is built upon the fact that the Darcy problem (2.8) solves the weak problem
(2.10) 
for test functions and with
(2.11)  
(2.12) 
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 MasudHughes weak formulation of the Darcy surface problem (2.8) is to seek satisfying
(2.13) 
for all , where
(2.14)  
(2.15) 
Expanding the forms, the bilinear form and linear form can be rewritten as
(2.16)  
(2.17) 
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
(2.18) 
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
(3.1)
Let be a quasiuniform 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
(3.2) 
see Figure 3.1 for an illustration. Finally, for the domain covered by we introduce the notation
(3.3) 
3.2. Stabilized Cut Finite Element Methods
On the active mesh we introduce the discrete space of continuous piecewise polynomials of order ,
(3.4) 
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
(3.5) 
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 ,
(3.6) 
where
(3.7)  
(3.8)  
(3.9)  
(3.10) 
For the stabilization form , two realizations will be proposed in this work. First, we consider a full gradient based stabilization originally introduced for LaplaceBeltrami surface problem in [8],
(3.11)  
Second, to devise a higher order approximation scheme, the normal gradient stabilization  
(3.12) 
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
(4.1) 
where denotes the seminorm induced by the stabilization form ,
(4.2) 
To show that actually defines a proper norm, we recall the following result from [6].
Lemma 4.1.
For , the following estimate holds
(4.3) 
for with small enough.
Remark 4.2.
Simple counter examples show that the sole expression defines only a seminorm 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 surfacebased discrete Poincaré estimate. For a proof we refer to [4].
Lemma 4.3.
Let , then it holds
(4.4) 
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
(4.5) 
for with small enough.
4.2. Trace Estimates and Inverse Inequalities
First, we recall the following trace inequality for
(4.6) 
while for the intersection the corresponding inequality
(4.7) 
holds whenever is small enough, see [26] for a proof. In the following, we will also need some wellknown inverse estimates for :
(4.8)  
(4.9) 
and the following “cut versions” when
(4.10) 
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
(4.11) 
the derivative of the closest point projection and of an extended function is given by
(4.12)  
(4.13) 
The selfadjointness of , , and , and the fact that and then leads to the identities
(4.14)  
(4.15) 
where the invertible linear mapping
(4.16) 
maps the tangential space of at to the tangential space of at . Setting and using the identity , we immediately get that
(4.17) 
for any elementwise differentiable function on lifted to . We recall from [22, Lemma 14.7] that for , the Hessian admits a representation
(4.18) 
where are the principal curvatures with corresponding principal curvature vectors . Thus
(4.19) 
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.
(4.20)  
(4.21) 
Proof.
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
(4.22)  
(4.23)  
(4.24) 
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
(4.25) 
Proof.
Next, for a subset , we have the change of variables formula
(4.26) 
with denoting the absolute value of the determinant of . The determinant satisfies the following estimates.
Lemma 4.7.
It holds
(4.27) 
Combining the various estimates for the norm and the determinant of shows that for
(4.28)  
(4.29) 
Next, we observe that thanks to the coareaformula (cf. Evans and Gariepy [17])
the extension operator defines a bounded operator satisfying the stability estimate
(4.30) 
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
(4.31) 
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
(4.32) 
Proof.
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
(4.33)  
(4.34)  
(4.35) 
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 Strangtype 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 twofold. First, it shows that the full gradient stabilization will not affect the expected convergence order when loworder 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
(5.1)  
(5.2) 
Proof.
A simple application of stability estimate (4.30) with shows that for ,
(5.3) 
Turning to , the pressure part of the normal gradient stabilization can be estimated by
(5.4) 
and similarly, for . ∎
Lemma 5.2.
Let be the solution to weak problem (2.10) and assume that . Then
(5.5) 
Furthermore, if and , we have the improved estimate
(5.6) 
Proof.
We start with the term . Unwinding the definition of the linear forms and , we get
(5.7)  
(5.8) 
For the first term, a change of variables together with estimate (4.27) for the determinant yields
(5.9)  
(5.10)  
(5.11) 
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
(5.12) 
Note that for the tangential field , the identities
(5.13) 
hold and thus using once more and the fact that allows us to rewrite as
(5.14)  
(5.15)  
(5.16)  