# Low-order divergence-free approximations for the Stokes problem on Worsey-
Farin and Powell-Sabin splits

Maurice Fabien † † Division of Applied Mathematics, Brown University,
Providence, RI 02912, USA<EMAIL_ADDRESS>,  Johnny Guzmán† † Division of
Applied Mathematics, Brown University, Providence, RI 02912, USA
<EMAIL_ADDRESS>,  Michael Neilan ‡ ‡Department of Mathematics,
University of Pittsburgh<EMAIL_ADDRESS>and  Ahmed Zytoon ‡ ‡Department of
Mathematics,University of Pittsburgh<EMAIL_ADDRESS>
###### Abstract.

We derive low-order, inf-sup stable and divergence-free finite element
approximations for the Stokes problem using Worsey-Farin splits in three
dimensions and Powell-Sabin splits in two dimensions. The velocity space
simply consists of continuous, piecewise linear polynomials, where as the
pressure space is a subspace of piecewise constants with weak continuity
properties at singular edges (3D) and singular vertices (2D). We discuss
implementation aspects that arise when coding the pressure space, and in
particular, show that the pressure constraints can be enforced at an algebraic
level.

###### Key words and phrases:

divergence–free, low–order, Worsey–Farin, Powell–Sabin

J. Guzmán was partially supported by NSF grant DMS-1913083,

M. Neilan and A. Zytoon was partially supported by NSF grant DMS-2011733,

##  1\. Introduction

The first inf-sup stable finite element spaces that yield divergence-free
approximations for incompressible fluid models on simplicial triangulations
were given in the classical paper by Scott and Vogelius [15]. Their lowest
order pair, defined on a two-dimensional domain, used quartic Lagrage finite
elements and piecewise cubic polynomials for the discrete velocity and
pressure spaces, respectively; since then, there has been several efforts to
find lower order inf-sup stable and conforming finite element spaces that
produce divergence-free approximations [13, 14, 1, 18, 19, 20, 10, 11, 21].

If we confine ourselves to Lagrange elements for the velocity space, then a
reduction of polynomial degree, while preserving stability, can be done on
certain splits (i.e., refinements) of a simplicial mesh. For example, Zhang
proved stability of the three-dimensional cubic-quadratic pair on Alfeld
splits [18], and stability of the three-dimensional quadratic-linear pair on
Worsey-Farin splits [21]. On the other hand, Zhang suggests that piecewise
linear Lagrange elements for the velocity are unstable on Worsey-Farin splits
[21, p.244]. One of the main objectives of this paper is to show that indeed
the linear case is stable if coupled with the correct pressure space.

The main goal of this paper is to construct inf-sup stable and divergence-free
pairs on simplicial (split) meshes using a linear-constant velocity-pressure
pair on Worsey-Farin splits (3D) and Powell-Sabin splits (2D). As far as we
are aware, this is the first three-dimensional pair with these properties
using linear Lagrange elements for the discrete velocity space. The main
driver of our construction is an explicit characterization of the divergence
operator acting on the Lagrange finite element space developed in [8, 9]. Such
characterizations naturally lead to the definitions of the discrete pressure
space. In particular, the pressure space is a subspace of piecewise constants
with weak continuity constraints on singular edges (Worsey-Farin) and singular
vertices (Powell-Sabin).

Our results in two dimensions are similar to those in [19], where Zhang showed
stability of a Stokes pair with linear Lagrange elements for the velocity
space on Powell–Sabin splits. However, the pressure space was not
characterized in [19] and as a result, a pressure basis is not explicitly
constructed. Due to the implicit definition of the pressure space, one is
forced to solve the resulting algebraic system by the iterated penalty method
(IPM). In contrast, the explicit characterization of the pressure space we
provide opens the door to a library of solvers based on the standard mixed
formulation.

An advantage of the proposed schemes is their low computational cost and
simple velocity spaces, which are supported in finite element software
libraries. On the other hand, the discrete pressure spaces are subspaces of
piecewise constants with weak continuity properties at singular edges and
singular vertices; these spaces are nonstandard, and in particular are not
readily available on computational software packages. Nonetheless, we show
that these weak continuity properties can be enforced entirely at the
algebraic level. One can simply form the algebraic saddle point problem using
the full (unstable) linear-constant pair, and then perform elementary row and
column operations to this system to enforce the weak continuity constraint. As
a result, the proposed discretizations can be implemented on standard finite
element software packages (e.g., FEniCS).

We note that degrees of freedom and commuting projections were given in [8,
9]. In fact, an entire exact sequence of spaces were presented. However,
stability (e.g., inf-sup stability) was not shown, which is something we carry
out here.

Moreover, we provide numerical experiments to validate our theoretical
results. In particular, we show that the computed solution is satisfying the
standard error estimates in mixed formulation. Also, we provide a time
comparison between the stiffness matrix modification method described above
and the iterated penalty method.

The paper is organized as follows. In the next section we give notation that
will be used in the rest of the paper. In Section 2 we provide the notation
used throughout the paper and introduce the Stokes problem. In Section 3 prove
the inf-sup stability of a low-order finite element pair on Powell-Sabin
splits. In Section 4 we prove stability of the analogous three-dimensional
pair on Worsey-Farin splits. In Section 5 we discuss implementation aspects on
Powell-Sabin splits and in Section 6 we do the same for Worsey-Farin splits.
Finally, in Section 7 we provide numerical experiments.

##  2\. Preliminaries

In this section we develop basic notation that we use throughout the paper. We
provide this in the following list:

  * •

$\mathcal{T}_{h}$ is a shape-regular, simplicial triangulation of a
contractible polytope $\Omega\subset\mathbb{R}^{d}\ (d=2,3)$.

  * •

$h_{T}={\rm diam}(T)$ for all $T\in\mathcal{T}_{h}$ and
$h=\max_{T\in\mathcal{T}_{h}}h_{T}$.

  * •

For an $n$-dimensional simplex $S$ ($n\leq d$) and $m\in\\{0,\ldots n\\}$,
denote by $\Delta_{m}(S)$ the set of $m$-dimensional simplices of $S$.

  * •

$\EuScript{P}_{r}(S)$ denotes the space of polynomials of degree $\leq r$ with
domain $S$. Analogous spaces of vector-valued spaces are given in boldface,
e.g., $\bm{\EuScript{P}}_{r}(S)=[\EuScript{P}_{r}(S)]^{d}$.

  * •

We define the following function spaces on $\Omega$:

$\displaystyle L^{2}(\Omega)$ $\displaystyle:=\\{w:\Omega\mapsto\mathbb{R}:\
\|w\|_{{\it L}^{2}(\Omega)}:=(\int_{\Omega}\ {|w|^{2}}\,dx)^{1/2}<\infty\\},$
$\displaystyle H^{m}(\Omega)$ $\displaystyle:=\\{w:\Omega\mapsto\mathbb{R}:\
\|w\|_{{\it H}^{m}(\Omega)}:=(\sum_{|\beta|\leq m}\|{\it D}^{\beta}w\|_{{\it
L}^{2}(\Omega)}^{2})^{1/2}<\infty\\},$

and the spaces with boundary conditions:

$\displaystyle L_{0}^{2}(\Omega)$ $\displaystyle:=\\{w\in{\it L}^{2}(\Omega):\
\int_{\Omega}\ {w}\,dx=0\\},$ $\displaystyle H_{0}^{m}(\Omega)$
$\displaystyle:=\\{{w}\in{\it H}^{m}(\Omega):\
{D}^{\beta}{w}|_{\partial\Omega}=0,\forall\beta:\ |\beta|\leq m-1\\}.$

Analogous vector-valued spaces are denoted in boldface, e.g.,
$\bm{H}^{1}(\Omega)=[H^{1}(\Omega)]^{d}$.

  * •

For a simplicial triangulation $\mathcal{Q}_{h}$ of $\Omega$ we define the
spaces of piecewise polynomials

$\displaystyle\EuScript{P}_{k}(\mathcal{Q}_{h})=\prod_{K\in\mathcal{Q}_{h}}\EuScript{P}_{k}(K),\qquad$
$\displaystyle\bm{\EuScript{P}}_{k}(\mathcal{Q}_{h})=\prod_{K\in\mathcal{Q}_{h}}\bm{\EuScript{P}}_{k}(K),$
$\displaystyle\mathring{\EuScript{P}}_{k}(\mathcal{Q}_{h})=\EuScript{P}_{k}(\mathcal{Q}_{h})\cap
L^{2}_{0}(\Omega),\qquad$
$\displaystyle\mathring{\bm{\EuScript{P}}}_{k}^{c}(\mathcal{Q}_{h})=\bm{\EuScript{P}}_{k}(\mathcal{Q}_{h})\cap\bm{H}^{1}_{0}(\Omega).$

Thus $\mathring{\EuScript{P}}_{k}(\mathcal{Q}_{h})$ consists of piecewise
polynomials of degree $\leq k$ with respect to the triangulation
$\mathcal{Q}_{h}$ with mean zero, and
$\mathring{\bm{\EuScript{P}}}_{k}^{c}(\mathcal{Q}_{h})$ is the space of
continuous, piecewise polynomials of degree $\leq k$ with vanishing trace
(i.e., the $k$th degree vector-valued Lagrange finite element space).

  * •

The constant $C$ denotes a generic positive constant, independent of the mesh
parameter $h$.

###  2.1. The Stokes problem

The Stokes equations defined on a polygonal domain
$\Omega\subset\mathbb{R}^{d}$ with Lipschitz continuous boundary
${\partial}\Omega$ is given by the system of equations

(2.1a) $\displaystyle-\nu\Delta{\bm{u}}+{\nabla}p$
$\displaystyle={{\bm{f}}}\qquad$ $\displaystyle\text{in }\Omega,$ (2.1b)
$\displaystyle{\rm div}\,\bm{u}$ $\displaystyle=0\qquad$
$\displaystyle\text{in }\Omega,$ (2.1c) $\displaystyle\bm{u}$
$\displaystyle=0\qquad$ $\displaystyle\text{on }{\partial}\Omega,$

where the velocity $\bm{u}=(u_{1},...,u_{d})^{\intercal}$ and pressure $p$ are
functions of $x=(x_{1},...,x_{d})^{\intercal}$, and $\nabla$, $\Delta$ denote
the gradient operator and vector Laplacian operator with respect to $x$
respectively. In (2.1a), $\nu$ is the viscosity.

The weak formulation for (2.1) reads: Find
$({\bm{u}},p)\in\bm{H}_{0}^{1}(\Omega)\times L_{0}^{2}(\Omega)$ such that
$\forall({\bm{v}},q)\in\bm{H}_{0}^{1}(\Omega)\times{L}_{0}^{2}(\Omega)$ we
have

(2.2a)
$\displaystyle\nu\int_{\Omega}{\nabla{\bm{u}}:\nabla{\bm{v}}}-\int_{\Omega}{({\rm
div}\,\bm{v})p}$ $\displaystyle=\int_{\Omega}{{{\bm{f}}\cdot{\bm{v}}}},$
(2.2b) $\displaystyle\int_{\Omega}{({\rm div}\,{\bm{u}})q}$ $\displaystyle=0.$

It is well known that the problem (2.2) has a unique solution [7].

Let
$\mathring{\bm{V}}_{h}\times\mathring{Y}_{h}\subset\bm{H}^{1}_{0}(\Omega)\times
L^{2}_{0}(\Omega)$ be a conforming and finite dimensional pair. Then the
finite element method for (2.1), based on the standard velocity-pressure
formulation, seeks
$(\bm{u}_{h},p_{h})\in\mathring{\bm{V}}_{h}\times\mathring{Y}_{h}$ such that
$\forall(\bm{v},q)\in\mathring{\bm{V}}_{h}\times\mathring{Y}_{h}$ we have

(2.3a)
$\displaystyle\nu\int_{\Omega}{\nabla{\bm{u}_{h}}:\nabla{\bm{v}}}-\int_{\Omega}{({\rm
div}\,\bm{v})p_{h}}$ $\displaystyle=\int_{\Omega}{{{\bm{f}}\cdot{\bm{v}}}},$
(2.3b) $\displaystyle\int_{\Omega}{({\rm div}\,{\bm{u}_{h}})q}$
$\displaystyle=0.$

The discrete problem (2.3) is well-posed if and only if the pair satisfies the
inf-sup condition

$\sup_{0\neq\bm{v}\in\mathring{\bm{V}}_{h}}\frac{\int_{\Omega}({\rm
div}\,\bm{v})q}{\|\nabla\bm{v}\|_{L^{2}(\Omega)}}\geq\beta\|q\|_{L^{2}(\Omega)},$

for some $\beta>0$.

##  3\. Inf-sup stability on Powell–Sabin Triangulations

Let $\mathcal{T}_{h}$ be a simplicial, shape-regular triangulation of
$\Omega\subset\mathbb{R}^{2}$. The Powell–Sabin refinement of
$\mathcal{T}_{h}$, denoted by $\mathcal{T}_{h}^{\rm PS}$ is obtained by the
following procedure. First connect the incenter of each triangle
$T\in\mathcal{T}_{h}$ with its three vertices. Next, the incenters of each
adjacent pair of triangles are connected with an edge. If $T$ has a boundary
edge, then the midpoint of the boundary edge is connected to the incenter of
$T$. Thus, this procedure splits each triangle in $\mathcal{T}_{h}$ into six
subtriangles; cf. Figure 1.

Figure 1.  Two triangles in the mesh $\mathcal{T}_{h}$ (left) and their
Powell–Sabin refinement (right). The singular vertices in the Powell–Sabin
refinement are depicted in red.

One feature of the Powell-Sabin triangulation is the presence of singular
vertices.

######  Definition 3.1.

We say that a vertex in a simplicial triangulation is singular if the edges
meeting at the vertex fall on exactly two straight lines.

Let $\mathcal{S}^{I}_{h}$ and $\mathcal{S}^{B}_{h}$ denote the sets of
interior and boundary singular vertices in $\mathcal{T}_{h}^{\rm PS}$,
respectively, and set
$\mathcal{S}_{h}=\mathcal{S}^{I}_{h}\cup\mathcal{S}_{h}^{B}$. Note that the
cardinalities of the sets $\mathcal{S}_{h}^{I}$ and $\mathcal{S}^{B}_{h}$
correspond to the number of interior and boundary edges in $\mathcal{T}_{h}$,
respectively. For $z\in\mathcal{S}_{h}^{I}$, we denote by
$\mathcal{T}_{z}\subset\mathcal{T}_{h}^{\rm PS}$ the set of four triangles
that have $z$ as a vertex. We write
$\mathcal{T}_{z}=\\{K_{z}^{(1)},K_{z}^{(2)},K_{z}^{(3)},K_{z}^{(4)}\\}$ with
$K_{z}^{(j)}\in\mathcal{T}_{h}^{\rm PS}$, labeled such that $K_{z}^{(j)}$ and
$K_{z}^{(j+1)}$ have a common edge. Likewise, for $z\in\mathcal{S}_{h}^{B}$,
we set
$\mathcal{T}_{z}=\\{K_{z}^{(1)},K_{z}^{(2)}\\}\subset\mathcal{T}_{h}^{\rm PS}$
to denote the two triangles that have $z$ as a vertex. We set
$n_{z}=|\mathcal{T}_{z}|$, i.e., $n_{z}=4$ if $z\in\mathcal{S}_{h}^{I}$ and
$n_{z}=2$ if $z\in\mathcal{S}_{h}^{B}$.

###  3.1. Finite Element Spaces on Powell–Sabin Triangulations

It is well known that, on singular vertices, the divergence operator acting on
the Lagrange finite element space has “weak continuity properties”; the
precise meaning of this statement is provided in the following lemma. Its
proof can be found in, e.g., [15, 12].

######  Lemma 3.2.

For $z\in\mathcal{S}_{h}$, and piecewise smooth function $q$, define

$\theta_{z}(q):=\left\\{\begin{array}[]{ll}q|_{K_{z}^{(1)}}(z)-q|_{K_{z}^{(2)}}(z)+q|_{K_{z}^{(3)}}(z)-q|_{K_{z}^{(4)}}(z)&z\in\mathcal{S}_{h}^{I},\\\
q|_{K_{z}^{(1)}}(z)-q|_{K_{z}^{(2)}}&z\in\mathcal{S}_{h}^{B}.\end{array}\right.$

Then there holds $\theta_{z}({\rm div}\,\bm{v})=0$ for all
$\bm{v}\in\mathring{\bm{\EuScript{P}}}_{k}^{c}(\mathcal{T}_{h}^{\rm PS})$.

Based on Lemma 3.2, and on a discrete de Rham complex, divergence–free finite
element pairs for the Stokes problem have been constructed and analyzed in [8,
19] on Powell–Sabin triangulations. Here, we focus on the lowest-order case,
where the velocity space is the linear Lagrange space, and the pressure space
is a subspace of piecewise constants with a weak continuity property at
singular vertices. We first define the spaces without boundary conditions and
then the ones with boundary conditions.

$\displaystyle\bm{V}_{h}^{\rm PS}:$
$\displaystyle=\bm{\EuScript{P}}^{c}_{1}(\mathcal{T}_{h}^{\rm PS}),$
$\displaystyle Y_{h}^{\rm PS}:$
$\displaystyle=\\{q\in\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm PS}):\
\theta_{z}(q)=0\ \forall z\in\mathcal{S}_{h}^{I}\\}.$

We now define an intermediate pressure space

$\displaystyle\hat{Y}_{h}^{\rm PS}:$
$\displaystyle=\\{q\in\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm PS}):\
\theta_{z}(q)=0\ \forall z\in\mathcal{S}_{h}\\}.$

The spaces with boundary conditions are

$\displaystyle\mathring{\bm{V}}_{h}^{\rm PS}:$ $\displaystyle=\bm{V}_{h}^{\rm
PS}\cap\bm{H}_{0}^{1}(\Omega),$ $\displaystyle\mathring{Y}_{h}^{\rm PS}:$
$\displaystyle=\hat{Y}_{h}^{\rm PS}\cap L_{0}^{2}(\Omega).$

###  3.2. Stability of $\mathring{\bm{V}}_{h}^{\rm
PS}\times\mathring{Y}_{h}^{\rm PS}$

The stability of the pair $\mathring{\bm{V}}_{h}^{\rm
PS}\times\mathring{Y}_{h}^{\rm PS}$ is implicitly shown in [2]. Here we give
more details.

######  Theorem 3.3.

There holds ${\rm div}\,\mathring{\bm{V}}_{h}^{\rm
PS}\subseteq\mathring{Y}_{h}^{\rm PS}$. Moreover, the pair
$\mathring{\bm{V}}_{h}^{\rm PS}\times\mathring{Y}_{h}^{\rm PS}$ is an inf-sup
stable pair for the Stokes problem, i.e.,

$\sup_{0\neq\bm{v}\in\mathring{\bm{V}}_{h}^{\rm PS}}\frac{\int_{\Omega}({\rm
div}\,\bm{v})q}{\|\nabla\bm{v}\|_{L^{2}(\Omega)}}\geq\beta\|q\|_{L^{2}(\Omega)}\quad\forall
q\in\mathring{Y}_{h}^{\rm PS},$

for some $\beta>0$ independent of the mesh parameter $h$.

###### Proof.

For fixed $q\in\mathring{Y}_{h}^{\rm PS}$, let
$\bm{w}\in\bm{H}^{1}_{0}(\Omega)$ satisfy ${\rm div}\,\bm{w}=q$ and
$\|\bm{w}\|_{H^{1}(\Omega)}\leq C\|q\|_{L^{2}(\Omega)}$. [2, Theorem 3.4]
gives the existence of an operator $\bm{\Pi}_{h}$ such that
$\bm{\Pi}_{h}\bm{w}\in\bm{V}_{h}^{\rm PS}$, ${\rm
div}\,\bm{\Pi}_{h}\bm{w}={\rm div}\,\bm{w}=q$, and

$\|\bm{w}-\bm{\Pi}_{h}\bm{w}\|_{L^{2}(T)}\leq
Ch_{T}\big{[}\|\bm{w}\|_{H^{1}(\omega(T))}+\|{\rm
div}\,\bm{w}\|_{L^{2}(T)}\big{]}\qquad\forall T\in\mathcal{T}_{h},$

where
$\omega(T)=\mathop{\cup_{T^{\prime}\in\mathcal{T}_{h}}}_{\bar{T}^{\prime}\cap\bar{T}\neq\emptyset}T^{\prime}$.
Standard arguments then show $\|\nabla\bm{\Pi}_{h}\bm{w}\|_{L^{2}(\Omega)}\leq
C\big{(}\|\bm{w}\|_{H^{1}(\Omega)}+\|{\rm
div}\,\bm{w}\|_{L^{2}(\Omega)}\big{)}\leq C\|q\|_{L^{2}(\Omega)}$. Setting
$\bm{v}=\bm{\Pi}_{h}\bm{w}$, we have

$\displaystyle\int_{\Omega}({\rm div}\,\bm{v})q\,dx=\int_{\Omega}({\rm
div}\,\bm{\Pi}_{h}\bm{w})q\,dx=\|q\|_{L^{2}(\Omega)}^{2}\geq
C^{-1}\|q\|_{L^{2}(\Omega)}\|\nabla\bm{v}\|_{L^{2}(\Omega)}.$

This yields the inf-sup stability result with $\beta=C^{-1}$. ∎

##  4\. Inf-sup stability on Worsey-Farin splits

Let $\mathcal{T}_{h}$ be a simplicial triangulation of a polyhedral domain
$\Omega\subset\mathbb{R}^{3}$. The Worsey-Farin triangulation
$\mathcal{T}_{h}^{\rm WF}$ is obtained by splitting each tetrahedron into
twelve sub-tetrahedra by the following procedure. Similar to the Powell–Sabin
case, for each $T\in\mathcal{T}_{h}$, we connect the incenter of $T$ to its
vertices. Next, the incenters of neighboring pairs of tetrahedra are connected
with an edge. This creates a face split point (a vertex) on each face of $T$.
If $T$ has a boundary face, then we connect the incenter of $T$ to the
barycenter of the face by a line. Finally, the face split points are connected
to the vertices of the face. For each $T\in\mathcal{T}_{h}$, we denote by
$T^{\rm WF}$ the triangulation resulting from local Worsey-Farin refinement of
$T$, i.e.,

$T^{\rm WF}=\\{K\in\mathcal{T}_{h}^{\rm WF}:\ K\subset\bar{T}\\}.$

######  Definition 4.1.

An edge in a 3D simplicial triangulation is called singular if the faces
meeting at the edge fall on exactly two planes.

By construction, the Worsey-Farin triangulation contains many singular edges;
for each face in the unrefined triangulation $\mathcal{T}_{h}$, there are
three associated singular edges in $\mathcal{T}_{h}^{\rm WF}$.

Let $\mathcal{E}^{\mathcal{S}}_{h}$ denote the set of singular edges in
$\mathcal{T}_{h}^{\rm WF}$, and let $\mathcal{E}^{\mathcal{S},I}_{h}$ and
$\mathcal{E}_{h}^{\mathcal{S},B}$ denote the sets of interior and boundary
singular edges, respectively. For each $e\in\mathcal{E}^{\mathcal{S}}_{h}$,
let $\mathcal{T}_{e}=\\{K_{e}^{(1)},\ldots,K_{e}^{(n_{e})}\\}$ denote the set
of tetrahedra that have $e$ as an edge. Here, $n_{e}=4$ if $e$ is an interior
edge, and $n_{e}=2$ if $e$ is a boundary edge. We assume the tetrahedra are
labeled such that $K_{e}^{(j)}$ and $K_{e}^{(j+1)}$ share a common face.

###  4.1. Finite Element Spaces on Worsey–Farin Triangulations

Similar to the two-dimensional case, the divergence operator acting on the
Lagrange finite element space has weak continuity properties on singular edges
(cf. [9]).

######  Lemma 4.2.

For $e\in\mathcal{E}^{\mathcal{S}}_{h}$, and a piecewise smooth function $q$,
define

$\theta_{e}(q)=\left\\{\begin{array}[]{ll}q_{e}^{(1)}|_{e}-q_{e}^{(2)}|_{e}+q_{e}^{(3)}|_{e}-q_{e}^{(4)}|_{e}&e\in\mathcal{E}^{\mathcal{S},I}_{h},\\\
q_{e}^{(1)}|_{e}-q_{e}^{(2)}|_{e}&e\in\mathcal{E}^{\mathcal{S},B}_{h},\end{array}\right.$

where $q_{e}^{(j)}=q|_{K_{e}^{(j)}}$. Then there holds $\theta_{e}({\rm
div}\,\bm{v})=0$ for all
$\bm{v}\in\mathring{\bm{\EuScript{P}}}_{k}^{c}(\mathcal{T}_{h}^{\rm WF})$.

Analogous to the Powell–Sabin case, we define the finite element spaces to
discretize the Stokes problem on Worsey-Farin splits. We first define the
spaces without boundary conditions

$\displaystyle\bm{V}_{h}^{\rm WF}$
$\displaystyle=\bm{\EuScript{P}}_{1}^{c}(\mathcal{T}_{h}^{\rm WF}),$
$\displaystyle Y_{h}^{\rm WF}$
$\displaystyle=\\{q\in\mathring{\EuScript{P}}_{0}(\mathcal{T}_{h}^{\rm WF}):\
\theta_{e}(q)=0\ \forall e\in\mathcal{E}^{\mathcal{S},I}_{h}\\}.$

Then, we define an intermediate pressure space

$\displaystyle\hat{Y}_{h}^{\rm WF}$
$\displaystyle=\\{q\in\mathring{\EuScript{P}}_{0}(\mathcal{T}_{h}^{\rm WF}):\
\theta_{e}(q)=0\ \forall e\in\mathcal{E}^{\mathcal{S}}_{h}\\}.$

We now define the spaces with boundary conditions

$\displaystyle\mathring{\bm{V}}_{h}^{\rm WF}$ $\displaystyle=\bm{V}_{h}^{\rm
WF}\cap\bm{H}_{0}^{1}(\Omega),$ $\displaystyle\mathring{Y}_{h}^{\rm WF}$
$\displaystyle=\hat{Y}_{h}^{\rm WF}\cap L_{0}^{2}(\Omega).$

###  4.2. Stability of $\mathring{\bm{V}}_{h}^{\rm
WF}\times\mathring{Y}_{h}^{\rm WF}$

In this section, we show the pair $\mathring{\bm{V}}_{h}^{\rm
WF}\times\mathring{Y}_{h}^{\rm WF}$ is inf-sup stable. First we introduce some
notation.

Let $T\in\mathcal{T}_{h}$, and let $T^{\rm A}$ denote the local triangulation
of $T$ consisting of four tetrahedra, obtained by connecting the vertices of
$T$ with its incenter, i.e., $T^{\rm A}$ denotes the Alfeld split of $T$. For
a face $F\subset T$, denote by $F^{\rm CT}$ the set of three triangles formed
from $F$ by the Worsey-Farin refinement, i.e., $F^{\rm CT}$ is the Clough-
Tocher refinement of $F$. We denote by $\Delta_{1}^{I}(F^{\rm CT})$ the set of
three interior edges in $F^{\rm CT}$, and let $e_{F}\in\Delta_{1}^{I}(F^{\rm
CT})$ denote an arbitrary, fixed interior edge of $F^{\rm CT}$.

######  Definition 4.3.

Consider the triangulation $F^{CT}$ of a face $F\in\Delta_{2}(T)$, and let the
three triangles of $F^{\rm CT}$ be labeled $Q_{1},Q_{2},Q_{3}$. Let
$e={\partial}Q_{1}\cap{\partial}Q_{2}$ be an internal edge, let $\bm{t}$ be
the unit vector tangent to $e$ pointing away from the split point $m_{F}$, and
let $\bm{s}$ be the unit vector orthogonal to $\bm{t}$ such that
$\bm{s}\times\bm{n}_{F}=\bm{t}$. Then the jump of a piecewise smooth function
$p$ across $e$ is defined as

$\left[\left[p\right]\right]_{e}=\big{(}p|_{Q_{1}}-p|_{Q_{2}}\big{)}\bm{s}.$

We now state the degrees of freedom for the spaces $\mathring{\bm{V}}_{h}^{\rm
WF}$ and $\mathring{Y}_{h}^{\rm WF}$. The proofs of the following two lemmas
are given in [9, Lemmas 5.11–5.12].

######  Lemma 4.4.

A function $\bm{v}\in\bm{V}_{h}^{\rm WF}$ is uniquely determined by the values

$\displaystyle\bm{v}(a)\qquad$ $\displaystyle\forall a\in\Delta_{0}(T),$
$\displaystyle\int_{F}(\bm{v}\cdot\bm{n}_{F})\qquad$ $\displaystyle\forall
F\in\Delta_{2}(T),$ $\displaystyle\int_{e}\left[\left[{\rm
div}\,\bm{v}\right]\right]_{e}\qquad$ $\displaystyle\forall
e\in\Delta_{1}^{I}(F^{\rm CT})\backslash\\{e_{F}\\},\ \forall
F\in\Delta_{2}(T),$ $\displaystyle\int_{T}({\rm div}\,\bm{v})p\qquad$
$\displaystyle\forall
p\in\mathring{\mathcal{V}}_{0}(T):=\EuScript{P}_{0}(T^{\rm A})\cap
L^{2}_{0}(T).$

for each $T\in\mathcal{T}_{h}$.

######  Lemma 4.5.

A function $q\in Y_{h}^{\rm WF}$ is uniquely determined by

$\displaystyle\int_{e}\left[\left[q\right]\right]_{e}\qquad$
$\displaystyle\forall e\in\Delta_{1}^{I}(F^{\rm CT})\backslash\\{e_{F}\\},\
\forall F\in\Delta_{2}(T),$ $\displaystyle\int_{T}qp\qquad$
$\displaystyle\forall p\in\EuScript{P}_{0}(T^{\rm A}).$

for all $T\in\mathcal{T}_{h}$.

If we restrict ourselves to $\hat{Y}_{h}^{\rm WF}$ then we only take interior
faces $F$ in the first set of degrees of freedom.

######  Proposition 4.6.

Let $\bm{v}\in\bm{V}_{h}^{\rm WF}$ and $T\in\mathcal{T}_{h}$. For $m=0,1$,
there holds

$\displaystyle|\bm{v}|_{H^{m}(T)}$ $\displaystyle\leq
Ch_{T}^{-1-2m}\Big{(}h_{T}^{4}\sum_{a\in\Delta_{0}(T)}|\bm{v}(a)|^{2}+\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}\bm{v}\cdot\bm{n}_{F}\Big{|}^{2}+h_{T}^{3}\|{\rm
div}\,\bm{v}\|_{L^{2}(T)}^{2}\Big{)}.$

###### Proof.

Let $\hat{T}$ be the reference tetrahedron, and let $F_{T}:\hat{T}\to T$ be an
affine bijection with $F_{T}(\hat{x})=A_{T}\hat{x}+b_{T}$ with
$A_{T}\in\mathbb{R}^{3\times 3}$ and $b_{T}\in\mathbb{R}^{3}$. We define
$\hat{\bm{v}}:\hat{T}\to\mathbb{R}^{3}$ via the Piola transform

$\bm{v}(x)=\frac{A_{T}\hat{\bm{v}}(\hat{x})}{\det(A_{T})},\qquad
x=F_{T}(\hat{x}).$

Let $\hat{T}^{\rm WF}$ be the split of $\hat{T}$ induced by $T^{\rm WF}$ and
the mapping $F_{T}^{-1}$, i.e.,

$\hat{T}^{\rm WF}=\\{F_{T}^{-1}(K):\ K\in T^{\rm WF}\\}.$

Then $\hat{\bm{v}}$ is a continuous piecewise linear polynomial with respect
to $\hat{T}^{\rm WF}$, and therefore by equivalence of norms, and Lemma 4.4,

$\displaystyle|\hat{\bm{v}}|_{H^{m}(\hat{T})}^{2}$ $\displaystyle\leq
C\Big{(}\sum_{\hat{a}\in\Delta_{0}(\hat{T})}|\hat{\bm{v}}(\hat{a})|^{2}+\sum_{\hat{F}\in\Delta_{2}(\hat{T})}\Big{|}\int_{\hat{F}}\hat{\bm{v}}\cdot\hat{\bm{n}}_{\hat{F}}\Big{|}^{2}$
$\displaystyle\qquad+\sum_{\hat{F}\in\Delta_{2}(\hat{T})}\sum_{\hat{e}\in\Delta^{I}_{1}(\hat{F}^{\rm
CT})\backslash\\{\hat{e}_{\hat{F}}\\}}\Big{|}\int_{\hat{e}}\left[\left[\widehat{{\rm
div}\,}\hat{\bm{v}}\right]\right]_{\hat{e}}\Big{|}^{2}+\mathop{\sup_{\hat{p}\in\mathring{\mathcal{V}}_{0}(\hat{T})}}_{\|\hat{p}\|_{L^{2}(\hat{T})}=1}\Big{|}\int_{\hat{T}}(\widehat{\rm
div}\,\hat{\bm{v}})\hat{p}\Big{|}^{2}\Big{)}.$

By well-known properties of the Piola transform, we have

$\displaystyle{\rm div}\,\bm{v}(x)=\frac{1}{\det(A_{T})}\widehat{\rm
div}\,\hat{\bm{v}}(\hat{x}),\qquad\int_{F}\bm{v}\cdot\bm{n}_{F}=\int_{\hat{F}}\hat{\bm{v}}\cdot\hat{\bm{n}}_{\hat{F}}.$

Thus, we have

$\displaystyle|\hat{\bm{v}}|_{H^{m}(\hat{T})}^{2}$ $\displaystyle\leq
C\Big{(}\sum_{a\in\Delta_{0}(T)}|\det(A_{T})A_{T}^{-1}\bm{v}(a)|^{2}+\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}\bm{v}\cdot\bm{n}_{F}\Big{|}^{2}$
$\displaystyle\qquad+|\det(A_{T})|^{2}\sum_{F\in\Delta_{2}(T)}\sum_{e\in\Delta_{1}(F)\backslash\\{e_{F}\\}}\Big{|}\frac{|\hat{e}|}{|e|}\int_{e}\left[\left[{{\rm
div}\,}\bm{v}\right]\right]_{e}\Big{|}^{2}+\mathop{\sup_{\hat{p}\in\mathring{\mathcal{V}}_{0}(\hat{T})}}_{\|\hat{p}\|_{L^{2}(\hat{T})}=1}\Big{|}\int_{\hat{T}}(\widehat{\rm
div}\,\hat{\bm{v}})\hat{p}\Big{|}^{2}\Big{)}.$

Next, for $\hat{p}\in\mathring{\mathcal{V}}_{0}(\hat{T})$ with
$\|\hat{p}\|_{L^{2}(\hat{T})}=1$, let $p:T\to\mathbb{R}$ be given by
$p(x)=\hat{p}(\hat{x})$. Then $p\in\mathring{\mathcal{V}}_{0}(T)$,
$\|p\|_{L^{2}(T)}=\sqrt{6|T|}$, and

$\int_{\hat{T}}(\widehat{\rm div}\,\hat{\bm{v}})\hat{p}=\int_{T}({\rm
div}\,\bm{v})p.$

We conclude

$\displaystyle\mathop{\sup_{\hat{p}\in\mathring{\mathcal{V}}_{0}(\hat{T})}}_{\|\hat{p}\|_{L^{2}(\hat{T})}=1}\Big{|}\int_{\hat{T}}(\widehat{\rm
div}\,\hat{\bm{v}})\hat{p}\Big{|}^{2}$
$\displaystyle\leq\mathop{\sup_{p\in\mathring{\mathcal{V}}_{0}(T)}}_{\|p\|_{L^{2}(T)}=\sqrt{6|T|}}\Big{|}\int_{T}({\rm
div}\,\bm{v})p\Big{|}^{2}\leq Ch_{T}^{3}\|{\rm div}\,\bm{v}\|_{L^{2}(T)}^{2}.$

Finally, we use $\|A_{T}^{-1}\|\leq Ch_{T}^{-1}$ and $|\det(A_{T})|=6|T|\leq
Ch_{T}^{3}$ to get

$\displaystyle|\hat{\bm{v}}|_{H^{m}(\hat{T})}^{2}$ $\displaystyle\leq
C\Big{(}h_{T}^{4}\sum_{a\in\Delta_{0}(T)}|\bm{v}(a)|^{2}+\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}\bm{v}\cdot\bm{n}_{F}\Big{|}^{2}$
$\displaystyle\qquad+h_{T}^{4}\sum_{F\in\Delta_{2}(T)}\sum_{e\in\Delta^{I}_{1}(F^{\rm
CT})\backslash\\{e_{F}\\}}\Big{|}\int_{e}\left[\left[{{\rm
div}\,}\bm{v}\right]\right]_{e}\Big{|}^{2}+h_{T}^{3}\|{\rm
div}\,\bm{v}\|_{L^{2}(T)}^{2}\Big{)},$

and therefore

$\displaystyle|\bm{v}|_{H^{m}(\hat{T})}^{2}\leq
Ch_{T}^{-1-2m}|\hat{\bm{v}}|_{H^{m}(\hat{T})}^{2}$ $\displaystyle\leq
Ch_{T}^{-1-2m}\Big{(}h_{T}^{4}\sum_{a\in\Delta_{0}(T)}|\bm{v}(a)|^{2}+\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}\bm{v}\cdot\bm{n}_{F}\Big{|}^{2}$
$\displaystyle\qquad+h_{T}^{4}\sum_{F\in\Delta_{2}(T)}\sum_{e\in\Delta_{1}(F)\backslash\\{e_{F}\\}}\Big{|}\int_{e}\left[\left[{{\rm
div}\,}\bm{v}\right]\right]_{e}\Big{|}^{2}+h_{T}^{3}\|{\rm
div}\,\bm{v}\|_{L^{2}(T)}^{2}\Big{)}$ $\displaystyle\leq
Ch_{T}^{-1-2m}\Big{(}h_{T}^{4}\sum_{a\in\Delta_{0}(T)}|\bm{v}(a)|^{2}+\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}\bm{v}\cdot\bm{n}_{F}\Big{|}^{2}+h_{T}^{3}\|{\rm
div}\,\bm{v}\|_{L^{2}(T)}^{2}\Big{)},$

where the last inequality comes from standard trace and inverse inequalities.
∎

######  Theorem 4.7.

The pair $\mathring{\bm{V}}_{h}^{\rm WF}\times\mathring{Y}_{h}^{\rm WF}$ is
inf-sup stable.

###### Proof.

Fix a $q\in\mathring{Y}_{h}^{\rm WF}$, and let
$\bm{w}\in\bm{H}^{1}_{0}(\Omega)$ satisfy ${\rm div}\,\bm{w}=q$ and
$\|\nabla\bm{w}\|_{L^{2}(\Omega)}\leq C\|q\|_{L^{2}(\Omega)}$. Let
$\bm{w}_{h}\in\mathring{\bm{\EuScript{P}}}_{1}(\mathcal{T}_{h})$ be the Scott-
Zhang interpolant of $\bm{w}$ with respect to $\mathcal{T}_{h}$. Define
$\bm{v}\in\mathring{\bm{\EuScript{P}}}_{1}(\mathcal{T}_{h}^{\rm WF})$ such
that

$\displaystyle\bm{v}(a)=\bm{w}_{h}(a)\qquad$ $\displaystyle\forall
a\in\Delta_{0}(T),$
$\displaystyle\int_{F}(\bm{v}\cdot\bm{n}_{F})=\int_{F}(\bm{w}\cdot\bm{n}_{F})\qquad$
$\displaystyle\forall F\in\Delta_{2}(T),$
$\displaystyle\int_{e}\left[\left[{\rm
div}\,\bm{v}\right]\right]_{e}=\int_{e}\left[\left[q\right]\right]_{e}\qquad$
$\displaystyle\forall e\in\Delta_{1}^{I}(F^{\rm CT})\backslash\\{e_{F}\\},\
\forall F\in\Delta_{2}(T),$ $\displaystyle\int_{T}({\rm
div}\,\bm{v})p=\int_{T}qp\qquad$ $\displaystyle\forall
p\in\mathring{\mathcal{V}}_{0}(T).$

Noting $({\rm div}\,\bm{v}-q)\in{\bm{Y}}_{h}^{\rm WF}$, and

$\displaystyle\int_{e}\left[\left[{\rm
div}\,\bm{v}-q\right]\right]_{e}=0\qquad$ $\displaystyle\forall
e\in\Delta_{1}^{I}(F^{\rm CT})\backslash\\{e_{F}\\},\ \forall
F\in\Delta_{2}(T),$ $\displaystyle\int_{T}({\rm div}\,\bm{v}-q)p\qquad$
$\displaystyle\forall p\in\EuScript{P}_{0}(T^{\rm A})$

for all $T\in\mathcal{T}_{h}$ (by the divergence theorem), we conclude ${\rm
div}\,\bm{v}=q$ by Lemma 4.5.

We apply Proposition 4.6 to $(\bm{v}-\bm{w}_{h})$ with $m=1$:

$\displaystyle\|\nabla(\bm{v}-\bm{w}_{h})\|_{L^{2}(T)}^{2}$ $\displaystyle\leq
Ch_{T}^{-3}\Big{(}h_{T}^{4}\sum_{a\in\Delta_{0}(T)}|(\bm{v}-\bm{w}_{h})(a)|^{2}+\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}(\bm{v}-\bm{w}_{h})\cdot\bm{n}_{F}\Big{|}^{2}$
$\displaystyle\qquad+h_{T}^{3}\|{\rm
div}\,(\bm{v}-\bm{w}_{h})\|_{L^{2}(T)}^{2}\Big{)}$
$\displaystyle=Ch_{T}^{-3}\Big{(}\sum_{F\in\Delta_{2}(T)}\Big{|}\int_{F}(\bm{w}-\bm{w}_{h})\cdot\bm{n}_{F}\Big{|}^{2}+h_{T}^{3}\|q-{\rm
div}\,\bm{w}_{h}\|_{L^{2}(T)}^{2}\Big{)}$ $\displaystyle\leq
Ch_{T}^{-3}\Big{(}h_{T}^{2}\sum_{F\in\Delta_{2}(T)}\|\bm{w}-\bm{w}_{h}\|_{L^{2}(F)}^{2}+h_{T}^{3}\|q-{\rm
div}\,\bm{w}_{h}\|_{L^{2}(T)}^{2}\Big{)}$ $\displaystyle\leq
C\Big{(}\|q\|_{L^{2}(T)}^{2}+h_{T}^{-2}\|\bm{w}-\bm{w}_{h}\|_{L^{2}(T)}^{2}+\|\nabla(\bm{w}-\bm{w}_{h})\|_{L^{2}(T)}^{2}+\|\nabla\bm{w}_{h}\|_{L^{2}(T)}^{2}\Big{)}.$

We then sum over $T\in\mathcal{T}_{h}$ and apply stability and approximation
properties of the Scott-Zhang interpolant to conclude
$\|\nabla\bm{v}\|_{L^{2}(\Omega)}\leq C\|q\|_{L^{2}(\Omega)}$.

∎

##  5\. Implementation Aspects for Powell-Sabin splits

The only tricky part to implement these finite elements is the pressure spaces
since they have non-standard constraints in their definitions. In this section
and the subsequent one, we give details to form the algebraic system for the
Stokes problem.

###  5.1. A basis for $\hat{Y}_{h}^{\rm PS}$ and the construction of the
algebraic system

Recall that the space $\hat{Y}_{h}^{\rm PS}$ consists of piecewise constant
functions satisfying a weak continuity property at vertices. Clearly, this is
a non-standard space, and in particular the space is not explicitly found in
current finite element software packages. Nonetheless, in this section, we
identify a basis of the space $\mathring{Y}_{h}^{\rm PS}$, and as a byproduct
show that the weak continuity property $\theta_{z}(q)=0$ can be imposed purely
at the algebraic level.

As a first step, we note that, by definition of the Powell–Sabin
triangulation,

$\mathcal{T}_{h}^{\rm PS}=\\{K_{z}^{(j)}:\ K_{z}^{(j)}\in\mathcal{T}_{z},\
z\in\mathcal{S}_{h}\\}.$

With this (implicit) labeling of the triangles in $\mathcal{T}_{h}^{\rm PS}$,
we can write the canonical basis of $\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm
PS})$ as the set
$\\{\varphi_{z}^{(j)}\\}\subset\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm PS})$
with

$\varphi_{z}^{(j)}|_{K_{v}^{(i)}}=\delta_{v,z}\delta_{i,j}\qquad\forall
z,v\in\mathcal{S}_{h},\ i=1,\ldots,n_{v},\ j=1,\ldots,n_{z}.$

The next proposition shows that a basis of $\mathring{Y}_{h}^{\rm PS}$ is
easily extracted from the basis of $\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm
PS})$ (see Figure 2).

######  Proposition 5.1.

For each $z\in\mathcal{S}_{h}$ and $j\in\\{2,\ldots,n_{z}\\}$, define

$\psi_{z}^{(j)}=\varphi_{z}^{(j)}+(-1)^{j}\varphi_{z}^{(1)}.$

Then $\\{\psi_{z}^{(j)}:\ z\in\mathcal{S}_{h},\ j=2,\ldots,n_{z}\\}$ forms a
basis of $\hat{Y}_{h}^{\rm PS}$.

###### Proof.

Note that the number of functions $\\{\psi_{z}^{(j)}\\}$ given is
$\sum_{z\in\mathcal{S}_{h}}(n_{z}-1)$, and

$\dim Y_{h}^{PS}=\dim\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm
PS})-|\mathcal{S}_{h}|=\sum_{z\in\mathcal{S}_{h}}n_{z}-|\mathcal{S}_{h}|=\sum_{z\in\mathcal{S}_{h}}(n_{z}-1).$

Because $\psi_{z}^{(j)}\in\mathring{Y}_{h}^{\rm PS}$, and they are clearly
linear independent, we conclude that $\\{\psi_{z}^{(j)}\\}$ form a basis of
$\mathring{Y}_{h}^{\rm PS}$. ∎

$1$$0$$0$$0$

$0$$1$$0$$0$

$0$$0$$0$$1$

$0$$0$$1$$0$

$1$$1$$0$$0$

$-1$$0$$0$$1$

$1$$0$$1$$0$

Figure 2. Local mesh $\mathcal{T}_{z}$ with $z\in\mathcal{S}_{h}^{I}$. Top
row: Values of canonical basis functions of piecewise constants
$\\{\varphi_{z}^{(j)}\\}_{j=1}^{n_{z}}$. Bottom row: Values of basis functions
of piecewise constants with weak continuity constraint
$\\{\psi_{z}^{(j)}\\}_{j=2}^{n_{z}}$.

We now explain how Proposition 5.1 provides a simple way to construct the
stiffness matrix for the Stokes problem using the $\mathring{\bm{V}}_{h}^{\rm
PS}\times\mathring{Y}_{h}^{\rm PS}$ pair. To explain the procedure, we require
some notation.

Let $A$ be the matrix associated with the bilinear form

$(\bm{v},\bm{w})\to\int_{\Omega}\nu\nabla\bm{v}:\nabla\bm{w}\,dx\quad\text{over
}\bm{v},\bm{w}\in\mathring{\bm{V}}_{h}^{\rm PS},$

and let $\tilde{B}$ is the matrix associated with the bilinear form

$(\bm{v},q)\to-\int_{\Omega}({\rm div}\,\bm{v})q\,dx\quad\text{over
}\bm{v}\in\mathring{\bm{V}}_{h}^{\rm PS},\
q\in\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm PS}).$

The stiffness matrix for the Stokes problem based on the (unstable)
$\mathring{\bm{V}}_{h}^{\rm PS}\times\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm
PS})$ pair is given by

$\displaystyle\begin{pmatrix}A&\tilde{B}\\\
\tilde{B}^{\intercal}&0\end{pmatrix}.$

We emphasize that this system can be easily constructed using standard finite
element software packages.

Let $\\{\bm{\phi}^{(i)}\\}_{k=1}^{N}$ denote a basis of
$\mathring{\bm{V}}_{h}^{\rm PS}$ with $N=\dim\mathring{\bm{V}}_{h}^{\rm PS}$
so that

$A_{i,j}=\nu\int_{\Omega}\nabla\bm{\phi}^{(j)}:\nabla\bm{\phi}^{(i)}\,dx.$

Let $M=\dim\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm PS})$, the number of
triangles in $\mathcal{T}_{h}^{\rm PS}$, and introduce the local-to-global
label mapping
$\sigma:\mathcal{S}\times\\{1,\ldots,n_{z}\\}\to\\{1,2,\ldots,M\\}$ such that

$\tilde{B}_{i,\sigma(z,j)}=-\int_{\Omega}({\rm
div}\,\bm{\phi}^{(i)})\varphi_{z}^{(j)}\,dx.$

Then by Proposition 5.1, we have for $z\in\mathcal{S}_{h}$ and
$j=2,\ldots,n_{z}$,

$\displaystyle-\int_{\Omega}({\rm div}\,\bm{\phi}^{(i)})\psi_{z}^{(j)}\,dx$
$\displaystyle=-\int_{\Omega}({\rm
div}\,\bm{\phi}^{(i)})\varphi_{z}^{(j)}\,dx-(-1)^{j}\int_{\Omega}({\rm
div}\,\bm{\phi}^{(i)})\varphi_{z}^{(1)}\,dx$
$\displaystyle=\tilde{B}_{i,\sigma(z,j)}+(-1)^{j}\tilde{B}_{i,\sigma(z,1)}.$

This identity leads to the following algorithm.

Algorithm

  1. (1)

Construct Powell–Sabin triangulation $\mathcal{T}_{h}^{\rm PS}$

  2. (2)

Construct $\tilde{B}\in\mathbb{R}^{N\times M}$ based on the
$\mathring{\bm{V}}_{h}^{\rm PS}\times\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm
PS})$ pair.

  3. (3)

Set $B=\tilde{B}$.

  4. (4)

For each $z\in\mathcal{S}_{h}$ and for each $j\in\\{2,\ldots,n_{z}\\}$, do the
elementary column operation

$B_{:,\sigma(z,j)}=B_{:,\sigma(z,j)}+(-1)^{j}B_{:,\sigma(z,1)}.$

  5. (5)

Delete column $B_{:,\sigma(z,1)}$ for each $z\in\mathcal{S}_{h}$.

The stiffness matrix for the Stokes problem based on the
$\mathring{\bm{V}}_{h}^{\rm PS}\times\hat{Y}_{h}^{\rm PS}$ pair is then given
by

(5.1) $\displaystyle\begin{pmatrix}A&B\\\ B^{\intercal}&0\end{pmatrix}.$

##  6\. Implemenation aspects for Worsey Farin Splits

###  6.1. A basis for $\hat{Y}_{h}^{\rm WF}$ and the construction of the
algebraic system

Notice that the collection of local triangulations $\mathcal{T}_{e}$ (with
$e\in\mathcal{E}^{\mathcal{S}}_{h}$) do not form a disjoint partition of the
global triangulation $\mathcal{T}_{h}^{\rm WF}$. In particular, there exists
$K\in\mathcal{T}_{h}^{\rm WF}$ such that $K\in\mathcal{T}_{e}$ and
$K\in\mathcal{T}_{s}$ with $e,s\in\mathcal{E}^{\mathcal{S}}_{h}$ and $e\neq
s$. As a result, the methodology used in the previous section for Powell–Sabin
meshes is not directly applicable.

Instead, we consider a geometric decomposition of the mesh based on the face
split points in $\mathcal{T}_{h}^{\rm WF}$. To this end, we denote by
$\mathcal{S}^{I}_{h}$ and $\mathcal{S}_{h}^{B}$ the sets of interior and
boundary face split points, respectively, and set
$\mathcal{S}_{h}=\mathcal{S}^{I}_{h}\cup\mathcal{S}_{h}^{B}$. For
$z\in\mathcal{S}_{h}$, let
$\mathcal{T}_{z}:=\\{K_{z}^{(1)},\ldots,K_{z}^{(n_{z})}\\}$ denote the set of
tetrahedra in $\mathcal{T}_{h}^{\rm WF}$ that have $z$ as a vertex. Here,
$n_{z}=6$ if $z$ is an interior vertex, and $n_{z}=3$ if $z$ is a boundary
vertex. For an interior face split point $z$, we assume the simplices in
$\mathcal{T}_{z}$ are labeled such that

$K_{z}^{(1)},K_{z}^{(2)},K_{z}^{(3)}\subset T^{(1)},\qquad
K_{z}^{(4)},K_{z}^{(5)},K_{z}^{(6)}\subset T^{(2)}$

for some $T^{(1)},T^{(2)}\in\mathcal{T}_{h}$, and that $K_{z}^{(j)}$ and
$K_{z}^{(j+3)}$ share a common face for $j\in\\{1,2,3\\}$. For a boundary
split point $z$, the set
$\mathcal{T}_{z}=\\{K_{z}^{(1)},K_{z}^{(2)},K_{z}^{(3)}\\}$ is labeled
arbitrarily.

We clearly have

(6.1) $\mathcal{T}_{h}^{\rm WF}=\\{K_{z}^{(j)}:\
z\in\mathcal{S}_{h},j=1,\ldots,n_{z}\\},$

and the map $(z,j)\to K_{z}^{(j)}$ is injective. Furthermore, each local
partition $\mathcal{T}_{z}$ contains three singular edges.

######  Proposition 6.1.

There holds

$\dim\hat{Y}_{h}^{\rm WF}=4|\mathcal{S}_{h}^{I}|+|\mathcal{S}_{h}^{B}|.$

###### Proof.

By Proposition 4.5, we have

$\dim\hat{Y}_{h}^{\rm WF}=4|\mathcal{T}_{h}|+2|\mathcal{F}^{I}_{h}|,$

where $|\mathcal{F}^{I}_{h}|$ is the number of interior faces in
$\mathcal{T}_{h}$. From (6.1), we have

$12|\mathcal{T}_{h}|=|\mathcal{T}_{h}^{\rm
WF}|=6|\mathcal{S}_{h}^{I}|+3|\mathcal{S}_{h}^{B}|,$

and by construction of the Worsey-Farin split, there holds

$|\mathcal{F}_{h}|=|\mathcal{S}^{I}_{h}|.$

Therefore,

$\displaystyle\dim\hat{Y}_{h}^{\rm WF}$
$\displaystyle=4|\mathcal{T}_{h}|+2|\mathcal{F}^{I}_{h}|=\frac{1}{3}\big{(}6|\mathcal{S}_{h}^{I}|+3|\mathcal{S}_{h}^{B}|\big{)}+2|\mathcal{S}_{h}^{I}|=4|\mathcal{S}_{h}^{I}|+|\mathcal{S}_{h}^{B}|.$

∎

For an interior split point $z$, and for a piecewise constant function $q$ on
$\mathcal{T}_{z}$, the three constraints $\theta_{e}(q)=0$ read

$\displaystyle q_{1}-q_{2}+q_{5}-q_{4}$ $\displaystyle=0,$ $\displaystyle
q_{2}-q_{3}+q_{6}-q_{5}$ $\displaystyle=0,$ $\displaystyle
q_{3}-q_{1}+q_{4}-q_{6}$ $\displaystyle=0,$

where $q_{j}=q|_{K^{(j)}_{z}}$

We write this as a $3\times 6$ linear system

$\displaystyle C\vec{q}:=\begin{pmatrix}1&-1&0&-1&1&0\\\ 0&1&-1&0&-1&1\\\
-1&0&1&1&0&-1\end{pmatrix}\begin{pmatrix}q_{1}\\\ q_{2}\\\ q_{3}\\\ q_{4}\\\
q_{5}\\\ q_{6}\end{pmatrix}=0.$

We clearly see that this matrix has rank $2$ (e.g., adding the first and third
rows gets the negation of the second row). We find that the nullspace of $C$
is given by

$\displaystyle{\rm null}(C)={\rm span}\left\\{\begin{pmatrix}1\\\ 1\\\ 1\\\
0\\\ 0\\\ 0\end{pmatrix},\begin{pmatrix}1\\\ 0\\\ 0\\\ 1\\\ 0\\\
0\end{pmatrix},\begin{pmatrix}0\\\ 1\\\ 0\\\ 0\\\ 1\\\
0\end{pmatrix},\begin{pmatrix}-1\\\ -1\\\ 0\\\ 0\\\ 0\\\
1\end{pmatrix}\right\\}.$

These four vectors implicitly give us a basis for $\hat{Y}_{h}^{\rm WF}$. In
particular, we have

######  Proposition 6.2.

For $z\in\mathcal{S}_{h}$ and $j\in\\{1,2,\ldots,n_{z}\\}$, let
$\varphi_{z}^{(j)}$ be the piecewise constant function

$\varphi_{z}^{(j)}|_{K_{v}^{(i)}}=\delta_{v,z}\delta_{i,j}\qquad\forall
v,z\in\mathcal{S}_{h},\ i=1,\ldots,n_{v},\ j=1,\ldots,n_{z}.$

For an interior face split point $z$, define

$\displaystyle\psi_{z}^{(3)}$
$\displaystyle=\varphi_{z}^{(3)}+\varphi_{z}^{(1)}+\varphi_{z}^{(2)},$
$\displaystyle\psi_{z}^{(4)}$
$\displaystyle=\varphi_{z}^{(4)}+\varphi_{z}^{(1)},$
$\displaystyle\psi_{z}^{(5)}$
$\displaystyle=\varphi_{z}^{(5)}+\varphi_{z}^{(2)},$
$\displaystyle\psi_{z}^{(6)}$
$\displaystyle=\varphi_{z}^{(6)}-\varphi_{z}^{(1)}-\varphi_{z}^{(2)}.$

For a boundary face split point $z$, define

$\displaystyle\psi_{z}^{(3)}$
$\displaystyle=\varphi_{z}^{(3)}+\varphi_{z}^{(1)}+\varphi_{z}^{(2)}.$

Then $\\{\psi_{z}^{(j)}\\}$ is a basis of $\hat{Y}_{h}^{\rm WF}$.

###### Proof.

The proof essentially follows from the same arguments as Proposition 5.1,
noting that the number of given $\psi_{z}^{(j)}$ is

$\displaystyle
4|\mathcal{S}_{h}^{I}|+|\mathcal{S}_{h}^{B}|=\dim\hat{Y}_{h}^{\rm WF}$

by Proposition 6.1. ∎

As in the two-dimensional case, Proposition 6.2 give an algorithm to construct
the stiffness matrix for the Stokes problem using the
$\mathring{\bm{V}}_{h}^{\rm WF}\times\hat{Y}_{h}^{\rm WF}$ pair. First, we
construct the stiffness matrix based on the $\mathring{\bm{V}}_{h}^{\rm
WF}\times\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm WF})$ pair:

$\displaystyle\begin{pmatrix}A&\tilde{B}\\\
\tilde{B}^{\intercal}&0\end{pmatrix},$

and then perform elementary column operations on the $\tilde{B}$.

Let $\\{\bm{\phi}^{(k)}\\}_{k=1}^{N}$ denote a basis of
$\mathring{\bm{V}}_{h}^{\rm WF}$ with $N=\dim\bm{V}_{h}^{\rm PS}$ Let
$M=\dim\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm PS})$, the number of tetrahedra
in $\mathcal{T}_{h}^{\rm WF}$, and introduce the local-to-global label mapping
$\sigma:\mathcal{S}_{h}\times\\{1,\ldots,n_{z}\\}\to\\{1,2,\ldots,M\\}$ such
that

$\tilde{B}_{k,\sigma(z,j)}=-\int_{\Omega}({\rm
div}\,\bm{\phi}^{(k)})\varphi_{z}^{(j)}.$

Algorithm 1

  1. (1)

Construct Worsey-Farin triangulation $\mathcal{T}_{h}^{\rm PS}$

  2. (2)

Construct $\tilde{B}\in\mathbb{R}^{N\times M}$ based on the
$\mathring{\bm{V}}_{h}^{\rm WF}\times\EuScript{P}_{0}(\mathcal{T}_{h}^{\rm
WF})$ pair.

  3. (3)

Set $B=\tilde{B}$.

  4. (4)

For each interior face split point $z\in\mathcal{S}_{h}$ do the elementary
column operations

$\displaystyle B_{:,\sigma(z,3)}$
$\displaystyle=B_{:,\sigma(z,3)}+B_{:,\sigma(z,1)}+B_{:,\sigma(z,2)},$
$\displaystyle B_{:,\sigma(z,4)}$
$\displaystyle=B_{:,\sigma(z,4)}+B_{:,\sigma(z,1)},$ $\displaystyle
B_{:,\sigma(z,5)}$ $\displaystyle=B_{:,\sigma(z,5)}+B_{:,\sigma(z,2)},$
$\displaystyle B_{:,\sigma(z,6)}$
$\displaystyle=B_{:,\sigma(z,6)}-B_{:,\sigma(z,1)}-B_{:,\sigma(z,2)}.$

  5. (5)

For each boundary face split point $z\in\mathcal{S}_{h}$ do the elementary
column operation

$\displaystyle B_{:,\sigma(z,3)}$
$\displaystyle=B_{:,\sigma(z,3)}+B_{:,\sigma(z,1)}+B_{:,\sigma(z,2)}.$

  6. (6)

Delete columns $B_{:,\sigma(z,1)}$ and $B_{:,\sigma(z,2)}$ for each
$z\in\mathcal{S}_{h}$.

The stiffness matrix for the Stokes problem based on the
$\mathring{\bm{V}}_{h}^{\rm WF}\times\hat{Y}_{h}^{\rm WF}$ pair is then given
by

(6.2) $\displaystyle\begin{pmatrix}A&B\\\ B^{\intercal}&0\end{pmatrix},$

##  7\. Numerical Experiments

In this section, we perform some simple numerical experiments for the Stokes
problem on Powell–Sabin and Worsey–Farin splits. We note standard theory shows
that the velocity and pressure errors satisfy

(7.1) $\displaystyle|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$
$\displaystyle\leq(1+\beta^{-1})\inf_{\bm{v}_{h}\in\bm{V}_{h}}|\bm{v}_{h}-\bm{u}|_{H^{1}(\Omega)},$
(7.2) $\displaystyle\|p-p_{h}\|_{L^{2}(\Omega)}$ $\displaystyle\leq\inf_{q\in
Y_{h}}\|p-q\|_{L^{2}(\Omega)}+\frac{\nu}{\beta}|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)},$

where either $\bm{V}_{h}\times Y_{h}=\bm{V}_{h}^{\rm
PS}\times\mathring{Y}_{h}^{\rm PS}$ or $\bm{V}_{h}\times
Y_{h}=\mathring{\bm{V}}_{h}^{\rm WF}\times\mathring{Y}_{h}^{\rm WF}$, $\nu>0$
is the viscosity, and $\beta$ is the inf-sup constant for the finite element
pair $\bm{V}_{h}\times Y_{h}$.

###  7.1. The Stokes pair on Powell-Sabin Splits

We consider the example such that the data is taken to be $\Omega=(0,1)^{2}$,
and the source function is chosen such that the exact velocity and pressure
solutions for (2.1) are given respectively as

(7.3) $\displaystyle\bm{u}=\begin{pmatrix}\pi\sin^{2}(\pi x_{1})\sin(2\pi
x_{2})\\\ -\pi\sin^{2}(\pi x_{2})\sin(2\pi x_{1})\end{pmatrix},\quad
p=\cos(\pi x_{1})\cos(\pi x_{2}).$

Let $\mathcal{T}_{h}$ be a Delaunay triangulation of $\Omega$ and
$\mathcal{T}_{h}^{PS}$ the corresponding Powell-Sabin global triangulation.

The resulting errors, rates of convergence, and inf-sup constants are listed
in Tables 1 and 2 for viscosities $\nu=1$ and $\nu=10^{-2}$, respectively. The
results show that the $L^{2}$ pressure error and the $H^{1}$ velocity error
converge with linear rate, the discrete velocity solution (and error) are
independent of the viscosity $\nu$, and the pressure error improves for small
viscosity. The experiments also show that the inf-sup constant does not
deteriorate as the mesh is refined with $\beta\approx 0.1$. These results are
in agreement with the theoretical estimates (7.1)–(7.2)

In Tables 3 and 4, we compute the right-hand side of (7.2) and (7.1),
respectively, and compare the data with the computed errors
$\|p-p_{h}\|_{L^{2}(\Omega)}$ and $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$.
Again, the results are consistent with (7.1)–(7.2), and they suggest that the
term $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$ is the dominant term in the
pressure error (7.2).

Table 1. Errors and rates of convergence for example (7.3) with $\nu=1$. $h$ | $\|\bm{u}-\bm{u}_{h}\|_{L^{2}(\Omega)}$ | rate | $\|p-p_{h}\|_{L^{2}(\Omega)}$ | rate | $\|\nabla\cdot\bm{u}_{h}\|_{L^{2}(\Omega)}$ | $\beta$
---|---|---|---|---|---|---
$2^{-2}$ | 1.70E-01 | – | 5.26E 00 | – | 2.70E-14 | 1.56E-01
$2^{-3}$ | 5.66E-02 | 1.587 | 3.77E 00 | 0.480 | 6.65E-14 | 1.38E-01
$2^{-4}$ | 1.35E-02 | 2.068 | 1.68E 00 | 1.166 | 2.38E-13 | 1.07E-01
$2^{-5}$ | 3.35E-03 | 2.011 | 8.28E-01 | 1.021 | 8.38E-12 | 1.06E-01
$2^{-6}$ | 8.77E-04 | 1.934 | 4.25E-01 | 0.962 | 4.05E-10 | 9.34E-02
Table 2. Errors and rates of convergence for example (7.3) with $\nu=10^{-2}$. $h$ | $\|\bm{u}-\bm{u}_{h}\|_{L^{2}(\Omega)}$ | rate | $\|p-p_{h}\|_{L^{2}(\Omega)}$ | rate | $\|\nabla\cdot\bm{u}_{h}\|_{L^{2}(\Omega)}$
---|---|---|---|---|---
$2^{-2}$ | 1.70E-01 | – | 1.02E-01 | – | 2.43E-14
$2^{-3}$ | 5.66E-02 | 1.587 | 5.79E-02 | 0.816 | 5.88E-14
$2^{-4}$ | 1.35E-02 | 2.068 | 2.76E-02 | 1.069 | 2.36E-13
$2^{-5}$ | 3.35E-03 | 2.011 | 1.37E-02 | 1.010 | 8.39E-12
$2^{-6}$ | 8.77E-04 | 1.934 | 6.96E-03 | 0.977 | 4.05E-10
Table 3. Errors for example (7.3) with $\nu=10^{-2}$ and the RHS of (7.2) . $h$ | $\|p-p_{h}\|_{L^{2}(\Omega)}$ | $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$ | $\beta$ | $\inf_{q\in\mathring{Y}_{h}^{\rm PS}}\|p-q\|_{L^{2}(\Omega)}$ | RHS of (7.2)
---|---|---|---|---|---
$2^{-2}$ | 1.02E-01 | 3.77E 00 | 1.56E-01 | 6.08E-02 | 3.02E-01
$2^{-3}$ | 5.79E-02 | 2.17E 00 | 1.38E-01 | 2.77E-02 | 1.84E-01
$2^{-4}$ | 2.76E-02 | 1.07E 00 | 1.07E-01 | 1.35E-02 | 1.13E-01
$2^{-5}$ | 1.37E-02 | 5.32E-01 | 1.06E-01 | 6.61E-03 | 5.67E-02
$2^{-6}$ | 6.96E-03 | 2.72E-01 | 9.34E-02 | 3.28E-03 | 3.24E-02
Table 4. Errors for example (7.3) with $\nu=10^{-2}$ and the RHS of (7.1) . $h$ | $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$ | $\beta$ | $\inf_{\bm{v}_{h}\in\bm{V}_{h}^{\rm PS}}|\bm{v}_{h}-\bm{u}|_{H^{1}(\Omega)}$ | RHS of (7.1)
---|---|---|---|---
$2^{-2}$ | 3.77E 00 | 1.56E-01 | 3.08E 00 | 2.28E+01
$2^{-3}$ | 2.17E 00 | 1.38E-01 | 1.63E 00 | 1.34E+01
$2^{-4}$ | 1.07E 00 | 1.07E-01 | 8.04E-01 | 8.31E 00
$2^{-5}$ | 5.32E-01 | 1.06E-01 | 4.06E-01 | 4.23E 00
$2^{-6}$ | 2.72E-01 | 9.34E-02 | 2.05E-01 | 2.39E 00

###  7.2. The Stokes pair on Worsey-Farin Splits

We consider the example such that the data is taken to be $\Omega=(0,1)^{3}$,
and the source function is chosen such that the exact velocity and pressure
solutions for (2.1) are given respectively as

(7.4) $\displaystyle\bm{u}=\begin{pmatrix}\pi\sin^{2}(\pi x_{1})\sin(2\pi
x_{2})\\\ -\pi\sin^{2}(\pi x_{2})\sin(2\pi x_{1})\\\ 0\end{pmatrix},\quad
p=\cos(\pi x_{1})\cos(\pi x_{2})\cos(\pi x_{3}).$

Let $\mathcal{T}_{h}$ be a Delaunay triangulation of $\Omega$ and
$\mathcal{T}_{h}^{WF}$ be the corresponding Worsey-Farin global triangulation.

The resulting rates of convergence of the numerical experiments for
viscosities $\nu=1$ and $\nu=10^{-3}$ are listed in Tables 5 and 6,
respectively. We also state the computed inf-sup constant on these meshes, and
the results show that that it stays uniformly bounded from below with
$\beta\approx 0.13$ on all meshes. The stated errors, especially those in
Table 5, indicate that the rates of convergence are still in the preasymptotic
regime. On the other hand, for small viscosity value $\nu=10^{-3}$, Table 6
shows that the pressure error converges with linear rate. This behavior
suggests that the velocity error is the dominating term in (7.2).

To verify this claim, we explicitly compute the right-hand side of (7.2) and
(7.1) and report the results in Tables 7 and 8, respectively. The results show
that indeed $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$ is the dominating term in
the pressure error (7.2).

Table 5. Errors and rates of convergence for example (7.4) with $\nu=1$. $h$ | $\|\bm{u}-\bm{u}_{h}\|_{L^{2}(\Omega)}$ | rate | $\|p-p_{h}\|_{L^{2}(\Omega)}$ | rate | $\|\nabla\cdot\bm{u}_{h}\|_{L^{2}(\Omega)}$ | $\beta$
---|---|---|---|---|---|---
$1/2$ | 1.29E 00 | – | 9.81E 00 | – | 5.07E-14 | 1.31E-01
$1/4$ | 8.58E-01 | 0.588 | 19.4E 00 | -0.98 | 5.20E-13 | 1.31E-01
$1/8$ | 3.93E-01 | 1.286 | 16.6E 00 | 0.414 | 2.68E-12 | 1.32E-01
$1/16$ | 1.32E-01 | 1.573 | 10.5E 00 | 0.667 | 4.10E-12 | 1.32E-01
$1/32$ | 3.69E-02 | 1.839 | 5.75E 00 | 0.872 | 4.32E-12 | 1.32E-01
$1/48$ | 1.68E-02 | 1.941 | 3.93E 00 | 0.936 | 6.07E-12 | 1.32E-01
Table 6. Errors and rates of convergence for example (7.4) with $\nu=10^{-3}$. $h$ | $\|\bm{u}-\bm{u}_{h}\|_{L^{2}(\Omega)}$ | rate | $\|p-p_{h}\|_{L^{2}(\Omega)}$ | rate | $\|\nabla\cdot\bm{u}_{h}\|_{L^{2}(\Omega)}$
---|---|---|---|---|---
$1/2$ | 1.29E 00 | – | 1.33E-01 | – | 1.28E-15
$1/4$ | 8.58E-01 | 0.588 | 6.97E-02 | 0.932 | 3.43E-14
$1/8$ | 3.93E-01 | 1.286 | 3.70E-02 | 0.911 | 3.22E-13
$1/16$ | 1.32E-01 | 1.574 | 1.91E-02 | 0.953 | 6.40E-13
$1/32$ | 3.69E-02 | 1.838 | 9.68E-03 | 0.980 | 9.52E-13
$1/48$ | 9.63E-03 | 1.940 | 4.89E-03 | 0.983 | 1.03E-12
Table 7. Errors for example (7.4) with $\nu=10^{-3}$ and the RHS of (7.2) . $h$ | $\|p-p_{h}\|_{L^{2}(\Omega)}$ | $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$ | $\beta$ | $\inf_{q\in\mathring{Y}_{h}^{\rm WF}}\|p-q\|_{L^{2}(\Omega)}$ | RHS of (7.2)
---|---|---|---|---|---
$1/2$ | 1.33E-01 | 1.07E+01 | 1.31E-01 | 5.00E-01 | 5.81E-01
$1/4$ | 6.97E-02 | 8.50E 00 | 1.31E-01 | 6.70E-02 | 1.31E-01
$1/6$ | 4.81E-02 | 6.82E 00 | 1.32E-01 | 4.43E-02 | 9.59E-02
$1/8$ | 3.70E-02 | 5.60E 00 | 1.32E-01 | 3.30E-02 | 7.54E-02
$1/10$ | 3.02E-02 | 4.71E 00 | 1.32E-01 | 2.63E-02 | 6.19E-02
$1/12$ | 2.55E-02 | 4.06E 00 | 1.32E-01 | 2.19E-02 | 5.26E-02
Table 8. Errors for example (7.4) with $\nu=10^{-3}$ and the RHS of (7.1) . $h$ | $|\bm{u}-\bm{u}_{h}|_{H^{1}(\Omega)}$ | $\beta$ | $\inf_{\bm{v}_{h}\in\mathring{\bm{V}}_{h}^{\rm WF}}|\bm{v}_{h}-\bm{u}|_{H^{1}(\Omega)}$ | RHS of (7.1)
---|---|---|---|---
$1/2$ | 1.07E+01 | 1.31E-01 | 9.99E 00 | 8.62E+01
$1/4$ | 8.50E 00 | 1.31E-01 | 5.85E 00 | 5.05E+01
$1/6$ | 6.82E 00 | 1.32E-01 | 4.08E 00 | 3.49E+01
$1/8$ | 5.60E 00 | 1.32E-01 | 3.10E 00 | 2.65E+01
$1/10$ | 4.71E 00 | 1.32E-01 | 2.49E 00 | 2.13E+01
$1/12$ | 4.06E 00 | 1.32E-01 | 2.08E 00 | 1.78E+01

###  7.3. Iterated Penalty Method for $\bf{(P_{1},P_{0})}$ pair on Worsey-
Farin Splits

We consider the example such that the data is taken to be $\Omega=(0,1)^{3}$,
and the source function is chosen such that the exact velocity and pressure
solutions for (2.1) are given respectively as

(7.5) $\displaystyle\bm{u}(x,y,z)=\nabla\times\begin{pmatrix}0\\\ g\\\
g\end{pmatrix},\quad p=\frac{1}{9}\frac{\partial^{2}g}{\partial x\partial y},$

where

$g=g(x,y,z)=2^{12}(x-x^{2})^{2}(y-y^{2})^{2}(z-z^{2})^{2}.$

Similar to the previous section, we let $\mathcal{T}_{h}$ be a Delaunay
triangulation of $\Omega$ and $\mathcal{T}_{h}^{\rm WF}$ be the corresponding
Worsey-Farin global triangulation.

The iterated penalty method [3] applied to the Stokes equations with
$\mathring{\bm{V}}_{h}=\mathring{\bm{V}}_{h}^{\rm WF}$ reads : Let
$\bm{u}_{h}^{0}=\bm{0}$ and $\rho,\gamma>0$ be parameters. For $n\geq 1$,
$\bm{u}_{h}^{n}$ is recursively defined to be the solution to the variational
formulation

(7.6)
$\nu(\nabla\bm{u}_{h}^{n},\nabla{\bm{v}})+\gamma(\nabla\cdot{\bm{v}},\nabla\cdot{\bm{u}_{h}^{n}})=({{\bm{f}},{\bm{v}}})-(\sum_{i=0}^{n-1}\rho\nabla\cdot{\bm{u}_{h}^{i}},\nabla\cdot{\bm{v}}),\quad\forall\bm{v}\in\mathring{\bm{V}}_{h}^{\rm
WF}.$

It is shown in [3] that $\lim_{n\to\infty}\bm{u}_{h}^{n}=\bm{u}_{h}$ and
$\lim_{n\to\infty}\sum_{i=0}^{n}\rho\nabla\cdot\bm{u}_{h}^{i}=p_{h}$. Also,
[3] suggests to use $\|\nabla\cdot\bm{u}_{h}^{n}\|_{L^{2}(\Omega)}$ as a
stopping criterion since the difference error between
$\bm{u}_{h}^{n},\bm{u}_{h}$ is given by

$\|\bm{u}_{h}^{n}-\bm{u}_{h}\|_{L^{2}(\Omega)}\leq
C\|\nabla\cdot\bm{u}_{h}^{n}\|_{L^{2}(\Omega)}.$

The resulting rates of convergence of the numerical experiment are listed in
Tables 9. The errors $\|\bm{u}-\bm{u}_{h}^{n}\|_{L^{2}(\Omega)}$ and
$\|p-p_{h}^{n}\|_{L^{2}(\Omega)}$ are computed with
$\|\nabla\cdot\bm{u}_{h}^{n}\|_{L^{2}(\Omega)}\leq 10^{-7}$ and
$\gamma=\rho=100$.

Table 9. Errors and rates of convergence for example (7.5) with $\nu=1$. $h$ | $\|\bm{u}-\bm{u}_{h}^{n}\|_{L^{2}(\Omega)}$ | rate | $|\bm{u}-\bm{u}_{h}^{n}|_{H^{1}(\Omega)}$ | rate | $\|p-p_{h}^{n}\|_{L^{2}(\Omega)}$ | rate
---|---|---|---|---|---|---
$2^{-2}$ | 1.11768 | - | 11.55063 | - | 25.32256 | -
$2^{-3}$ | 0.48896 | 1.19273 | 7.53829 | 0.61566 | 22.35349 | 0.17992
$2^{-4}$ | 0.15482 | 1.65908 | 4.15598 | 0.85905 | 13.67635 | 0.70882
$2^{-5}$ | 0.04176 | 1.89040 | 2.13224 | 0.96282 | 7.24129 | 0.91736
$1/48$ | 0.01881 | 1.96680 | 1.42643 | 0.99145 | 4.88909 | 0.96875

In Table 10 we provide a time comparison between the (IPM) and the method
described in Algorithm 1. All of the timings were done on a machine with a
single 3.60 GHz Intel Core i9-9900K processor with 128 GB of 2400 MHz DDR4
memory.

For the iterated penalty method, we select $\gamma=\rho=100$ and terminate
iterations once $\|\nabla\cdot\bm{u}_{h}^{n}\|_{L^{2}(\Omega)}\leq 10^{-7}$.
At each iteration of IPM, equation (7.6) must be resolved. In our work, this
vector Poisson problem type problem is solved via a conjugate gradient method
with an algebraic multigrid (AMG) preconditioner. It is well known that
optimal multigrid methods can be designed for this class of problems [17].

For Algorithm 1, we instead work with the full discretization matrix (6.2),
which results in a symmetric indefinite linear system. To efficiently solve
this saddle point problem, a block preconditioned Krylov subspace method is
used [16]. In particular, we precondition the stiffness matrix (6.2) by the
block diagonal matrix

$\begin{pmatrix}A&0\\\ 0&S\end{pmatrix},$

where $S=B^{\intercal}A^{-1}B$ is the Schur complement, and the flexible GMRES
method as an outer iteration. The solver is terminated once the Euclidean norm
of the residual is less than or equal to $10^{-8}.$ For simplicity, the inner
preconditioners are as follows: we use a preconditioner to $A$ for
approximating $A^{-1}$, and the Schur complement is approximated as $S\approx
B^{\intercal}\mathrm{diag}(A^{-1})B$. It should be noted that this is not the
only choice for block-type preconditioning of the Stokes problem (e.g., see
[4]); however, our numerical experiments indicate that we can repurpose
existing preconditioners to efficiently solve linear systems that arise from
the proposed finite element discretization.

From Table 10, it is evident that for a fixed mesh spacing, the total run
times for both methods are competitive, with the IPM technique being slightly
slower for larger meshes. We find that in the case of the linear Stokes
problem, by utilizing preexisting preconditioners, both approaches provide
similar accuracy and time to solution metrics.

To obtain a Reynolds robust, optimally scaling preconditioner, special
multigrid techniques may have to be investigated [6, 5].

Table 10. Time comparison between (IPM) and Algorithm 1 for example (7.5) with $\nu=1$. $h$ | (IPM) | Algorithm 1
---|---|---
$2^{-2}$ | 5.08E+00 | 3.35E+00
$2^{-3}$ | 1.68E+01 | 2.93E+01
$2^{-4}$ | 4.80E+02 | 2.60E+02
$2^{-5}$ | 2.39E+03 | 9.37E+02
$1/48$ | 7.31E+03 | 6.45E+03

## References

  * [1] D. N. Arnold and J. Qin.  Quadratic velocity/linear pressure Stokes elements.  Advances in computer methods for partial differential equations, 7:28–34, 1992.
  * [2] D. Boffi, J. Guzman, and M. Neilan.  Convergence of Lagrange finite elements for the Maxwell eigenvalue problem in 2d.  arXiv preprint arXiv:2003.08381, 2020.
  * [3] S. C. Brenner and L. Scott.  The Mathematical Theory of Finite Elements Methods, volume 3.  Springer, 2007.
  * [4] H. Elman, V. Howle, J. Shadid, R. Shuttleworth, and R. Tuminaro.  A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations.  Journal of Computational Physics, 227(3):1790–1808, 2008.
  * [5] P. E. Farrell, M. G. Knepley, L. Mitchell, and F. Wechsung.  Pcpatch: software for the topological construction of multigrid relaxation methods.  Transactions on Mathematical Software, 2021.
  * [6] P. E. Farrell, L. Mitchell, L. R. Scott, and F. Wechsung.  A Reynolds-robust preconditioner for the Scott-Vogelius discretization of the stationary incompressible Navier-Stokes equations.  The SMAI journal of computational mathematics, 7:75–96, 2021.
  * [7] V. Girault and P.-A. Raviart.  Finite element methods for Navier-Stokes equations. Theory and algorithms, volume 5.  Springer, 1986.
  * [8] J. Guzmán, A. Lischke, and M. Neilan.  Exact sequences on Powell–Sabin splits.  Calcolo, 57(2):1–25, 2020.
  * [9] J. Guzman, A. Lischke, and M. Neilan.  Exact sequences on Worsey-Farin splits.  arXiv preprint arXiv:2008.05431, 2020.
  * [10] J. Guzmán and M. Neilan.  Conforming and divergence-free Stokes elements in three dimensions.  IMA Journal of Numerical Analysis, 34(4):1489–1508, 2014.
  * [11] J. Guzmán and M. Neilan.  Inf-sup stable finite elements on barycentric refinements producing divergence–free approximations in arbitrary dimensions.  SIAM Journal on Numerical Analysis, 56(5):2826–2844, 2018.
  * [12] J. Guzmán and L. Scott.  The Scott-Vogelius finite elements revisited.  Mathematics of Computation, 88(316):515–529, 2019.
  * [13] J. Qin.  On the convergence of some low order mixed finite elements for incompressible fluids.  PhD thesis, Citeseer, 1994.
  * [14] J. Qin and S. Zhang.  Stability and approximability of the P1-P0 element for Stokes equations.  International journal for numerical methods in fluids, 54(5):497–515, 2007.
  * [15] L. R. Scott and M. Vogelius.  Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials.  RAIRO Modél. Math. Anal. Numér., 19(1):111–143, 1985.
  * [16] D. Silvester and A. Wathen.  Fast iterative solution of stabilised Stokes systems part ii: Using general block preconditioners.  SIAM Journal on Numerical Analysis, 31(5):1352–1367, 1994.
  * [17] U. Trottenberg, C. Ulrich Trottenberg, C. Oosterlee, A. Schuller, A. Brandt, P. Oswald, and K. Stüben.  Multigrid.  Elsevier Science, 2001.
  * [18] S. Zhang.  A new family of stable mixed finite elements for the 3D Stokes equations.  Mathematics of computation, 74(250):543–554, 2005.
  * [19] S. Zhang.  On the P1 Powell-Sabin divergence-free finite element for the Stokes equations.  Journal of Computational Mathematics, pages 456–470, 2008.
  * [20] S. Zhang.  A family of 3D continuously differentiable finite elements on tetrahedral grids.  Appl. Numer. Math., 59(1):219–233, 2009.
  * [21] S. Zhang.  Quadratic divergence-free finite elements on Powell-Sabin tetrahedral grids.  Calcolo, 48(3):211–244, 2011.