A tangential and penalty-free finite element method for the surface Stokes problem

In this paper, we consider the surface Stokes problem:

Here, $\gamma\subset\mathbb{R}^{3}$ is a smooth and connected two-dimensional surface with outward unit normal ${\bm{\nu}}$, $\bm{\Pi}={\bf I}-{\bm{\nu}}\otimes{\bm{\nu}}$ is the orthogonal projection onto the tangent space of $\gamma$, and $\nabla_{\gamma}$ and ${\rm div}_{\gamma}$ are the surface gradient and surface divergence operators, respectively. Furthermore, ${\rm Def}_{\gamma}$ is the tangential deformation operator, and the forcing function ${\bm{f}}$ is assumed to be tangential to the surface to ensure well-posedness. Further assumptions and notation are given in Section 2; cf. [20] for derivation of the surface Stokes problem and related models and further discussion of their properties. The system of equations (1.1) is subject to the tangential velocity constraint $\bm{u}\cdot{\bm{\nu}}=0$. To address degeneracies related to Killing fields, i.e., non-trivial tangential vector fields in the kernel of ${\rm Def}_{\gamma}$, we include a zeroth-order mass term in the momentum equations (1.1a) (cf. Remark 4.1).

We consider surface finite element methods (SFEMs), a natural methodology mimicking the variational formulation and built upon classical Galerkin principles. In this approach the domain $\gamma$ is approximated by a polyhedral (or higher-order) surface $\Gamma_{h}$ whose faces constitute the finite element mesh. Similar to the Euclidean setting, SFEMs for the surface Stokes problem based on the standard velocity-pressure formulation must use compatible discrete spaces. Specifically, a discrete inf-sup condition must be satisfied. Given that SFEMs utilize the same framework as their Euclidean counterparts, employing mappings via affine or polynomial diffeomorphisms, one may anticipate that numerous classical inf-sup stable Stokes pairs can be adapted to their surface analogues, readily enabling the construction of stable SFEMs for (1.1).

However, the tangential velocity constraint poses a significant hurdle to constructing stable and convergent SFEMs. As $\Gamma_{h}$ is merely Lipschitz continuous, its outward unit normal is discontinuous at mesh edges and vertices. As a result, the tangential projection of continuous, piecewise smooth functions does not lead to $\bm{H}^{1}$-conforming functions. Moreover, there do not exist canonical, degree-of-freedom-preserving pullbacks for tangential $\bm{H}^{1}$ vector fields, in particular, the Piola transform preserves tangentiality and in-plane normal continuity, but not in-plane tangential continuity. Finally, a continuous, tangential, and piecewise smooth vector field on $\Gamma_{h}$ must necessarily vanish on mesh corners except in exceptional cases where all incident triangles are coplanar. Indeed, at a mesh corner there are at least three faces emanating from a common vertex, whose outward unit normal vectors span $\mathbb{R}^{3}$. Therefore tangentiality of a continuous vector field with respect to each of the three planes implies that it vanishes at the vertex. Thus any piecewise polynomial space simultaneously satisfying both tangentiality and continuity exhibits a locking-type phenomenon with poor approximation properties.

There is a substantial recent literature on numerical approximation of the surface Stokes and related problems such as the surface vector Laplace equation. Most of these circumvent the difficulties described above in one of three ways: by relaxing the pointwise tangential constraint, by relaxing $\bm{H}^{1}$-conformity of the finite element space, or by using a different formulation of the surface Stokes problem. For the former, one can weakly impose the tangential constraint via penalization or Lagrange multipliers [16, 17, 25, 26, 18, 21, 4]. In principle, this allows one to use inf-sup stable Euclidean Stokes pairs to solve the analogous surface problem. However, this methodology requires superfluous degrees of freedom, as the velocity space is approximated by arbitrary vectors in $\mathbb{R}^{3}$ rather than tangential vectors. In addition an unnatural high-order geometric approximation of the unit normal of the true surface is needed to obtain optimal-order approximations. Therefore for problems in which full information of the exact surface is unknown (e.g., free-boundary problem), these penalization schemes lead to SFEMs with sub-optimal convergence properties. However, it was shown recently in [19] that the tangential component of the solution converges optimally for a standard isoparametric geometry approximation in most cases assuming a correct choice of penalty parameters. The only exception is the case where tangential $L_{2}$ errors are considered along with affine (polyhedral) surface approximations. Alternatively, one may relax $\bm{H}^{1}$-conformity and use finite element trial and test functions that are not continuous on the discrete surface $\Gamma_{h}$. In this direction, SFEMs utilizing tangentially- and $\bm{H}({\rm div})$-conforming finite element spaces such as Raviart-Thomas and Brezzi-Douglas-Marini combined with discontinuous Galerkin techniques are proposed and analyzed in [3, 23]; cf. [10] for similar methods for Euclidean Stokes equations. Here, additional consistency, symmetry, and stability terms are added to the method. These terms add some complexity to the implementation, especially for higher-order surface approximations, but are standard in the context of discontinuous Galerkin methods. Optimal-order convergence is observed experimentally for a standard SFEM formulation that does not require higher-order approximations of any geometric information. Discretizations of stream function formulations of the surface Stokes equations have also appeared in the literature [24, 28, 5, 4]. However, as with methods weakly enforcing tangentiality, they require higher-order approximation to the surface normal and in addition require computation of curvature information which can in and of itself be a challenging problem. These methods are also restricted to simply connected surfaces. As a final note, trace SFEMs, in which discretizations of surface PDE are formulated with respect to a background 3D mesh and a corresponding 3D finite element space, are especially important in the context of dynamic surface fluid computations. Trace formulations are well-developed for $H^{1}$ conforming/tangentially nonconforming methods and stream function formulations, but have not yet appeared for $\bm{H}({\rm div})$-conforming methods.

In this paper, we design a SFEM for the surface Stokes problem (1.1) using a strongly tangential finite element space that is based on a conforming, inf-sup stable Euclidean pair. The method is based on the standard variational formulation for the Stokes problem and does not require additional consistency terms or extrinsic penalization. As far as we are aware, this is the first SFEM for the surface Stokes problem with these properties. The key issue that we address is the assignment of degrees of freedom (DOFs) of tangential vector fields at Lagrange nodes, in particular, at vertices of the surface triangulation.

To expand on this last point and to describe our proposed approach, consider a vertex/DOF, call $a$, of the triangulation of the discrete geometry approximation $\Gamma_{h}$, and let $\mathcal{T}_{a}$ denote the set of faces in the triangulation that have $a$ as a vertex. We wish to interpret and define the values of tangential vector fields forming our finite element space at this vertex in a way that ensures the resulting discrete spaces have desirable approximation and weak-continuity properties. As the mesh elements in $\mathcal{T}_{h}$ generally lie in different planes, it is immediate that such vector fields are generally multi-valued at $a$.

Let $\bm{p}=\bm{p}_{\gamma}$ denote the closest point projection onto $\gamma$, and note that, because $\gamma$ is smooth, continuous and tangential vector fields are well-defined and single-valued at $\bm{p}(a)$. Thus, as the Piola transform preserves tangentiality, a natural assignment is to construct finite element functions $\bm{v}$ with the property $\bm{v}|_{K}(a)=\mathcal{P}_{\bm{p}^{-1}}\wideparen{\bm{v}}\big{|}_{\bm{p}(K)}\ \forall K\in\mathcal{T}_{a}$ for some vector field $\wideparen{\bm{v}}$ tangent at $\bm{p}(a)$, where $\mathcal{P}_{\bm{p}^{-1}}$ is the Piola transform of the inverse mapping $\bm{p}^{-1}:\gamma\to\Gamma_{h}$; see Figure 1. Imposing this condition on Lagrange finite element DOFs likely leads to the sought out approximation and weak-continuity properties, and thus, conceptually may lead to convergent SFEMs for (1.1). However, the implementation of the resulting finite element method requires explicit information about the exact surface $\gamma$ and its closest point projection. Therefore this construction is of little practical value.

Instead of this idealized construction, we fix an arbitrary face $K_{a}\in\mathcal{T}_{a}$. Given the value $\bm{v}|_{K_{a}}(a)$ and $K\in\mathcal{T}_{a}$, we then assign $\bm{v}|_{K}(a)=\EuScript{P}_{\bm{p}^{-1}_{K_{a}}}\bm{v}|_{K_{a}}(a)$, the Piola transform of $\bm{v}|_{K_{a}}(a)$ with respect to the inverse of the closest point projection onto the plane containing $K_{a}$; see Figure 2. This transform is linear with a relatively simple formula (cf. Definition 2.3), and it only uses geometric information from $\Gamma_{h}$. Moreover, we show that this construction is only an $O(h^{2})$ perturbation from the idealized setting. As a result, the constructed finite element spaces possess sufficient weak continuity properties to ensure that the resulting scheme is convergent for the surface Stokes problem (1.1).

To clearly communicate the main ideas and to keep technicalities at minimum, we focus on a polyhedral approximation to $\gamma$ and on the lowest-order MINI pair, which in the Euclidean setting takes the discrete velocity space to be the (vector-valued) linear Lagrange space enriched with cubic bubbles, and the discrete pressure space to be the (scalar) Lagrange space. We expect the main ideas to be applicable to other finite element pairs (e.g, Taylor–Hood, Scott-Vogelius [29, 22, 31, 2]), although the stability must be shown on a case-by-case basis. Below we present numerical experiments demonstrating the viability of our approach for $\mathbb{P}^{2}$ surface approximations paired with a $\mathbb{P}^{2}-\mathbb{P}^{1}$ Taylor-Hood finite element space and plan to address generalizations of our approach more fully in future works.

The rest of the paper is organized as follows. In the next section, we introduce the notation and provide some preliminary results. In Section 3, we define the surface finite element spaces based on the classical MINI pair. Here, we show that the spaces have optimal-order approximation properties and are inf-sup stable. We also establish weak continuity properties of the discrete velocity space via an $H^{1}$-conforming relative on the true surface. In Section 4, we define the finite element space and prove optimal-order estimates in the energy norm. Finally in Section 5 we provide numerical experiments which support the theoretical results.

We assume $\gamma\subset\mathbb{R}^{3}$ is a smooth, connected, and orientable two-dimensional surface without boundary. The signed distance function of $\gamma$ is denoted by $d$, which satisfies $d<0$ in the interior of $\gamma$ and $d>0$ in the exterior. We set ${\bm{\nu}}(x)=\nabla d(x)$ to be the outward-pointing unit normal (where the gradient is understood as a column vector) and ${\bf H}(x)=D^{2}d(x)$ the Weingarten map. The tangential projection operator is $\bm{\Pi}={\bf I}-{\bm{\nu}}\otimes{\bm{\nu}}$, where ${\bf I}$ is the $3\times 3$ identity matrix, and the outer product of two vectors ${\bm{a}}$ and ${\bm{b}}$ satisfies $({\bm{a}}\otimes{\bm{b}})_{i,j}=a_{i}b_{j}$. The smoothness of $\gamma$ ensures the existence of $\delta>0$ sufficiently small such that the closest point projection

is well defined in the tubular region $U=\{x\in\mathbb{R}^{3}:\ {\rm dist}(x,\gamma)\leq\delta\}$.

For a scalar function $q:\gamma\to\mathbb{R}$ we define its extension $q^{e}:U\to\mathbb{R}$ via $q^{e}=q\circ\bm{p}$. Likewise, for $\bm{v}=(v_{1},v_{2},v_{3})^{\intercal}:\gamma\to\mathbb{R}^{3}$ its extension $\bm{v}^{e}:U\to\mathbb{R}^{3}$ satisfies $(\bm{v}^{e})_{i}=v_{i}^{e}$ for $i=1,2,3$. Define the surface gradient $\nabla_{\gamma}q=\bm{\Pi}\nabla q^{e}$, and for a (column) vector field $\bm{v}=(v_{1},v_{2},v_{3})^{\intercal}:\gamma\to\mathbb{R}^{3}$, we let $\nabla\bm{v}^{e}=(\nabla v_{1}^{e},\nabla v_{2}^{e},\nabla v_{3}^{e})^{\intercal}$ denote the Jacobian matrix of $\bm{v}^{e}$. We then see $(\nabla\bm{v}^{e}\bm{\Pi})_{i,:}=((\nabla v^{e}_{i})^{\intercal}\bm{\Pi})=(\bm{\Pi}\nabla v^{e}_{i})^{\intercal}=(\nabla_{\gamma}v_{i})^{\intercal}$, i.e., the $ith$ row of $\nabla\bm{v}^{e}\bm{\Pi}$ coincides with $(\nabla_{\gamma}v_{i})^{\intercal}$. The tangential surface gradient (covariant derivative) of $\bm{v}$ is defined by $\nabla_{\gamma}\bm{v}=\bm{\Pi}\nabla\bm{v}^{e}\bm{\Pi}$, and the surface divergence operator of $\bm{v}$ is ${\rm div}_{\gamma}\bm{v}={\rm tr}(\nabla_{\gamma}\bm{v})$. The deformation of a tangential vector field is defined as the symmetric part of its surface gradient, i.e.,

For a matrix field ${\bf A}:\mathbb{R}^{3\times 3}$, the divergence ${\rm div}_{\gamma}{\bf A}$ is understood to act row-wise.

Let $L_{2}(\gamma)$ denote the space of square-integrable functions on $\gamma$ and let $\mathring{L}_{2}(\gamma)$ be the subspace of $L_{2}(\gamma)$ consisting of $L_{2}$-functions with vanishing mean. We let $W_{p}^{m}(\gamma)$ be the Sobolev space of order $m$ and exponent $p$ on $\gamma$ with corresponding norm $\|\cdot\|_{W^{m}_{p}(\gamma)}$. We use the notation $H^{m}(\gamma)=W_{2}^{m}(\gamma)$ with $\|\cdot\|_{H^{m}(\gamma)}=\|\cdot\|_{W^{m}_{2}(\gamma)}$, and the convention $|\cdot|_{H^{0}}=\|\cdot\|_{L_{2}}$, $|\cdot|_{W^{0}_{p}}=\|\cdot\|_{L_{p}}$. Analogous vector-valued spaces are denoted in boldface (e.g., $\bm{L}_{2}(\gamma)=(L_{2}(\gamma))^{3}$ and $\bm{H}^{1}(\gamma)=(H^{1}(\gamma))^{3}$). We let $\bm{H}_{T}^{1}(\gamma)$ be the subspace of $\bm{H}^{1}(\gamma)$ whose members are tangent to $\gamma$, and set

Let $\Gamma_{h}$ be a polyhedral surface approximation of $\gamma$ with triangular faces. We assume that $\Gamma_{h}$ is an $O(h^{2})$ approximation in the sense that $d(x)=O(h^{2})$ for all $x\in\Gamma_{h}$. We further assume $h$ is sufficiently small to ensure $\Gamma_{h}\subset U$, in particular, the closest point projection is well-defined on $\Gamma_{h}$. We denote by $\mathcal{T}_{h}$ the set of faces of $\Gamma_{h}$, and assume this triangulation is shape-regular (i.e., the ratio of the diameters of the inscribed and circumscribed circles of each face is uniformly bounded). For simplicity and to ease the presentation, we further assume that $\mathcal{T}_{h}$ is quasi-uniform, i.e., $h:=\max_{K^{\prime}}{\rm diam}(K^{\prime})\approx{\rm diam}(K)$ for all $K\in\mathcal{T}_{h}$. The image of the mesh elements and the resulting set on the exact surface are given, respectively, by

We use the notation $a\lesssim b$ (resp., $a\gtrsim b$) if there exists a constant $C>0$ independent of the mesh parameter $h$ such that $a\leq Cb$ (resp., $a\geq Cb$). The statement $a\approx b$ means $a\lesssim b$ and $a\gtrsim b$.

Set $\mathcal{V}_{h}$ to be the set of vertices in $\mathcal{T}_{h}$, and for each $K\in\mathcal{T}_{h}$, let $\mathcal{V}_{K}$ denote the set of three vertices of $K$. For each $a\in\mathcal{V}_{h}$, let $\mathcal{T}_{a}\subset\mathcal{T}_{h}$ denote the set of faces having $a$ as a vertex. For $K\in\mathcal{T}_{h}$, we define the patches

so that $\omega_{K}\subset\omega_{K}^{\prime}\subset\Gamma_{h}$. The patches $\omega_{K^{\gamma}}$ and $\omega^{\prime}_{K^{\gamma}}$ associated with $K^{\gamma}=\bm{p}(K)$ are defined analogously.

The (piecewise constant) outward unit normal of $\Gamma_{h}$ is denoted by ${\bm{\nu}}_{h}$, and we shall use the notation ${\bm{\nu}}_{K}={\bm{\nu}}_{h}|_{K}\in\mathbb{R}^{3}$, its restriction to $K\in\mathcal{T}_{h}$. We assume that $|{\bm{\nu}}\circ\bm{p}-{\bm{\nu}}_{h}|\lesssim h$. The tangential projection with respect to $\Gamma_{h}$ is $\bm{\Pi}_{h}={\bf I}-{\bm{\nu}}_{h}\otimes{\bm{\nu}}_{h}$, and we assume there exists $c>0$ independent of $h$ such that ${\bm{\nu}}\cdot{\bm{\nu}}_{h}\geq c>0$ on $\Gamma_{h}$. We let $\mu_{h}(x)$ satisfy $\mu_{h}d\sigma_{h}(x)=d\sigma(\bm{p}(x))$, where $d\sigma$ and $d\sigma_{h}$ are surface measures of $\gamma$ and $\Gamma_{h}$, respectively. In particular,

From [11, Proposition 2.5], we have

and

where $\{\kappa_{1},\kappa_{2}\}$ are the eigenvalues of ${\bf H}$, whose corresponding eigenvectors are orthogonal to ${\bm{\nu}}$. We set $\mu_{K}=\mu_{h}|_{K}$ to be the restriction of $\mu_{h}$ to $K\in\mathcal{T}_{h}$.

Surface differential operators with respect to $\Gamma_{h}$ are denoted and defined analogously to those on $\gamma$. We also set ($m\in\mathbb{N}$)

to be the piecewise $H^{m}$ Sobolev space and norm, respectively. Likewise, $\bm{H}^{m}_{h}(\gamma)$ is the piecewise Sobolev space with respect to $\mathcal{T}_{h}^{\gamma}$ with corresponding norm $\|\bm{v}\|_{H^{m}_{h}(\gamma)}^{2}=\sum_{K\in\mathcal{T}_{h}}\|\bm{v}\|_{H^{m}(K^{\gamma})}^{2}$, and ${\rm Def}_{\gamma,h}$ denotes the piecewise deformation operator with respect to $\mathcal{T}_{h}^{\gamma}$.

We end this section by stating a well-known characterization of $\bm{H}({\rm div}_{\Gamma_{h}};\Gamma_{h})=\{\bm{v}\in\bm{L}_{2}(\Gamma_{h}):\ {\rm div}_{\Gamma_{h}}\bm{v}\in L_{2}(\Gamma_{h})\}$. For each edge $e$ of the mesh, denote by $K_{1}^{e},K_{2}^{e}\in\mathcal{T}_{h}$ the two triangles in the mesh such that $e={\partial}K_{1}^{e}\cap{\partial}K_{2}^{e}$. Let $\bm{n}_{j}^{e}$ denote the outward in-plane normal to ${\partial}K_{j}^{e}$, and note that in general, $\bm{n}^{e}_{1}\neq-\bm{n}^{e}_{2}$ on $e$. Then a vector field $\bm{v}\in\bm{H}^{1}_{h}(\Gamma_{h})$ satisfies $\bm{v}\in\bm{H}({\rm div}_{\Gamma_{h}};\Gamma_{h})$ if and only if [23]

where $\bm{v}_{j}=\bm{v}|_{K_{j}^{e}}$.

By utilizing Lemma 2.5, we can construct tangential finite element spaces on the surface approximation $\Gamma_{h}$ using nodal (Lagrange) basis functions. The essential idea is to enforce continuity at nodal degrees of freedom in a weak sense through the mapping $\mathcal{M}_{a}^{K}$ given in Definition 2.3. Although this procedure does not yield a globally continuous finite element space, it preserves in-plane normal continuity and exhibits weak continuity properties. These properties are generally sufficient for achieving convergence in second-order elliptic problems.

In the following discussion, we focus on the construction of the lowest-order MINI Stokes pair for simplicity [1]. However, we expect that Definition 2.3 and Lemma 2.5 provide a general framework for constructing convergent finite element schemes based on classical and conforming finite element pairs such as Taylor-Hood and Scott-Vogelius [2].