# An $H^{1}$-Conforming Solenoidal Basis for Velocity Computation on Powell-
Sabin Splits for the Stokes Problem

Jeffrey M. Connors111University of Connecticut, 341 Mansfield Road, Storrs,
CT, USA Michael Gaiewski††footnotemark:  222Corresponding author. E-mail:
<EMAIL_ADDRESS>
###### Abstract

A solenoidal basis is constructed to compute velocities using a certain finite
element method for the Stokes problem. The method is conforming, with
piecewise linear velocity and piecewise constant pressure on the Powell-Sabin
split of a triangulation. Inhomogeneous Dirichlet conditions are supported by
constructing an interpolating operator into the solenoidal velocity space. The
solenoidal basis reduces the problem size and eliminates the pressure variable
from the linear system for the velocity. A basis of the pressure space is also
constructed that can be used to compute the pressure after the velocity, if it
is desired to compute the pressure. All basis functions have local support and
lead to sparse linear systems. The basis construction is confirmed through
rigorous analysis. Velocity and pressure system matrices are both symmetric,
positive definite, which can be exploited to solve their corresponding linear
systems. Significant efficiency gains over the usual saddle-point formulation
are demonstrated computationally.

##  1 Introduction

This paper relates to finite element computations for the incompressible
Stokes problem in two dimensions. Given the real, simply connected and
polygonal domain $\Omega\subset{\mathbb{R}}^{2}$, we consider the Dirichlet
problem to solve for velocity
$\textbf{{u}}:\overline{\Omega}\to{\mathbb{R}}^{2}$ and pressure
$p:\Omega\to{\mathbb{R}}$ such that

$\displaystyle-\nu\Delta\textbf{{u}}+\nabla p$
$\displaystyle=\textbf{{f}}\quad\text{on}\ \Omega,$ (1)
$\displaystyle\nabla\cdot\textbf{{u}}$ $\displaystyle=0\quad\text{on}\
\Omega,$ (2) u $\displaystyle=\textbf{{g}}\quad\text{on}\ \partial\Omega,$ (3)
$\displaystyle\text{and}\ \int_{\Omega}p\,d\mathbf{x}$ $\displaystyle=0,$ (4)

where $\nu>0$ is a constant viscosity parameter and g is a target set of
boundary values for the velocity field. Bold font will be reserved for vectors
and spaces of vector-valued functions.

There are countless research expositions related to the Stokes problem since
it connects to many scientific models and problems in mathematics. Here, we
are primarily interested in the computation of velocity and pressure variables
using a certain finite element method where the incompressibility constraint
(2) is ultimately satisfied pointwise over the domain, yielding a true
solenoidal velocity field. In contrast, many methods only satisfy this
condition in some weak (integrated) sense, and special mixed pairings of
elements for the velocity and pressure are generally needed for strong
incompressibility. A motivation is that solenoidal velocities break a coupling
in the consistency error between velocity and pressure variables, so that
velocity computation is not polluted by errors that should only affect the
pressure accuracy. In the literature, this property is called pressure
robustness. A review of methods may be found in [SIAMRev2017].

These methods turn out to carry an additional potential benefit: to construct
a locally-supported and solenoidal basis directly for the velocity space.
Whereas a certain saddle-point problem is typically solved for the velocity
and pressure variables, a solenoidal basis can be used to express the linear
system in a block-triangular fashion that allows the velocity to be computed
without computing the pressure. If the pressure is desired, it can then be
calculated after via a separate solve, but in both cases the saddle-point
solve is replaced by smaller, symmetric positive-definite systems. Few methods
of this type exist at present, and it is not clear that a solenoidal basis can
always have a local support. There are methods of discontinuous-Galerkin (DG)
type [HL2008, MFH2008, MX2017]. Also, for Raviart-Thomas (RT), Brezzi-
Douglass-Marini (BDM) and hybridized locally-DG mixed elements that possess a
weak divergence but not a weak gradient, see [CCS2006, WWY2009, RCSR2018,
YE2021]. In the $H^{1}$-conforming case there is method with fourth-order
polynomial velocities [ParkARXIV] and a method with first-order velocities
[QS2007].

The purpose of this paper is to develop the decoupled velocity and pressure
computations for the mixed pair in [Fabien_2022]. The method enforces
solenoidal, first-order polynomial velocities and is $H^{1}$-conforming. The
element order is the same as the method in [QS2007], but the finite element
mesh is quite different. We use Powell-Sabin splits that subdivide triangles
of a fairly general mesh into six subtriangles (detailed later), whereas the
meshes of [QS2007] use rectangular meshes and subdivide each rectangle into
four triangles. This latter meshing approach may be less convenient for some
applications. Besides the geometry, the method of this paper allows for a
local macro-element assembly requiring only the six local triangles of the
Powell-Sabin split at once, grouping nodal basis functions for four nodes. The
method of [QS2007] refers to macro-element basis functions over a nine-
rectangle grid of four triangles per rectangle, hence thirty-six triangles,
and ultimately groups thirteen nodal functions together.

We summarize the paper contents as follows. Details of the Powell-Sabin mesh
and finite element method for Stokes are given in Section 2, along with some
preliminary technical lemmas. This includes a discussion about handling
Dirichlet boundary conditions. In Section 3 we construct a solenoidal basis
for the finite element velocity space with local support that can accommodate
Dirichlet conditions. The properties are proved, and we include details of the
local construction for implementation. In Section 4 we construct a local basis
of the pressure space such that each basis function is the divergence of a
known basis function in the (non-solenoidal) discrete velocity space, with
proof. This pressure basis requires a subset of vertices to be marked using a
graph-theoretic spanning tree, for which purpose there is already an
inexpensive algorithm due to Kruskal [kruskal1956]. The pressure basis is
fairly simple to implement once this is done. Computational examples are given
in Section 5 comparing the velocity and pressure computations to the classical
saddle-point system, the latter using the usual non-solenoidal basis for
velocity. A summary discussion is provided in Section 6.

##  2 A finite element method on Powell-Sabin splits

This section provides some mathematical preliminaries and then outlines the
finite element spaces studied, focusing on a specific instance from amongst
those discussed in [Fabien_2022].

###  2.1 Preliminaries

Given any countable set $S$, we denote by $|S|$ its cardinality, whereas for
sets $S\subset\Omega$, $|S|$ denotes the geometric area. Also, for
$S\subset\Omega$ we refer to the function spaces

$\displaystyle L^{2}(S)$
$\displaystyle:=\left\\{f:S\to{\mathbb{R}}\,|\,\|f\|^{2}=\int_{S}|f|^{2}\,d\textbf{{x}}<\infty\right\\}$
$\displaystyle L^{2}_{0}(S)$ $\displaystyle:=\left\\{f\in
L^{2}(S)\,|\,\int_{S}f\,d\textbf{{x}}=0\right\\}$ $\displaystyle H^{k}(S)$
$\displaystyle:=\left\\{f:S\to{\mathbb{R}}\,|\,\|f\|_{k}^{2}=\sum_{|\alpha|\leq
k}\int_{S}|D^{\alpha}f|^{2}\,d\textbf{{x}}<\infty\right\\}$ $\displaystyle
H^{k}_{0}(S)$ $\displaystyle:=\left\\{f\in
H^{k}(S)\,|\,D^{\alpha}f|_{\partial\Omega}=0,\
\forall\alpha\,\text{with}\,|\alpha|<k\right\\}$

Here, $\textbf{{x}}=(x,y)$ are coordinates in $S\subset\Omega$. Generally,
vectors and their components have the form $\textbf{{f}}=(f_{1},f_{2})$.
Vector-valued spaces use bold font as well:

$\textbf{{L}}^{2}(S)=\left[L^{2}(S)\right]^{2},\
\textbf{{H}}^{k}(S)=\left[H^{k}(S)\right]^{2},\ \text{and}\
\textbf{{H}}^{k}_{0}(S)=\left[H^{k}_{0}(S)\right]^{2}.$

In order to treat the Dirichlet boundary conditions (3), let
$\gamma_{0}:\textbf{{H}}^{1}(\Omega)\to\textbf{{H}}^{1/2}(\partial\Omega)$ be
the classical trace operator. We also use the divergence-free subspaces

$\displaystyle\textbf{{V}}(S)$
$\displaystyle:=\left\\{\textbf{{v}}\in\textbf{{H}}^{1}(S)\,|\,\nabla\cdot\textbf{{v}}=0\,\text{on}\,\Omega\right\\}$
$\displaystyle\textbf{{V}}_{0}(S)$
$\displaystyle:=\left\\{\textbf{{v}}\in\textbf{{H}}^{1}_{0}(S)\,|\,\nabla\cdot\textbf{{v}}=0\,\text{on}\,\Omega\right\\}.$

Dot-product and inner-product notations:

$\displaystyle\textbf{{f}}\cdot\textbf{{g}}$
$\displaystyle=f_{1}g_{1}+f_{2}g_{2},\
\forall\textbf{{f}}\in{\mathbb{R}}^{2},\ \textbf{{g}}\in{\mathbb{R}}^{2}$
$\displaystyle\left(f,g\right)$
$\displaystyle=\int_{\Omega}f(\textbf{{x}})g(\textbf{{x}})\,d\textbf{{x}},\
\forall f\in L^{2}(\Omega),\ g\in L^{2}(\Omega)$
$\displaystyle\left(\textbf{{f}},\textbf{{g}}\right)$
$\displaystyle=\int_{\Omega}\textbf{{f}}(\textbf{{x}})\cdot\textbf{{g}}(\textbf{{x}})\,d\textbf{{x}},\
\forall\textbf{{f}}\in\textbf{{L}}^{2}(\Omega),\
\textbf{{g}}\in\textbf{{L}}^{2}(\Omega)$
$\displaystyle\nabla\textbf{{f}}(\textbf{{x}}):\nabla\textbf{{g}}(\textbf{{x}})$
$\displaystyle=\sum_{i=1,2}\nabla f_{i}(\textbf{{x}})\cdot\nabla
g_{i}(\textbf{{x}}),\ \forall\textbf{{f}}\in\textbf{{H}}^{1}(\Omega),\
\textbf{{g}}\in\textbf{{H}}^{1}(\Omega)$
$\displaystyle\left(\nabla\textbf{{f}},\nabla\textbf{{g}}\right)$
$\displaystyle=\int_{\Omega}\nabla\textbf{{f}}(\textbf{{x}}):\nabla\textbf{{g}}(\textbf{{x}})\,d\textbf{{x}},\
\forall\textbf{{f}}\in\textbf{{H}}^{1}(\Omega),\
\textbf{{g}}\in\textbf{{H}}^{1}(\Omega).$

###  2.2 Finite element method

Let $\tau_{h}^{M}=\\{T_{k}^{M}\\}$ be a family of conforming triangulations of
$\Omega$, denoting the triangular mesh elements by $T_{k}$ for $1\leq
k\leq|\tau_{h}^{M}|$. Here, $h=\max_{k}\\{h_{k}\\}$ where
$h_{k}=\text{diam}(T_{k}^{M})$ is the maximum diameter of a triangle in the
mesh. The superscript $M$ is used to distinguish this macro-mesh from its
Powell-Sabin split, the latter denoted here by $\tau_{h}$. The mesh $\tau_{h}$
is formed via a subdivision of each triangle $T_{k}^{M}$ into six smaller
triangles.

The Powell-Sabin mesh $\tau_{h}$ is constructed via the following procedure.
Find the incenter of each $T_{k}^{M}$. This is connected via line segments to
the vertices of $T_{k}^{M}$. Given any two adjacent triangles $T_{k}^{M}$ and
$T_{j}^{M}$, their incenters are then connected via a line segment, which
necessarily intersects their shared edge. If a triangle $T_{k}^{M}$ has an
exterior edge (on $\partial\Omega$), the midpoint of the exterior edge is
connected via a line segment to the incenter. An important feature of the
Powell-Sabin split is the creation of singular points in the mesh. A singular
point is any vertex such that all edges meeting at that vertex fall on
precisely two straight lines. By construction, the mesh $\tau_{h}$ has
precisely one singular vertex along each edge of $\tau_{h}^{M}$, and no
others.

Given a vertex z, it is convenient for us to use a local indexing for the set
of triangles that share this vertex. Such a set is denoted by

$\tau^{M}_{\textbf{{z}}}=\left\\{T^{M}_{1},\ldots,T^{M}_{|\tau^{M}_{\textbf{{z}}}|}\right\\}\subset\tau_{h}^{M}$

for triangles in the macro-mesh, or else

$\tau_{\textbf{{z}}}=\left\\{T_{1},\ldots,T_{|\tau_{\textbf{{z}}}|}\right\\}\subset\tau_{h}$

for the Powell-Sabin mesh. Throughout this paper we assume an ordering of the
local triangles so that $T^{M}_{j}$ and $T^{M}_{j+1}$ (or $T_{j}$ and
$T_{j+1}$) share an edge. Moreover, we assume that a point rotated counter-
clockwise around z that crosses the shared edge will move out of triangle $j$
and into $j+1$. Graphically, the Powell-Sabin splitting is summarized in
Figure 1, along with an example of the set $\tau_{\textbf{{z}}}$ with z an
incenter of a triangle in $\tau^{M}_{h}$.

$\partial\Omega$$T_{1}$$T_{2}$$T_{3}$$T_{4}$$T_{5}$$T_{6}$ Figure 1: The
Powell-Sabin split illustrated for a triangle of the macro-mesh
$\tau_{h}^{M}$. The incenter z (black dot) is connected to the vertices (green
dots, edges) and neighboring incenters (red edges). In case of an edge on the
boundary, the incenter is connected to the midpoint of that edge (blue edge,
dot). The split results in six triangles,
$\tau_{\textbf{{z}}}=\\{T_{1},\ldots,T_{6}\\}$, and a set of singular vertices
in correspondence with the edges of $\tau_{h}^{M}$ (red and blue dots).

Some notation of use follows here, where an edge or vertex is called exterior
if it lies on $\partial\Omega$, otherwise it is interior.

  * •

$E_{int}^{M}$ (and $E_{ext}^{M}$) – interior (and exterior) edges of
$\tau_{h}^{M}$

  * •

$V_{int}^{M}$ (and $V_{ext}^{M}$) – interior (and exterior) vertices of
$\tau_{h}^{M}$

  * •

$S_{int}$ (and $S_{ext}$) – interior (and exterior) singular vertices of
$\tau_{h}$

  * •

$E^{M}=E_{int}^{M}\cup E_{ext}^{M}$, $V^{M}=V_{int}^{M}\cup V_{ext}^{M}$,
$S=S_{int}\cup S_{ext}$,

Note the cardinality relationships

$|E_{int}^{M}|=|S_{int}|\quad\text{and}\quad|E_{ext}^{M}|=|S_{ext}|.$ (5)

Given any region $S\subset\Omega$, let ${\mathbb{P}}_{k}(S)$ be the space of
polynomials of order $k$ or less on $S$, with corresponding vector-valued
space $\boldsymbol{{{\mathbb{P}}}}_{k}(S)=[{\mathbb{P}}_{k}(S)]^{2}$. Spaces
needed to discuss the discrete velocity are

$\displaystyle\textbf{{X}}_{h}$
$\displaystyle=\left\\{\textbf{{v}}_{h}\in[\mathcal{C}(\overline{\Omega})]^{2}\
|\ \textbf{{v}}_{h}\in\boldsymbol{{{\mathbb{P}}}}_{1}(T),\ \forall
T\in\tau_{h}\right\\}\subset\textbf{{H}}^{1}(\Omega)$
$\displaystyle\textbf{{X}}_{h}^{0}$
$\displaystyle=\left\\{\textbf{{v}}_{h}\in\textbf{{X}}_{h}\ |\
\textbf{{v}}_{h}=0\ \text{on}\
\partial\Omega\right\\}\subset\textbf{{H}}^{1}_{0}(\Omega)$
$\displaystyle\textbf{{V}}_{h}$
$\displaystyle=\textbf{{X}}_{h}\cap\textbf{{V}}(\Omega)$
$\displaystyle\textbf{{V}}_{h}^{0}$
$\displaystyle=\textbf{{V}}_{h}\cap\textbf{{V}}_{0}(\Omega)$

The corresponding pressure space is identified after observing a key property
of the divergence of the discrete velocity (cf. [Fabien_2022]; see [SV1985]
for a proof).

######  Lemma 1.

Let $q\in\mathcal{C}(T)$ for all $T\in\tau_{h}$ and define a function

$\theta_{\textbf{{z}}}(q):=\left\\{\begin{array}[]{ll}q|_{T_{1}}(\textbf{{z}})-q|_{T_{2}}(\textbf{{z}})+q|_{T_{3}}(\textbf{{z}})-q|_{T_{4}}(\textbf{{z}}),&\
\text{if}\ \textbf{{z}}\in S_{int},\\\
q|_{T_{1}}(\textbf{{z}})-q|_{T_{2}}(\textbf{{z}}),&\ \text{if}\
\textbf{{z}}\in S_{ext},\end{array}\right.$

where $\tau_{\textbf{{z}}}=\\{T_{1},\ldots,T_{4}\\}$ if $\textbf{{z}}\in
S_{int}$ or $\tau_{\textbf{{z}}}=\\{T_{1},T_{2}\\}$ if $\textbf{{z}}\in
S_{ext}$. Then for any $\textbf{{v}}_{h}\in\textbf{{X}}_{h}^{0}$ and singular
vertex $\textbf{{z}}\in S$ there holds
$\theta_{\textbf{{z}}}(\nabla\cdot\textbf{{v}}_{h})=0$.

The pressure space is

$P_{h}=\left\\{q_{h}\in L^{2}_{0}(\Omega)\ |\
q_{h}|_{T}\in{\mathbb{P}}_{0}(T),\,\forall
T\in\tau_{h},\,\theta_{\textbf{{z}}}(q_{h})=0,\,\forall\textbf{{z}}\in
S\right\\}.$

Before the finite element problem can be stated formally, we must address a
critical property of the discrete, solenoidal velocity space
$\textbf{{V}}_{h}$. This affects the boundary conditions and is also
referenced later to construct a solenoidal basis in Section 3. First, we
require a technical lemma.

######  Lemma 2.

Let $T\subset{\mathbb{R}}^{2}$ be a triangle with vertices $\textbf{{z}}_{i}$,
$1\leq i\leq 3$, ordered counter-clockwise going from $\textbf{{z}}_{i}$ to
$\textbf{{z}}_{i+1}$ along $\partial T$. Each triangle edge, $e_{i}$, lies
opposite the vertex $\textbf{{z}}_{i}$, with outward unit normal vector
$\textbf{{n}}_{i}$. Note the identity

$\sum_{1\leq i\leq 3}|e_{i}|\textbf{{n}}_{i}=\boldsymbol{{0}}.$ (6)

Given any $\textbf{{v}}\in\textbf{{P}}_{1}(T)$, it holds that

$\int_{T}\nabla\cdot\textbf{{v}}\,d\textbf{{x}}=-\frac{1}{2}\sum_{1\leq i\leq
3}|e_{i}|\textbf{{v}}(\textbf{{z}}_{i})\cdot\textbf{{n}}_{i}.$ (7)

###### Proof.

Let $\mathcal{A}:{\mathbb{R}}^{2}\to{\mathbb{R}}^{2}$ be the linear operator
that rotates a vector clockwise $90$ degrees. Then

$\sum_{1\leq i\leq
3}|e_{i}|\textbf{{n}}_{i}=|e_{1}|\mathcal{A}\left(\frac{\textbf{{z}}_{3}-\textbf{{z}}_{2}}{|e_{1}|}\right)+|e_{2}|\mathcal{A}\left(\frac{\textbf{{z}}_{1}-\textbf{{z}}_{3}}{|e_{2}|}\right)+|e_{3}|\mathcal{A}\left(\frac{\textbf{{z}}_{2}-\textbf{{z}}_{1}}{|e_{3}|}\right)=\textbf{{0}}.$

Then (7) follows from (6) after applying the divergence form of Green’s
Theorem. ∎

We will use $ds$ for a measure along edges.

######  Lemma 3.

Let $e^{M}\in E^{M}$ have endpoints $\textbf{{z}}_{1}$ and $\textbf{{z}}_{2}$.
The trace of any $\textbf{{v}}_{h}\in\textbf{{V}}_{h}$ along $e^{M}$ may be
uniquely determined from the values

$\textbf{{v}}_{h}(\textbf{{z}}_{1}),\ \textbf{{v}}_{h}(\textbf{{z}}_{2}),\
\text{and}\ \int_{e^{M}}\textbf{{v}}_{h}\cdot\textbf{{n}}\,ds,$ (8)

where n is a unit normal along $e^{M}$.

###### Proof.

Refer to Figure 2 to visualize the vertices, edges and normal vectors used in
this proof. Let $\textbf{{z}}_{s}\in S$ be the singular vertex of $\tau_{h}$
on $e^{M}$. Since $\textbf{{v}}_{h}$ is piecewise linear and continuous along
$e^{M}$, it suffices to show that if the values in (8) are all zero, then also
$\textbf{{v}}_{h}(\textbf{{z}}_{s})=\textbf{{0}}$. Direct integration shows
that if $\textbf{{v}}_{h}(\textbf{{z}}_{i})=\textbf{{0}}$, $i=1,2$, then

$\int_{e^{M}}\textbf{{v}}_{h}\cdot\textbf{{n}}\,ds=0\Rightarrow\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot\textbf{{n}}=0.$
(9)

It remains only to show that the tangential component is zero; equivalently
$\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot(\textbf{{z}}_{2}-\textbf{{z}}_{1})=0$.
Let $\textbf{{z}}_{c}$ be the incenter of a triangle in $\tau_{h}^{M}$ that
has edge $e^{M}$. There are two triangles, say $T_{1}$ and $T_{2}$, in
$\tau_{\textbf{{z}}_{c}}$ with edges $e_{1}$ and $e_{2}$, respectively, that
lie on $e^{M}$, opposite $\textbf{{z}}_{c}$. They have a shared edge lying
opposite $\textbf{{z}}_{i}$ on $T_{i}$. Let the edge of $T_{i}$ opposite
$\textbf{{z}}_{s}$ be $e_{i}^{s}$ and have outward unit normal
$\textbf{{n}}_{i}^{s}$ for $i=1,2$.

Since $\textbf{{v}}_{h}$ is solenoidal, we may apply (7) on each triangle
$T_{i}$ to find

$\displaystyle|e_{1}^{s}|\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot\textbf{{n}}_{1}^{s}+|e_{1}|\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}$
$\displaystyle=0$ (10)
$\displaystyle\text{and}\quad|e_{2}^{s}|\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot\textbf{{n}}_{2}^{s}+|e_{2}|\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}$
$\displaystyle=0.$ (11)

Note that edges $e_{1}^{s}$, $e_{2}^{s}$ and $e^{M}$ form a perimeter of
triangle $T_{1}\cup T_{2}$, hence the vectors $\textbf{{n}}_{1}^{s}$ and
$\textbf{{n}}_{2}^{s}$ are equivalent to a certain basis of
${\mathbb{R}}^{2}$. Then (10)-(11) show that
$\textbf{{v}}_{h}(\textbf{{z}}_{s})=\textbf{{0}}$ if
$\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}=\textbf{{0}}$. Indeed,
now sum (10)-(11) and apply (6), (9) to show that

$(|e_{1}|+|e_{2}|)\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}=-|e_{1}^{s}|\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot\textbf{{n}}_{1}^{s}-|e_{2}^{s}|\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot\textbf{{n}}_{2}^{s}\\\
=(|e_{1}|+|e_{2}|)\textbf{{v}}_{h}(\textbf{{z}}_{s})\cdot\textbf{{n}}=\textbf{{0}}.$

∎

$\textbf{{z}}_{2}$$\textbf{{z}}_{1}$$\textbf{{z}}_{c}$$\textbf{{z}}_{s}$$e^{M}=e_{1}\cup
e_{2}$$e_{1}$$e_{2}$$T_{1}$$T_{2}$$\textbf{{n}}_{1}^{s}$$e_{1}^{s}$$\textbf{{n}}_{2}^{s}$$e_{2}^{s}$
Figure 2: Labels referenced in the proof of Lemma 3.

A consequence of Lemma 3 is that we cannot generally interpolate boundary data
$\textbf{{g}}\in[\mathcal{C}(\partial\Omega)]^{2}$ to compute a solenoidal
Stokes velocity $\textbf{{u}}_{h}\in\textbf{{V}}_{h}$ by setting
$\textbf{{u}}_{h}(\textbf{{z}})=\textbf{{g}}(\textbf{{z}})$ at all the
exterior nodes z of the Powell-Sabin mesh. The following definition provides
an appropriate interpolant, where the function $\textbf{{G}}_{h}$ will be
explicitly constructed in Section 3.

######  Definition 1 (Boundary data interpolation).

Let
$\textbf{{g}}\in\textbf{{H}}^{1/2}(\partial\Omega)\cap[\mathcal{C}(\partial\Omega)]^{2}$.
We denote by $\textbf{{g}}_{h}$ the interpolant of g, where
$\textbf{{g}}_{h}=\gamma_{0}(\textbf{{G}}_{h})$ for some
$\textbf{{G}}_{h}\in\textbf{{V}}_{h}$ and

$\displaystyle\textbf{{g}}_{h}(\textbf{{z}})$
$\displaystyle=\textbf{{g}}(\textbf{{z}}),\ \forall\textbf{{z}}\in
V_{ext}^{M}$ (12)
$\displaystyle\text{and}\quad\int_{e}\textbf{{g}}_{h}\cdot\textbf{{n}}\,ds$
$\displaystyle=\int_{e}\textbf{{g}}\cdot\textbf{{n}}\,ds,\ \forall e\in
E_{ext}^{M},$ (13)

where n is the outward unit normal vector on $\partial\Omega$.

We may now state the finite element problem. We solve for
$\textbf{{u}}_{h}\in\textbf{{X}}_{h}$ and $p_{h}\in P_{h}$ by decomposing
$\textbf{{u}}_{h}=\textbf{{w}}_{h}+\textbf{{G}}_{h}$,
$\textbf{{G}}_{h}\in\textbf{{V}}_{h}$ as in Definition 1 and
$\textbf{{w}}_{h}\in\textbf{{X}}_{h}^{0}$ such that

$\displaystyle
a(\textbf{{w}}_{h},\textbf{{v}}_{h})-(p_{h},\nabla\cdot\textbf{{v}}_{h})$
$\displaystyle=L(\textbf{{v}}_{h}),\
\forall\textbf{{v}}_{h}\in\textbf{{X}}_{h}^{0}$ (14)
$\displaystyle(\nabla\cdot\textbf{{w}}_{h},q_{h})$ $\displaystyle=0,\ \forall
q_{h}\in P_{h}$ (15) $\displaystyle a(\textbf{{w}}_{h},\textbf{{v}}_{h})$
$\displaystyle:=\nu\left(\nabla\textbf{{w}}_{h},\textbf{{v}}_{h}\right)$
$\displaystyle L(\textbf{{v}}_{h})$
$\displaystyle:=(\textbf{{f}},\textbf{{v}}_{h})-a(\textbf{{G}}_{h},\textbf{{v}}_{h})$

However, it is shown in [Fabien_2022] that the solution of (14)-(15) will in
fact have $\textbf{{w}}_{h}\in\textbf{{V}}_{h}^{0}$. Thus
$\textbf{{u}}_{h}\in\textbf{{V}}_{h}$, and we focus in this paper on directly
solving the problem on the solenoidal subspace; find
$\textbf{{w}}_{h}\in\textbf{{V}}_{h}^{0}$ satisfying

$a(\textbf{{w}}_{h},\textbf{{v}}_{h})=L(\textbf{{v}}_{h}),\
\forall\textbf{{v}}_{h}\in\textbf{{V}}_{h}^{0}.$ (16)

This is accomplished in Section 3. Note that the pressure variable need not be
computed, but if it is desired we also provide a way to compute it after the
velocity is found. One solves the problem: find $p_{h}\in P_{h}$ such that

$(p_{h},\nabla\cdot\textbf{{v}}_{h})=a(\textbf{{u}}_{h},\textbf{{v}}_{h})-(\textbf{{f}},\textbf{{v}}_{h}),\
\forall\textbf{{v}}_{h}\in\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0}.$
(17)

This method is detailed in Section 4.

##  3 A solenoidal velocity basis

Explicit constructions of basis functions are provided here. We also show how
to implement boundary conditions, and a localized implementation of the basis
using a macro-element technique.

###  3.1 Bases of $\textbf{{V}}_{h}$ and $\textbf{{V}}_{h}^{0}$

Given a vertex $\textbf{{z}}\in V^{M}$, the support of each basis function
will be $\tau_{\textbf{{z}}}^{M}$. Consider the subspaces with this support:

$\textbf{{V}}_{\textbf{{z}}}=\left\\{\textbf{{v}}_{h}\in\textbf{{V}}_{h}\ |\
\textbf{{v}}_{h}(\textbf{{x}})=\textbf{{0}}\ \forall\textbf{{x}}\in T^{M},\
\forall T^{M}\subset\tau_{h}^{M}\setminus\tau_{\textbf{{z}}}^{M}\right\\}.$
(18)

We will first construct a basis of $\textbf{{V}}_{\textbf{{z}}}$. A certain
union of these across vertices z will later form the desired basis for
$\textbf{{V}}_{h}$.

######  Lemma 4.

Given any $\textbf{{z}}\in V^{M}$, it holds that
$\text{dim}(\textbf{{V}}_{\textbf{{z}}})=3$. Moreover, a basis
$\\{\boldsymbol{{\Phi}}_{1},\boldsymbol{{\Phi}}_{2},\boldsymbol{{\Phi}}_{3}\\}$
is found by taking $(\alpha_{1},\beta_{1},\delta_{1})=(1,0,0)$,
$(\alpha_{2},\beta_{2},\delta_{2})=(0,1,0)$,
$(\alpha_{3},\beta_{3},\delta_{3})=(0,0,1)$, and setting

$\displaystyle\boldsymbol{{\Phi}}_{i}(\textbf{{z}})$
$\displaystyle=(\alpha_{i},\beta_{i})$ (19)
$\displaystyle\text{and}\quad\int_{e_{j}^{M}}\boldsymbol{{\Phi}}_{i}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=\delta_{i},$ (20)

for $1\leq i\leq 3$, where, for each $T_{j}^{M}\in\tau_{\textbf{{z}}}^{M}$,
$\\{e_{j-1}^{M},e_{j}^{M}\\}\subset E^{M}$ are its edges that have z as one
endpoint, and $\textbf{{n}}_{j}$ is the unit normal on $e_{j}^{M}$ that points
in the direction of counter-clockwise rotation around z.

###### Proof.

Note that we must have

$0=\int_{T_{j}^{M}}\nabla\cdot\boldsymbol{{\Phi}}_{i}\,d\textbf{{x}}=\pm\left(\int_{e_{j}^{M}}\boldsymbol{{\Phi}}_{i}\cdot\textbf{{n}}_{j}\,ds-\int_{e_{j-1}^{M}}\boldsymbol{{\Phi}}_{i}\cdot\textbf{{n}}_{j-1}\,ds\right)$

so that the value of $\delta_{i}$ is $j$-independent in (20). Since z is the
common vertex of all $T_{j}^{M}\in\tau_{\textbf{{z}}}^{M}$, the conditions
(19)-(20) impose the same constraints on each $T_{j}^{M}$, so it suffices to
look at the restriction to one $T_{j}^{M}$. Denote by $\textbf{{z}}_{c}$ the
incenter of $T_{j}^{M}$ and $e^{M}$ the edge that lies opposite z.

Consider the vector space

$\textbf{{W}}_{j}=\left\\{\textbf{{v}}_{h}\in\mathcal{C}(T_{j}^{M})\ |\
\textbf{{v}}_{h}\in\textbf{{P}}_{1}(T),\forall T\in\tau_{\textbf{{z}}_{c}},\
\text{and}\ \textbf{{v}}_{h}=\textbf{{0}}\ \text{on}\ e^{M}\right\\}.$

The six triangles in $\tau_{\textbf{{z}}_{c}}$ together have seven distinct
vertices (see e.g. Figure 1, with $\textbf{{v}}_{h}=\textbf{{0}}$ at the three
vertices on $e^{M}$. Any $\textbf{{v}}_{h}\in\textbf{{W}}_{j}$ is uniquely
determined by its values at the other four vertices, so
$\text{dim}(\textbf{{W}}_{j})=8$. Two triangles
$T_{i}\in\tau_{\textbf{{z}}_{c}}$ have an edge lying on $e^{M}$, say $T_{1}$
and $T_{2}$ (see Figure 2). Eight linear constraints may be introduced by
requiring

$\displaystyle\textbf{{v}}_{h}(\textbf{{z}})$ $\displaystyle=(\alpha,\beta),$
$\displaystyle\quad\int_{e_{j}^{M}}\textbf{{v}}_{h}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=\delta,$

for any desired values $\alpha$, $\beta$ and $\delta$, and also setting
$\nabla\cdot\textbf{{v}}_{h}=0$ on triangles $T_{i}$, $2\leq i\leq 6$. In
fact, then $\nabla\cdot\textbf{{v}}_{h}=\textbf{{0}}$ will hold also on
$T_{1}$. This is verified as follows. Note that $T_{1}$ and $T_{2}$ each have
two vertices on $e^{M}$, where $\textbf{{v}}_{h}=\textbf{{0}}$, and their
third vertex is $\textbf{{z}}_{c}$. Apply Lemma 2 on $T_{2}$ (where we impose
$\nabla\cdot\textbf{{v}}_{h}=\textbf{{0}}$) first, which reduces to

$0=\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}},$

where n is the outward unit normal for $T_{j}^{M}$ along $e^{M}$. Application
of Lemma 2 on $T_{1}$ shows that $\nabla\cdot\textbf{{v}}_{h}$ is proportional
to $\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}$, thus
$\nabla\cdot\textbf{{v}}_{h}=\textbf{{0}}$ on $T_{1}$.

The proof reduces to showing that if $\alpha=\beta=\delta=0$ then
$\textbf{{v}}_{h}=\textbf{{0}}$ at all seven vertices. We already have
$\textbf{{v}}_{h}=\textbf{{0}}$ at the three vertices on $e^{M}$, plus at z as
well under these conditions. Since $\delta=0$, we may apply Lemma 3 to
conclude that $\textbf{{v}}_{h}=\textbf{{0}}$ along the other two edges of the
macro-triangle $T_{j}^{M}$, hence at the two singular vertices on those edges,
leaving only the incenter vertex $\textbf{{z}}_{c}$. Above we showed
$\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}=0$; repeating the
arguments (apply Lemma 2) on any of the triangles $T_{i}$ with $3\leq i\leq 6$
shows that $\textbf{{v}}_{h}(\textbf{{z}}_{c})\cdot\textbf{{n}}^{*}=0$ for
$\textbf{{n}}^{*}=\textbf{{n}}_{j}$ or $\textbf{{n}}^{*}=\textbf{{n}}_{j-1}$.
Since n and $\textbf{{n}}^{*}$ cannot be parallel, it follows
$\textbf{{v}}_{h}(\textbf{{z}}_{c})=0$. ∎

The solenoidal basis is discussed in the next result. It requires the
exclusion of a single function to form a basis, but the association made with
the boundary by choosing a vertex $\textbf{{z}}_{0}\in V_{ext}^{M}$ is for
convenience in implementing boundary conditions; see Section 3.2. Also, in
Section 4 a single boundary vertex must be chosen to extract a special basis
for the pressure space, so in practice it is natural to use the same vertex
$\textbf{{z}}_{0}$.

######  Theorem 1.

Given an indexing $\textbf{{z}}_{k}\in V^{M}$, $0\leq k\leq|V^{M}|-1$ such
that $\textbf{{z}}_{0}\in V_{ext}^{M}$ (lies on $\partial\Omega$), we denote
by $\boldsymbol{{\Phi^{(k)}_{i}}}$, $1\leq i\leq 3$, the basis functions of
$\textbf{{V}}_{\textbf{{z}}}$ defined in Lemma 4 for
$\textbf{{z}}=\textbf{{z}}_{k}$. We define sets

$\displaystyle\mathcal{B}$
$\displaystyle=\left[\bigcup_{i=1}^{3}\bigcup_{k=0}^{|V^{M}|-1}\boldsymbol{{\Phi^{(k)}_{i}}}\right]\setminus\left\\{\boldsymbol{{\Phi^{(0)}_{3}}}\right\\}$
$\displaystyle\mathcal{B}^{0}$
$\displaystyle=\left\\{\boldsymbol{{\Phi^{(k)}_{i}}}\in\mathcal{B}\,|\,\textbf{{z}}_{k}\in
V_{int}^{M}\right\\}.$

It holds that $\mathcal{B}$ is a basis of $\mathbf{V}_{h}$ and
$\mathcal{B}^{0}$ is a basis of $\mathbf{V}_{h}^{0}$.

###### Proof.

We prove that $\mathcal{B}$ is a basis of $\mathbf{V}_{h}$ by showing that
$\text{dim}(\textbf{{V}}_{h})=3|V^{M}|-1$ along with the linear independence
of $\mathcal{B}$. The former claim is a counting argument. Since
$\textbf{{V}}_{h}$ is the kernel of the divergence operator on
$\textbf{{X}}_{h}$, we may use the rank-nullity theorem. It is proved in (
[GLN2020], Theorem 4) that
$\text{dim}(\nabla\cdot\textbf{{X}}_{h})=|E^{M}|+3|\tau_{h}^{M}|$. Thus

$\text{dim}(\textbf{{V}}_{h})=\text{dim}(\textbf{{X}}_{h})-|E^{M}|-3|\tau_{h}^{M}|=2|V|-|E^{M}|-3|\tau_{h}^{M}|\\\
=2(|V^{M}|+|E^{M}|+|\tau_{h}^{M}|)-|E^{M}|-3|\tau_{h}^{M}|\\\
=2|V^{M}|+|E^{M}|-|\tau_{h}^{M}|,$ (21)

and we apply the Euler identity $|E^{M}|-|\tau_{h}^{M}|=|V^{M}|-1$ to conclude
$\text{dim}(\textbf{{V}}_{h})=3|V^{M}|-1=|\mathcal{B}|$. Next, assume
$\textbf{{v}}\in\text{span}(\mathcal{B})$ satisfies

$\textbf{{v}}=\sum_{i=1}^{3}\sum_{k=0}^{|V^{M}|-1}c_{i}^{(k)}\boldsymbol{{\Phi}}_{i}^{(k)}=\textbf{{0}}=(0,0),$

where $c_{3}^{(0)}=0$ (recall this basis function is excluded from
$\mathcal{B}$). By construction, we have at the vertices in
$\textbf{{z}}_{j}\in V^{M}$

$\textbf{{v}}(\textbf{{z}}_{j})=(c_{1}^{(j)},c_{2}^{(j)})=(0,0),\ 0\leq
j\leq|V^{M}|-1,$

and we may reduce

$\textbf{{v}}=\sum_{k=0}^{|V^{M}|-1}c_{3}^{(k)}\boldsymbol{{\Phi}}_{3}^{(k)}=\textbf{{0}}=(0,0).$

Note that each coefficient $c_{3}^{(k)}$ may be associated with a point
$\textbf{{z}}_{k}\in V^{M}$; these are joined by the edges in $E^{M}$ to form
a connected graph. Given any edge $e\in E^{M}$ with unit normal n and
endpoints $\textbf{{z}}_{n}$ and $\textbf{{z}}_{m}$, we have

$\int_{e}\textbf{{v}}\cdot\textbf{{n}}\,ds=\pm\left(c_{3}^{(n)}-c_{3}^{(m)}\right)=0,$

and therefore $c_{3}^{(n)}=c_{3}^{(m)}$ holds for all values $n$, $m$. Since
$c_{3}^{(0)}=0$, we have $c_{3}^{(k)}=0$ for all $k$. This proves linear
independence of $\mathcal{B}$. Finally, the fact that $\mathcal{B}^{0}$ is a
basis of $\textbf{{V}}_{h}^{0}$ now follows by simply restricting the boundary
values to zero. ∎

###  3.2 Implementation details

The correct representation of Dirichlet conditions is now easy to show. Recall
the decomposition $\textbf{{u}}_{h}=\textbf{{G}}_{h}+\textbf{{w}}_{h}$, where
$\textbf{{w}}_{h}\in\textbf{{V}}_{h}^{0}$ solves (16) and
$\textbf{{G}}_{h}\in\textbf{{V}}_{h}$ satisfies Definition 1. The component
$\textbf{{w}}_{h}$ is computed using the basis $\mathcal{B}^{0}$. We may
expand out $\textbf{{G}}_{h}$ using $\mathcal{B}\setminus\mathcal{B}^{0}$; for
convenience in exposition, assume the exterior edge ordering $e_{k}^{M}\in
E_{ext}^{M}$, $1\leq k\leq|E_{ext}^{M}|$, such that $e_{k}^{M}$ and
$e_{k+1}^{M}$ share a vertex. This supports a vertex labelling
$\textbf{{z}}_{k}\in V^{M}$, $0\leq k\leq|V^{M}|-1$, such that each edge
$e_{k}^{M}\in E_{ext}^{M}$ has endpoints $\textbf{{z}}_{k-1}$ and
$\textbf{{z}}_{k}$.

We may expand out

$\textbf{{G}}_{h}=\sum_{i=1}^{3}\sum_{k=0}^{|V_{ext}^{M}|-1}c_{i}^{(k)}\boldsymbol{{\Phi}}_{i}^{(k)}.$

Boundary values are still interpolated at the points $\textbf{{z}}_{k}\in
V_{ext}^{M}$, corresponding to (12):

$\textbf{{G}}_{h}\left(\textbf{{z}}_{k}\right)=(c_{1}^{(k)},c_{2}^{(k)})=\textbf{{g}}\left(\textbf{{z}}_{k}\right),\
0\leq k\leq|V_{ext}^{M}|-1.$

At split points in the Powell-Sabin mesh along the boundary, the values of
$\textbf{{G}}_{h}$ are determined implicitly through the normal moments (13).
Specifically, we have the system of equations

$\int_{e_{k}^{M}}\textbf{{G}}_{h}\cdot\textbf{{n}}\,ds=\pm(c_{3}^{(k)}-c_{3}^{(k-1)})=\int_{e_{k}^{M}}\textbf{{g}}\cdot\textbf{{n}}\,ds,\
1\leq k\leq|V_{ext}^{M}|-1,$

where n is the outward unit normal on $\partial\Omega$ and the sign depends on
whether the edges are ordered clockwise or counter-clockwise around the
boundary. It is not necessary to include case
$k=|V_{ext}^{M}|(=|E_{ext}^{M}|)$ here because the boundary data must already
satisfy the compatability condition

$\nabla\cdot\textbf{{u}}_{h}=0\Rightarrow\int_{\partial\Omega}\textbf{{G}}_{h}\cdot\textbf{{n}}\,ds=0.$

We have $c_{3}^{(0)}=0$ by convention in Theorem 1, resulting in a square,
lower-triangular system of equations here that is trivial to solve.

A localized construction of basis functions is possible using a reference map.
The reference triangle is denoted by $\hat{T}^{M}$, which has vertices
$\hat{\textbf{{z}}}_{1}=(0,0)$, $\hat{\textbf{{z}}}_{2}=(1,0)$ and
$\hat{\textbf{{z}}}_{3}=(0,1)$. Let $T_{k}^{M}\in\tau_{\textbf{{z}}}^{M}$. We
denote by $F_{k}:\hat{T}^{M}\to T_{k}^{M}$ an affine map onto $T_{k}^{M}$ with
Jacobian $J_{k}$. The reference triangle is then divided into six
subtriangles, say $\hat{T}_{i}$ for $1\leq i\leq 6$. This is done by adding
four points that are the preimages under $F_{k}$ of the singular nodes and
incenter on $T_{k}^{M}$; see Figure 3.

$(0,0)$$(1,0)$$(0,1)$$\hat{T}^{M}$$T_{k}^{M}$ Figure 3: A localized
construction of basis functions can be made explicit using a reference mapping
$F_{k}:\hat{T}^{M}\to T_{k}^{M}$. The mapping $F_{k}$ is defined as an affine
function that maps the vertices of $\hat{T}^{M}$ to those of $T_{k}^{M}$
(black). The additional points used to subdivide $\hat{T}^{M}$ are the
preimages under $F_{k}$ of the singular vertices and incenter on $T_{k}^{M}$
(red).

Functions $\textbf{{v}}_{h}\in\textbf{{V}}_{h}$ will be represented locally on
$T_{k}^{M}$ by leveraging the Piola transformation

$\textbf{{v}}_{h}=\frac{1}{|\text{det}(J_{k})|}J_{k}\hat{\textbf{{v}}}_{h}.$
(22)

As shown in the proof of Lemma 4, $\mathbf{v}_{h}$ can be determined uniquely
on $T_{k}^{M}$ via three degrees of freedom using the basis (19)-(20). Given
that $F_{k}$ maps reference vertex $\hat{\textbf{{z}}}_{j}$ to z, we can
associate three basis functions on the reference element, say
$\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}$, $1\leq i\leq 3$, such that the
divergence (in reference coordinates) is
$\hat{\nabla}\cdot\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}=0$, with

$\displaystyle\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}(\hat{\textbf{{z}}}_{j})$
$\displaystyle=(\hat{\alpha}_{i},\hat{\beta}_{i})$ (23)
$\displaystyle\text{and}\quad\int_{\hat{e}_{j}^{M}}\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}\cdot\hat{\textbf{{n}}}_{j}\,d\hat{s}$
$\displaystyle=\hat{\delta}_{i}.$ (24)

Here, $\hat{e}_{j}^{M}$ is the edge of $\hat{T}^{M}$ between
$\hat{\textbf{{z}}}_{j-1}$ and $\hat{\textbf{{z}}}_{j}$ (modulo $3$), and
$\hat{\textbf{{n}}}_{j}$ is the outward-pointing unit normal vector along that
edge (see Figure 4). We choose
$(\hat{\alpha}_{1},\hat{\beta}_{1},\hat{\delta}_{1})=(1,0,0)$,
$(\hat{\alpha}_{2},\hat{\beta}_{2},\hat{\delta}_{2})=(0,1,0)$, and
$(\hat{\alpha}_{3},\hat{\beta}_{3},\hat{\delta}_{3})=(0,0,1)$.

Since there are three such functions for each reference vertex, we have a
total of nine local basis functions on $\hat{T}^{M}$. These correspond to the
three vertices of $T_{k}^{M}$. However, care is needed to map between the
basis functions when changing variables. Given the values

$(a_{i},b_{i})=|\text{det}(J_{k})|(J_{k})^{-1}\boldsymbol{{\Phi}}_{i}(\textbf{{z}}),\
i=1,\,2,$

the correct mappings are

$\displaystyle\boldsymbol{{\Phi}}_{1}|_{T_{k}^{M}}$
$\displaystyle=\frac{a_{1}}{|\text{det}(J_{k})|}J_{k}\hat{\boldsymbol{{\Phi}}}_{1}^{(j)}+\frac{b_{1}}{|\text{det}(J_{k})|}J_{k}\hat{\boldsymbol{{\Phi}}}_{2}^{(j)},$
$\displaystyle\boldsymbol{{\Phi}}_{2}|_{T_{k}^{M}}$
$\displaystyle=\frac{a_{2}}{|\text{det}(J_{k})|}J_{k}\hat{\boldsymbol{{\Phi}}}_{1}^{(j)}+\frac{b_{2}}{|\text{det}(J_{k})|}J_{k}\hat{\boldsymbol{{\Phi}}}_{2}^{(j)},$
$\displaystyle\text{and}\ \boldsymbol{{\Phi}}_{3}|_{T_{k}^{M}}$
$\displaystyle=\frac{1}{|\text{det}(J_{k})|}J_{k}\hat{\boldsymbol{{\Phi}}}_{3}^{(j)}.$

The Piola transformation is well-known to preserve the divergence and the
integral of the normal component along edges, so it is easy to verify the
relationships above.

Let $\hat{\textbf{{z}}}_{i}^{s}$, $1\leq i\leq 3$, satisfy the requirements
that $\hat{\textbf{{z}}}_{i}^{s}$ lies on $\hat{e}_{i}^{M}$ and
$\hat{\textbf{{z}}}_{i}^{s}=(F_{k})^{-1}(\textbf{{z}}_{i}^{s})$ for some
singular vertex $\textbf{{z}}_{i}^{s}\in S$ that lies on an edge of
$T_{k}^{M}$. Also, let
$\hat{\textbf{{z}}}^{c}=(F_{k})^{-1}(\textbf{{z}}^{c})$, where
$\textbf{{z}}^{c}$ is the incenter for $T_{k}^{M}$. Labels for points on
$\hat{T}^{M}$ are shown in Figure 4. In order to describe each
$\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}$ explicitly, it suffices to list its
values at the seven points $\hat{\textbf{{z}}}_{i}$ and
$\hat{\textbf{{z}}}_{i}^{s}$, $1\leq i\leq 3$, and $\hat{\textbf{{z}}}^{c}$.
Given that
$\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}=(\hat{u}_{i}^{(j)},\hat{v}_{i}^{(j)})$,
$\hat{\textbf{{z}}}^{c}=(\hat{x}^{c},\hat{y}^{c})$ and
$\hat{\textbf{{z}}}_{i}^{s}=(\hat{x}_{i}^{s},\hat{y}_{i}^{s})$, these values
are provided in Table 1 ($j=1$), Table 2 ($j=2$), and Table 3 ($j=3$).

$\hat{\textbf{{z}}}_{1}$$\hat{\textbf{{z}}}_{2}$$\hat{\textbf{{z}}}_{3}$$\hat{\textbf{{z}}}_{3}^{s}$$\hat{\textbf{{z}}}_{1}^{s}$$\hat{e}_{1}^{M}$$\hat{e}_{2}^{M}$$\hat{e}_{3}^{M}$$\hat{\textbf{{z}}}^{c}$$\hat{\textbf{{z}}}_{2}^{s}$ Figure 4: Labels for some edges (red) and points on $\hat{T}^{M}$. | $\hat{\textbf{{z}}}^{c}$ | $\hat{\textbf{{z}}}_{1}^{s}$ | $\hat{\textbf{{z}}}_{2}^{s}$ | $\hat{\textbf{{z}}}_{3}^{s}$ | $\hat{\textbf{{z}}}_{1}$ | $\hat{\textbf{{z}}}_{2}$ | $\hat{\textbf{{z}}}_{3}$
---|---|---|---|---|---|---|---
$\hat{u}_{1}^{(1)}$ | $-\hat{y}^{c}$ | $-\hat{y}_{1}^{s}$ | $1-\hat{x}_{2}^{s}$ | 0 | 1 | 0 | 0
$\hat{v}_{1}^{(1)}$ | $\hat{y}^{c}$ | $\frac{\hat{y}_{1}^{s}-\hat{y}^{c}}{\hat{x}^{c}}$ | 0 | 0 | 0 | 0 | 0
$\hat{u}_{2}^{(1)}$ | $\hat{x}^{c}$ | 0 | $\frac{\hat{x}_{2}^{s}-\hat{x}^{c}}{\hat{y}^{c}}$ | 0 | 0 | 0 | 0
$\hat{v}_{2}^{(1)}$ | $-\hat{x}^{c}$ | $1-\hat{y}_{1}^{s}$ | $-\hat{x}_{2}^{s}$ | 0 | 1 | 0 | 0
$\hat{u}_{3}^{(1)}$ | -2 | -2 | 2$\frac{\hat{x}^{c}-\hat{x}_{2}^{s}}{\hat{y}^{c}}$ | 0 | 0 | 0 | 0
$\hat{v}_{3}^{(1)}$ | 2 | 2$\frac{\hat{y}_{1}^{s}-\hat{y}^{c}}{\hat{x}^{c}}$ | 2 | 0 | 0 | 0 | 0
Table 1:  Tabulated values of $(\hat{u}_{i}^{(1)},\hat{v}_{i}^{(1)})=\hat{\boldsymbol{{\Phi}}}_{i}^{(1)}$, $1\leq i\leq 3$, at the points shown in Figure 4. | $\hat{\textbf{{z}}}^{c}$ | $\hat{\textbf{{z}}}_{1}^{s}$ | $\hat{\textbf{{z}}}_{2}^{s}$ | $\hat{\textbf{{z}}}_{3}^{s}$ | $\hat{\textbf{{z}}}_{1}$ | $\hat{\textbf{{z}}}_{2}$ | $\hat{\textbf{{z}}}_{3}$
---|---|---|---|---|---|---|---
$\hat{u}_{1}^{(1)}$ | 0 | 0 | $\hat{x}_{2}^{s}$ | $\hat{x}_{3}^{s}+\frac{\hat{x}^{c}-\hat{x}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 1 | 0
$\hat{v}_{1}^{(1)}$ | $-\hat{y}^{c}$ | 0 | 0 | $\frac{\hat{y}^{c}-\hat{y}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
$\hat{u}_{2}^{(1)}$ | 0 | 0 | $\frac{\hat{x}_{2}^{s}-\hat{x}^{c}}{\hat{y}^{c}}$ | $\frac{\hat{x}^{c}-\hat{x}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
$\hat{v}_{2}^{(1)}$ | $\hat{x}^{c}-1$ | 0 | $\hat{x}_{2}^{s}-1$ | $\hat{x}_{3}^{s}+\frac{\hat{y}^{c}-\hat{y}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 1 | 0
$\hat{u}_{3}^{(1)}$ | 0 | 0 | 2$\frac{\hat{x}_{2}^{s}-\hat{x}^{c}}{\hat{y}^{c}}$ | 2$\frac{\hat{x}^{c}-\hat{x}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
$\hat{v}_{3}^{(1)}$ | -2 | 0 | -2 | 2$\frac{\hat{y}^{c}-\hat{y}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
Table 2:  Tabulated values of $(\hat{u}_{i}^{(2)},\hat{v}_{i}^{(2)})=\hat{\boldsymbol{{\Phi}}}_{i}^{(2)}$, $1\leq i\leq 3$, at the points shown in Figure 4. | $\hat{\textbf{{z}}}^{c}$ | $\hat{\textbf{{z}}}_{1}^{s}$ | $\hat{\textbf{{z}}}_{2}^{s}$ | $\hat{\textbf{{z}}}_{3}^{s}$ | $\hat{\textbf{{z}}}_{1}$ | $\hat{\textbf{{z}}}_{2}$ | $\hat{\textbf{{z}}}_{3}$
---|---|---|---|---|---|---|---
$\hat{u}_{1}^{(1)}$ | $\hat{y}^{c}-1$ | $\hat{y}_{1}^{s}-1$ | 0 | $\frac{\hat{y}_{3}^{s}-\hat{y}^{c}}{1-\hat{x}^{c}-\hat{y}^{c}}-\hat{x}_{3}^{s}$ | 0 | 0 | 1
$\hat{v}_{1}^{(1)}$ | 0 | $\frac{\hat{y}_{1}^{s}-\hat{y}^{c}}{\hat{x}^{c}}$ | 0 | $\frac{\hat{y}^{c}-\hat{y}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
$\hat{u}_{2}^{(1)}$ | $-\hat{x}^{c}$ | 0 | 0 | $\frac{\hat{x}^{c}-\hat{x}_{3}^{s}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
$\hat{v}_{2}^{(1)}$ | 0 | $\hat{y}_{1}^{s}$ | 0 | $\frac{\hat{x}_{3}^{s}-\hat{x}^{c}}{1-\hat{x}^{c}-\hat{y}^{c}}-\hat{x}_{3}^{s}$ | 0 | 0 | 1
$\hat{u}_{3}^{(1)}$ | 2 | 2 | 0 | 2$\frac{\hat{x}_{3}^{s}-\hat{x}^{c}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
$\hat{v}_{3}^{(1)}$ | 0 | 2$\frac{\hat{y}^{c}-\hat{y}_{1}^{s}}{\hat{x}^{c}}$ | 0 | 2$\frac{\hat{y}_{3}^{s}-\hat{y}^{c}}{1-\hat{x}^{c}-\hat{y}^{c}}$ | 0 | 0 | 0
Table 3:  Tabulated values of
$(\hat{u}_{i}^{(3)},\hat{v}_{i}^{(3)})=\hat{\boldsymbol{{\Phi}}}_{i}^{(3)}$,
$1\leq i\leq 3$, at the points shown in Figure 4.

While the construction has been localized to $T_{k}^{M}$, it is important to
note that the basis functions $\hat{\boldsymbol{{\Phi}}}_{i}^{(j)}$ depend on
the locations of the singular vertices and incenter of $T_{k}^{M}$. In
implementation, assembly may be performed by looping over triangles in
$\tau_{h}^{M}$ with a subloop over the six subtriangles of the Powell-Sabin
split, computing nodal values locally for all basis functions. In turn, once
these nodal values are found, each $\boldsymbol{{\Phi}}_{i}$ has an obvious
representation in the standard nodal basis used in continuous, piecewise
linear finite element methods, so integration on each triangle can leverage
classical techniques. Also, since the values of $\boldsymbol{{\Phi}}_{i}$ at z
are defined already by Lemma 4, it is only necessary to compute basis values
at incenters and singular vertices. Values at singular vertices can be stored
to avoid recomputation when looping over triangles in $\tau_{h}^{M}$. The
ultimate cost to assemble the linear system and solve is investigated in
Section 5 via comparison with the classical saddle-point approach, which does
not use the solenoidal basis.

##  4 Computation of pressure from velocity

It is possible to compute $p_{h}$ after $\textbf{{u}}_{h}$ has been computed,
if desired, by solving the problem posed in (17). Bases for
$\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0}$ and $P_{h}$ will be
constructed for this purpose in this section. The next result is used later to
verify the correct number of basis functions, and also shows that the linear
system for the pressure solve is square.

######  Lemma 5.

It holds that
$\dim(\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0})=\dim(P_{h})=2|\tau_{h}^{M}|+2|E_{int}^{M}|-|V_{i}^{M}|$.

###### Proof.

The pressure space is piecewise constant on $\tau_{h}$, but with the
constraints imposed by Lemma 1, which remove a degree of freedom for each
singular point. As proved in [Fabien_2022], an equivalent way to count is
three basis functions per interior singular point and one per exterior
singular point. Along with the constraint $P_{h}\subset L^{2}_{0}(\Omega)$, it
follows from the correspondence between singular points and edges of
$\tau_{h}^{M}$ that $\dim(P_{h})=3|E_{int}^{M}|+|E_{ext}^{M}|-1$. We count
three edges per triangle, so including multiplicity for interior edges we find
$3|\tau_{h}^{M}|=2|E_{int}^{M}|+|E_{ext}^{M}|$, thus
$\dim(P_{h})=3|\tau_{h}^{M}|+|E_{int}^{M}|-1$. Insert the Euler identity in
the form

$1=|\tau_{h}^{M}|+|V^{M}|-|E^{M}|=|\tau_{h}^{M}|+|V_{i}^{M}|-|E_{int}^{M}|$

and we find $\dim(P_{h})=2|\tau_{h}^{M}|+2|E_{int}^{M}|-|V_{i}^{M}|$.

Now refer to Theorem 1:
$|\mathcal{B}^{0}|=\dim(\textbf{{V}}_{h}^{0})=3|V_{i}^{M}|$. All interior
vertices on the Powell-Sabin split contribute toward the dimension of
$\textbf{{X}}_{h}^{0}$; we may group the vertices as $V_{i}^{M}$ plus interior
singular vertices and incenters of triangles in $\tau_{h}^{M}$, so that

$\dim(\textbf{{X}}_{h}^{0})=2(|V_{i}^{M}|+|E_{int}^{M}|+|\tau_{h}^{M}|)\\\
\Rightarrow\dim(\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0})=2|\tau_{h}^{M}|+2|E_{int}^{M}|-|V_{i}^{M}|=\dim(P_{h}).$

∎

The construction of basis functions builds upon the scalar finite element
space

$S_{h}=\left\\{\lambda\in\mathcal{C}(\overline{\Omega})\,|\,\lambda|_{T}\in\mathbb{P}_{1}(T),\
\forall T\in\tau_{h}\right\\}.$

Let $\\{\lambda_{\textbf{{z}}}\\}$ for all vertices $\textbf{{z}}\in V$ be the
usual Lagrangian (nodal) basis for $S_{h}$, such that for a given indexing
$\\{\textbf{{z}}_{i}\\}$ of $V$ we have
$\lambda_{\textbf{{z}}_{j}}(\textbf{{z}}_{k})=\delta_{jk}$. Next, we denote
the singular points on the edges $e_{j}^{M}\in E_{int}^{M}$ as
$\textbf{{z}}_{s}^{(j)}$ for $1\leq j\leq|E_{int}^{M}|$. On each edge
$e_{j}^{M}$, associate a unit normal vector $\textbf{{n}}_{j}$ and a unit
tangent vector $\textbf{{t}}_{j}$. The basis construction will make use of
some functions

$\boldsymbol{{\psi}}_{s}^{(j)}=\lambda_{\textbf{{z}}_{s}^{(j)}}\,\textbf{{n}}_{j}\quad\text{and}\quad\boldsymbol{{\phi}}_{s}^{(j)}=\lambda_{\textbf{{z}}_{s}^{(j)}}\,\textbf{{t}}_{j}.$

The support of these functions is the set of four triangles $T\in\tau_{h}$
with the common interior singular vertex $\textbf{{z}}_{s}^{(j)}$. If we
denote the incenters of each triangle $T_{k}^{M}\in\tau_{h}^{M}$ by
$\textbf{{z}}_{c}^{(k)}$ for $1\leq k\leq|\tau_{h}^{M}|$, then the following
are also needed:

$\displaystyle\boldsymbol{{\psi}}_{c}^{(k)}=\lambda_{\textbf{{z}}_{c}^{(k)}}\,(1,0)\quad\text{and}\quad\boldsymbol{{\phi}}_{c}^{(k)}=\lambda_{\textbf{{z}}_{c}^{(k)}}\,(0,1).$

The support of these functions is the set of six triangles $T\in\tau_{h}$ with
the common vertex $\textbf{{z}}_{c}^{(k)}$.

Some functions above must be excluded to extract a basis. This must be done
carefully. First, identify one vertex $\textbf{{z}}_{0}\in V_{ext}^{M}$ (that
lies on $\partial\Omega$; see also Theorem 1), and group it together with the
interior vertices $\textbf{{v}}_{j}\in V_{i}^{M}$, $1\leq j\leq|V_{i}^{M}|$.
There exists some subset of edges in $E_{int}^{M}$, say $e_{k}^{M}$ that we
may index (without loss of generality) as $1\leq k\leq N$, such that these
edges form a spanning tree for the graph created by the vertex set
$\\{\textbf{{z}}_{0},\ldots,\textbf{{z}}_{|V_{i}^{M}|}\\}$ and all of the
edges that have both endpoints
$\boldsymbol{{\zeta}}_{1},\boldsymbol{{\zeta}}_{2}\in\\{\textbf{{z}}_{0},\ldots,\textbf{{z}}_{|V_{i}^{M}|}\\}$.
A spanning tree of a connected graph is itself a connected subgraph that
touches every vertex in the graph without creating any cycles. These
particular vertices and edges form a finite connected graph, so a spanning
tree exists. It is well-known that $N=|V_{i}^{M}|$. We compute a spanning tree
using Kruskal’s algorithm [kruskal1956] for implementation.

######  Theorem 2.

Given the indexing above, the following set $\mathcal{S}$ is a basis of
$\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0}$:

$\mathcal{S}=\left\\{\boldsymbol{{\psi}}_{s}^{(j)}\right\\}_{j=|V_{i}^{M}|+1}^{|E_{int}^{M}|}\bigcup\left\\{\boldsymbol{{\phi}}_{s}^{(j)}\right\\}_{j=1}^{|E_{int}^{M}|}\bigcup\left\\{\boldsymbol{{\psi}}_{c}^{(k)}\right\\}_{k=1}^{|\tau_{h}^{M}|}\bigcup\left\\{\boldsymbol{{\phi}}_{c}^{(k)}\right\\}_{k=1}^{|\tau_{h}^{M}|}.$

Furthermore, $\nabla\cdot\mathcal{S}$ is a basis of $P_{h}$.

###### Proof.

Indeed, by Lemma 5 we have the correct number of basis functions. Note that
the functions in $\mathcal{S}$ are equivalent to a subset of the standard
nodal basis functions typically used for the globally continuous, piecewise
linear velocity space on $\tau_{h}$ except that at singular points the vectors
are expressed in the local basis $\\{\textbf{{n}}_{j},\textbf{{t}}_{j}\\}$ of
${\mathbb{R}}^{2}$. Since the properties of these basis functions are well-
known, it is immediately clear that
$\mathcal{S}\subset\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0}$ by
construction and that $\mathcal{S}$ is linearly independent. Therefore, it is
a basis for $\textbf{{X}}_{h}^{0}\setminus\textbf{{V}}_{h}^{0}$.

Note that each function in $\nabla\cdot\mathcal{S}$ is also in
$L^{2}_{0}(\Omega)$ (apply the divergence form of Green’s theorem), thus each
function is in $P_{h}$ by Lemma 1. It remains to show that if

$\sum_{j=|V_{i}^{M}|+1}^{|E_{int}^{M}|}a_{j}\nabla\cdot\boldsymbol{{\psi}}_{e}^{(j)}+\sum_{j=1}^{|E_{int}^{M}|}b_{j}\nabla\cdot\boldsymbol{{\phi}}_{e}^{(j)}+\sum_{k=1}^{|\tau_{h}^{M}|}c_{k}\nabla\cdot\boldsymbol{{\psi}}_{c}^{(k)}+\sum_{k=1}^{|\tau_{h}^{M}|}d_{k}\nabla\cdot\boldsymbol{{\phi}}_{c}^{(k)}=0,$

then the coefficients shown are all zero. Equivalently,
$0=\nabla\cdot\boldsymbol{{\phi}}$ where

$\boldsymbol{{\phi}}=\sum_{j=|V_{i}^{M}|+1}^{|E_{int}^{M}|}a_{j}\boldsymbol{{\psi}}_{e}^{(j)}+\sum_{j=1}^{|E_{int}^{M}|}b_{j}\boldsymbol{{\phi}}_{e}^{(j)}+\sum_{k=1}^{|\tau_{h}^{M}|}c_{k}\boldsymbol{{\psi}}_{c}^{(k)}+\sum_{k=1}^{|\tau_{h}^{M}|}d_{k}\boldsymbol{{\phi}}_{c}^{(k)}.$

Thus $\boldsymbol{{\phi}}\in\textbf{{V}}_{h}^{0}$. Since $\boldsymbol{{\phi}}$
is already expressed as a linear combination of basis functions, it suffices
to show that $\boldsymbol{{\phi}}=\textbf{{0}}$. The basis functions are non-
zero at (specific) singular points and incenters but vanish at neighboring
vertices, so $\boldsymbol{{\phi}}(\textbf{{z}})=\textbf{{0}}$ for any
$\textbf{{z}}\in V^{M}$. Upon changing to the basis of $\textbf{{V}}_{h}^{0}$
in Theorem 1, this latter property implies

$\boldsymbol{{\phi}}=\sum_{k=1}^{|V_{i}^{M}|}\beta_{k}\boldsymbol{{\Phi}}_{3}^{(k)}.$

Given any edge $e_{j}^{M}$ with $1\leq j\leq|V_{i}^{M}|$, note that

$\displaystyle\int_{e_{j}^{M}}\boldsymbol{{\psi}}_{e}^{(m)}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=0,\ |V_{i}^{M}|+1\leq m\leq|E_{int}^{M}|,$
$\displaystyle\int_{e_{j}^{M}}\boldsymbol{{\phi}}_{e}^{(m)}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=0,\ 1\leq m\leq|E_{int}^{M}|,$
$\displaystyle\int_{e_{j}^{M}}\boldsymbol{{\psi}}_{c}^{(m)}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=0,\ 1\leq m\leq|\tau_{h}^{M}|,$
$\displaystyle\int_{e_{j}^{M}}\boldsymbol{{\phi}}_{c}^{(m)}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=0,\ 1\leq m\leq|\tau_{h}^{M}|,$
$\displaystyle\Rightarrow\int_{e_{j}^{M}}\boldsymbol{{\phi}}\cdot\textbf{{n}}_{j}\,ds$
$\displaystyle=0$

In turn, we have a square, homogeneous linear system of equations for the
coefficients $\beta_{k}$:

$\sum_{k=1}^{|V_{i}^{M}|}\beta_{k}\int_{e_{j}^{M}}\boldsymbol{{\Phi}}_{3}^{(k)}\cdot\textbf{{n}}_{j}\,ds=\sum_{k=1}^{|V_{i}^{M}|}\beta_{k}\left(\pm\delta_{kj}\right)=0,\
1\leq j\leq|V_{i}^{M}|.$ (25)

Recall that the edges $e_{j}^{M}$ form a spanning tree for the vertices
$\\{\textbf{{z}}_{0},\ldots,\textbf{{z}}_{|V_{i}^{M}|}\\}$. Each basis
function $\boldsymbol{{\Phi}}_{3}^{(k)}$ is associated with vertex
$\textbf{{z}}_{k}$, and they can be put into one-to-one correspondence with
the edges $e_{j}^{M}$, $1\leq j\leq|V_{i}^{M}|$, say $\textbf{{z}}_{k}$ is an
endpoint of edge $e_{j_{k}}^{M}$. Since the labeling of interior vertices is
arbitrary, for convenience we may assume that edge $e_{j_{k}}^{M}$ has
endpoints $\textbf{{z}}_{k-1}$ and $\textbf{{z}}_{k}$. As $\textbf{{z}}_{0}$
lies on $\partial\Omega$, only basis function $\boldsymbol{{\Phi}}_{3}^{(1)}$
is supported on edge $e_{j_{1}}^{M}$, so $\beta_{1}=0$ by (25). Inductively,
given that $\beta_{k-1}=0$ we find that only $\boldsymbol{{\Phi}}_{3}^{(k)}$
is supported on edge $e_{j_{k}}^{M}$, so $\beta_{k}=0$ also for $2\leq
k\leq|V_{i}^{M}|$. Thus, $\boldsymbol{{\phi}}=\textbf{{0}}$, completing the
proof. ∎

##  5 Computations

On the domain $\Omega=(0,1)\times(0,1)\subset{\mathbb{R}}^{2}$ we will use the
manufactured solution

$\textbf{{u}}=\begin{pmatrix}\sin(x)\cos(y)\\\
-\cos(x)\sin(y)\end{pmatrix},\quad p=xy-\frac{1}{4}$ (26)

with parameter choice $\nu=1$. The usual norms for the spaces $L^{2}(\Omega)$
and $[H^{1}(\Omega)]^{2}$ are denoted by $\|\cdot\|$ and $\|\cdot\|_{1}$,
respectively. Computations that use the solenoidal basis for velocity are
labelled as SOL, as opposed to SP for the saddle-point formulation (without
the solenoidal basis). All computations are done in MATLAB using the
multifrontal direct solver package invoked by the command mldivide, commonly
referred to as the “backslash” operator. The Choleski solver is called by
mldivide for all SOL calculations in order to exploit the symmetric, positive
definite matrix properties, whereas for the SP system the solver called is
MA57. The MA57 solver is optimized for sparse, symmetric indefinite systems.
The macro-meshes $\tau_{h}^{M}$ were generated using code from [gockenbach,
yonkermesh].

\begin{overpic}[width=312.9803pt]{Graphs/Error.eps} \put(21.0,-7.0){$h$}
\put(-10.0,21.0){\rotatebox{90.0}{$||\textbf{{u}}-\textbf{{u}}_{h}||_{1}$}}
\put(-2.0,-2.5){\scriptsize$10^{-2}$} \put(20.0,-2.5){\scriptsize$10^{-1}$}
\put(42.0,-2.5){\scriptsize$10^{0}$} \put(73.0,-7.0){$h$}
\put(54.0,-2.5){\scriptsize$10^{-2}$} \put(76.0,-2.5){\scriptsize$10^{-1}$}
\put(99.0,-2.5){\scriptsize$10^{0}$} \put(-6.0,6.0){\scriptsize$10^{-2}$}
\put(51.0,1.0){\scriptsize$10^{-2}$} \put(-6.0,29.5){\scriptsize$10^{-1}$}
\put(51.0,27.0){\scriptsize$10^{-1}$} \put(-5.0,53.5){\scriptsize$10^{0}$}
\put(52.0,53.5){\scriptsize$10^{0}$} \put(7.0,49.5){SOL} \put(7.0,44.0){SP}
\put(7.0,39.25){$r=1$} \put(63.25,49.5){SOL} \put(63.25,44.0){SP}
\put(63.25,39.25){$r=1$} \put(103.0,21.0){\rotatebox{90.0}{$||p-p_{h}||$}}
\end{overpic} Figure 5: Velocity and pressure errors converge with the optimal
scaling $\|\textbf{{u}}-\textbf{{u}}_{h}\|_{1}=\mathcal{O}(h^{r})$ and
$\|p-p_{h}\|=\mathcal{O}(h^{r})$ with $r=1$. The SOL and SP velocity
calculations agree to round-off precision, so the plotted data overlaps.

Errors for velocity and pressure are shown in Figure 5 in their respective
norms. The horizontal axis shows the mesh size $h$ for the Powell-Sabin splits
$\tau_{h}$. The optimal convergence rates were observed for both velocity and
pressure variables. Between SOL and SP, the computed states agreed to within
the limits of round-off error. We also observed (not shown) that the condition
number of the matrix for the SOL velocity system was consistently less than 1%
of the corresponding SP condition number.

\begin{overpic}[width=312.9803pt]{Graphs/VelocitySOLvsSP.eps}
\put(19.0,-7.0){$|V|$} \put(-10.0,11.0){\rotatebox{90.0}{Velocity SOL
$<\text{Time}>$}} \put(103.0,19.0){\rotatebox{90.0}{SP $<\text{Time}>$}}
\put(76.0,-7.0){$|V|$} \put(7.0,48.5){Assembly} \put(7.0,43.5){Solve}
\put(7.0,38.5){Total} \put(62.0,48.5){Assembly} \put(62.0,43.5){Solve}
\put(62.0,38.5){Total} \put(-0.5,-2.5){\scriptsize 0}
\put(4.0,-2.5){\scriptsize 0.5} \put(10.75,-2.5){\scriptsize 1}
\put(15.25,-2.5){\scriptsize 1.5} \put(22.25,-2.5){\scriptsize 2}
\put(27.0,-2.5){\scriptsize 2.5} \put(33.75,-2.5){\scriptsize 3}
\put(38.5,-2.5){\scriptsize 3.5} \put(39.0,-7.0){\scriptsize$\times 10^{4}$}
\put(55.0,-2.5){\scriptsize 0} \put(60.0,-2.5){\scriptsize 0.5}
\put(66.5,-2.5){\scriptsize 1} \put(71.0,-2.5){\scriptsize 1.5}
\put(78.0,-2.5){\scriptsize 2} \put(82.75,-2.5){\scriptsize 2.5}
\put(89.5,-2.5){\scriptsize 3} \put(94.25,-2.5){\scriptsize 3.5}
\put(95.0,-7.0){\scriptsize$\times 10^{4}$} \put(-4.0,0.75){\scriptsize 0}
\put(-5.0,6.0){\scriptsize 0.2} \put(-5.0,11.25){\scriptsize 0.4}
\put(-5.0,16.5){\scriptsize 0.6} \put(-5.0,21.5){\scriptsize 0.8}
\put(-5.0,26.75){\scriptsize 1.0} \put(-5.0,32.0){\scriptsize 1.2}
\put(-5.0,37.0){\scriptsize 1.4} \put(-5.0,42.25){\scriptsize 1.6}
\put(-5.0,47.5){\scriptsize 1.8} \put(-5.0,52.5){\scriptsize 2.0}
\put(52.0,0.75){\scriptsize 0} \put(51.0,6.0){\scriptsize 0.2}
\put(51.0,11.25){\scriptsize 0.4} \put(51.0,16.5){\scriptsize 0.6}
\put(51.0,21.5){\scriptsize 0.8} \put(51.0,26.75){\scriptsize 1.0}
\put(51.0,32.0){\scriptsize 1.2} \put(51.0,37.0){\scriptsize 1.4}
\put(51.0,42.25){\scriptsize 1.6} \put(51.0,47.5){\scriptsize 1.8}
\put(51.0,52.5){\scriptsize 2.0} \end{overpic} Figure 6: Time trials comparing
velocity computation with SOL versus SP.

Run times were computed to assemble and to solve the systems. In Figure 6 we
plot $<\text{Time}>$ versus $|V|$ , where $|V|$ is the number of vertices in
the Powell-Sabin mesh and $<\text{Time}>$ is the averaged timing over 20 runs.
The timings for SOL are only for computing velocity. A jump up in solver time
was observed on the two finest meshes for the SP case; output from the solver
package indicated this was due to an increase in the cost to correctly factor
the system. Although assembly times using SOL are significant, there is an
overall time savings to compute velocity alone using SOL versus SP on fine
meshes. Solver timings are much smaller for SOL even on fairly coarse meshes,
which would translate to further efficiency gains for applications that
require repeated solves without significant reassembly.

\begin{overpic}[width=312.9803pt]{Graphs/PressureSOL.eps}
\put(-10.0,12.0){\rotatebox{90.0}{Pressure SOL $<\text{Time}>$}}
\put(50.0,-7.0){$|V|$} \put(7.0,49.0){Assembly} \put(7.0,44.0){Solve}
\put(7.0,39.0){Total} \put(-3.0,0.5){\scriptsize 0}
\put(-5.0,11.25){\scriptsize 0.05} \put(-4.5,21.75){\scriptsize 0.1}
\put(-5.0,32.0){\scriptsize 0.15} \put(-4.5,43.0){\scriptsize 0.2}
\put(-5.0,53.0){\scriptsize 0.25} \put(-2.0,-2.5){\scriptsize 0.5}
\put(13.75,-2.5){\scriptsize 1} \put(27.0,-2.5){\scriptsize 1.5}
\put(42.25,-2.5){\scriptsize 2} \put(55.5,-2.5){\scriptsize 2.5}
\put(70.75,-2.5){\scriptsize 3} \put(84.25,-2.5){\scriptsize 3.5}
\put(99.0,-2.5){\scriptsize 4} \put(99.0,-7.0){\scriptsize$\times 10^{4}$}
\end{overpic} Figure 7: SOL Pressure Time
\begin{overpic}[width=312.9803pt]{Graphs/EverythingTotal.eps}
\put(-10.0,16.0){\rotatebox{90.0}{Total $<\text{Time}>$}}
\put(50.0,-7.0){$|V|$} \put(7.0,48.0){SOL} \put(7.5,42.5){SP}
\put(-4.0,-0.25){\scriptsize 0} \put(-5.0,5.25){\scriptsize 0.2}
\put(-5.0,11.0){\scriptsize 0.4} \put(-5.0,17.0){\scriptsize 0.6}
\put(-5.0,23.0){\scriptsize 0.8} \put(-5.0,28.75){\scriptsize 1.0}
\put(-5.0,34.5){\scriptsize 1.2} \put(-5.0,40.5){\scriptsize 1.4}
\put(-5.0,46.25){\scriptsize 1.6} \put(-5.0,52.0){\scriptsize 1.8}
\put(-2.0,-3.0){\scriptsize 0.5} \put(13.75,-3.0){\scriptsize 1}
\put(27.0,-3.0){\scriptsize 1.5} \put(42.25,-3.0){\scriptsize 2}
\put(55.5,-3.0){\scriptsize 2.5} \put(70.75,-3.0){\scriptsize 3}
\put(84.25,-3.0){\scriptsize 3.5} \put(99.0,-3.0){\scriptsize 4}
\put(99.0,-7.0){\scriptsize$\times 10^{4}$} \end{overpic} Figure 8: Total time
to assemble and solve for pressure and velocity.

Timings are displayed for the computation of pressure via the SOL method in
Figure 7. These times are somewhat less than the corresponding times for the
SOL velocity computations. The total time required to compute both velocity
and pressure is compared in Figure 8. In conclusion, the SOL method supports a
more efficient system solve on fairly coarse meshes, and may be used to
compute velocity alone more quickly. Due to assembly costs, the SOL and SP
methods are about even in computing the full velocity-pressure system until
the mesh is fine enough, when SOL is again faster.

##  6 Summary

A solenoidal basis of a first-order polynomial velocity space was explicitly
constructed in this paper for the conforming, divergence-free velocity-
pressure finite element pair studied in [Fabien_2022] for the Stokes problem.
These elements are defined on the Powell-Sabin split of a triangulation for a
two-dimensional domain. Inhomogeneous Dirichlet conditions were enabled by
constructing an interpolation operator for boundary data into the trace of the
solenoidal velocity subspace. A basis of the pressure space was also derived
using divergences of non-solenoidal functions from the full discrete velocity
space. Its construction used a graph-theoretic argument to discard functions
associated with a certain subset of vertices in the mesh, retaining a spanning
set. The implementation can be accomplished by generating a spanning tree
related to the mesh, for which purpose their is an efficient algorithm due to
Kruskal [kruskal1956]. All basis functions were shown to have local support.
Rigorous proofs were included that verify the construction of bases. Also,
some implementation details for the solenoidal velocity basis were discussed
that enable a localization of the computations to assemble the linear system.

The use of the bases studied enables a block-triangular structure for the
velocity-pressure system, wherein the velocity can be isolated and computed
without computing the pressure. The computed velocity can then be used after
to compute the pressure, if desired. The non-solenoidal component of the
discrete velocity space is eliminated from the velocity computation, which
further reduces the problem size.

The system matrices for the velocity and pressure solves are both symmetric,
positive-definite and support the use of the Choleski direct solver.
Computations using MATLAB verified a significant decrease in time to compute
velocity alone compared to solving the usual saddle-point system. On fine
enough meshes, even computing both pressure and velocity together was faster
using the solenoidal velocity basis.

Other benefits are expected for situations not included in the scope of this
paper. The conjugate gradient method and associated preconditioners could be
applied for larger sized problems when a direct solve is not feasible. Also,
for applications that require repeated linear solves with a reduced cost of
system reassembly compared to the initial assembly, there could be a larger
efficiency gain using the solenoidal basis than was demonstrated in this
paper.

## References

  * [1] A. Alonso Rodríguez, J. Camaño, E. De Los Santos, and F. Rapetti.  A graph approach for the construction of high order divergence-free Raviart-Thomas finite elements.  Calcolo, 55:1–28, 2018.
  * [2] J. Carrero, B. Cockburn, and D. Schötzau.  Hybridized globally divergence-free LDG methods. Part I: the Stokes problem.  Mathematics of Computation, 75:533–563, 2006.
  * [3] Maurice Fabien, Johnny Guzmán, Michael Neilan, and Ahmed Zytoon.  Low-order divergence-free approximations for the Stokes problem on Worsey–Farin and Powell–Sabin splits.  Computer Methods in Applied Mechanics and Engineering, 390:1–24, Feb 2022.
  * [4] Gockenbach, Mark.  Understanding and Implementing the Finite Element Method.  2006\.
  * [5] J. Guzmán, A. Lischke, and M. Neilan.  Exact sequences on Powell-Sabin splits.  Calcolo, 57:1–25, 2020.
  * [6] Peter Hansbo and Mats G. Larson.  Piecewise divergence-free discontinuous Galerkin methods for Stokes flow.  Communications in Numerical Methods in Engineering, 24(5):355–366, 2008.
  * [7] Volker John, Alexander Linke, Christian Merdon, Michael Neilan, and Leo G. Rebholz.  On the divergence constraint in mixed finite element methods for incompressible flows.  SIAM Review, 59(3):492–544, 2017.
  * [8] Joseph B Kruskal.  On the shortest spanning subtree of a graph and the traveling salesman problem.  Proceedings of the American Mathematical society, 7(1):48–50, 1956\.
  * [9] A. Montlaur, S. Fernandez-Mendez, and A. Huerta.  Discontinuous Galerkin methods for the Stokes equations using divergence-free approximations.  International Journal for Numerical Methods in Fluids, 57(9):1071–1092, 2008.
  * [10] Lin Mu and Xiu Ye.  A simple finite element method for the Stokes equations.  Advances in Computational Mathematics, 43:1305–1324, 2017.
  * [11] Chunjae Park.  A locally calculable $P^{3}$-pressure in a decoupled method for incompressible Stokes equations, arXiv:2107.08472, 2021.
  * [12] Jinshui Qin and Shangyou Zhang.  Stability and approximability of the P1-P0 element for Stokes equations.  International Journal for Numerical Methods in Fluids, 54(5):497–515, 2007.
  * [13] L. R. Scott and M. Vogelius.  Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials.  ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 19(1):111–143, 1985.
  * [14] Junping Wang, Yanqiu Wang, and Xiu Ye.  A robust numerical method for stokes equations based on divergence-free H(div) finite element methods.  SIAM Journal on Scientific Computing, 31(4):2784–2802, 2009.
  * [15] Xiu Ye and Shangyou Zhang.  A numerical scheme with divergence free H-div triangular finite element for the Stokes equations.  Applied Numerical Mathematics, 167:211–217, 2021.
  * [16] Yonker, Shea.  Triangulation manipulation program.  https://github.com/syonker/Triangulation-Manipulation-Program, April 2018\.