# Extraction and application of super-smooth cubic B-splines over
triangulations

Jan Grošelj  Hendrik Speleers

###### Abstract

The space of $C^{1}$ cubic Clough–Tocher splines is a classical finite element
approximation space over triangulations for solving partial differential
equations. However, for such a space there is no B-spline basis available,
which is a preferred choice in computer aided geometric design and
isogeometric analysis. A B-spline basis is a locally supported basis that
forms a convex partition of unity. In this paper, we explore several
alternative $C^{1}$ cubic spline spaces over triangulations equipped with a
B-spline basis. They are defined over a Powell–Sabin refined triangulation and
present different types of $C^{2}$ super-smoothness. The super-smooth
B-splines are obtained through an extraction process, i.e., they are expressed
in terms of less smooth basis functions. These alternative spline spaces
maintain the same optimal approximation power as Clough–Tocher splines. This
is illustrated with a selection of numerical examples in the context of least
squares approximation and finite element approximation for second and fourth
order boundary value problems.

Keywords: Triangular finite elements, $C^{1}$ cubic splines, B-spline basis,
Super-smoothness

##  1 Introduction

Splines over triangulations are a common tool for bivariate approximation.
They are inextricably linked to the finite element analysis, where they appear
as triangular elements and are used in the process of finding an approximate
solution to a partial differential equation. For many standard problems, it is
sufficient to work with continuous (or even discontinuous) splines, which are
easy to represent. Higher order problems, however, require smoother elements,
a prominent example of which is the $C^{1}$ cubic Clough–Tocher element
(Clough and Tocher, 1965).

The advent of isogeometric analysis has put a particular focus on the
representation of approximations in terms of spline basis functions. The
driving idea of this method is to use the same spline model for design and
analysis (Cottrell et al., 2009). In computer aided geometric design, B-spline
functions are a standard tool for representing geometric objects, and here we
understand these functions in a general sense, as functions that are locally
supported and form a convex partition of unity. Furthermore, B-splines
commonly have a high order of smoothness, which has proven to be beneficial in
achieving an excellent accuracy with respect to the number of degrees of
freedom and in obtaining a better spectral behavior in comparison with less
smooth finite elements; see, e.g., Hughes et al. (2008); Evans et al. (2009);
Bressan and Sande (2019); Sande et al. (2019); Manni et al. (2022).

In the context of isogeometric analysis, the potential of smooth spline spaces
over triangulations has been tested in several papers (Speleers et al., 2012;
Jaxon and Qian, 2014; Beirão da Veiga et al., 2015; Speleers and Manni, 2015;
Liu and Jeffers, 2018; Zareh and Qian, 2019), but is not fully exploited yet.
Compared with tensor-product spline spaces, using triangulations for tiling
the domain offers more flexibility in describing the geometry and performing
local refinements, however, this comes at the expense of more complicated
spline constructions. Over a general unstructured triangulation, the dimension
of $C^{1}$ spline spaces depends on the triangulation’s geometry (Lai and
Schumaker, 2007), and so even identifying the right number of basis functions
may be a tedious endeavor, possibly subjected to numerical instabilities. For
this reason, it is customary to refine triangulations with particular splits
(e.g., the Clough–Tocher split of each triangle into three smaller ones) which
provide some structure to the spline partition and allow the association of
degrees of freedom with topological properties, namely the vertices, edges,
and triangles of the triangulation.

The next step in an attempt to bring triangular elements closer to the ideas
of isogeometric analysis is the establishment of suitable basis functions.
Despite recent developments in B-spline techniques for finite elements
(Dierckx, 1997; Speleers, 2010a, b, 2013; Grošelj and Knez, 2021, 2022), the
standard $C^{1}$ cubic Clough–Tocher element still lacks a B-spline
representation, and it is not even clear whether such a representation exists.
In the absence of it, we explore in this paper alternative $C^{1}$ cubic
spline spaces over triangulations that are readily equipped with B-spline
basis functions (Grošelj and Speleers, 2017). They are defined over a
Powell–Sabin refined triangulation, which may be seen as a continuation of the
Clough–Tocher split. Moreover, we focus on $C^{2}$ super-smoothness properties
and develop an extraction process that enables a simple transition from less
to more smooth B-spline basis functions.

It has been observed by Grošelj and Speleers (2021) that the general $C^{1}$
cubic Powell–Sabin B-splines admit additional $C^{2}$ super-smoothness
properties when applied to three-directional meshes and that further
smoothness properties can be obtained by combining and reducing certain
B-splines. The aim of this paper is to provide an extraction procedure that
applies those techniques to unstructured triangulations wherever possible and
to numerically demonstrate the benefits of the procedure in approximation
methods. As the main framework, we consider B-spline representations over an
initial triangulation, which can be, for approximation purposes, uniformly
refined. In this way, we acquire a three-directional triangulation on each
initial triangle, where the symmetries can be exploited to reduce the number
of degrees of freedom and to impose local $C^{2}$ super-smoothness.

In general, obtaining maximal smoothness for splines over triangulations with
respect to the polynomial degree has proven to be very difficult. The
techniques that come the closest to achieving this goal are based on simplex
splines. This generalization of univariate B-splines is currently undergoing
an exciting development (for some recent advances, see Neamtu, 2007; Liu and
Snoeyink, 2007; Schmitt, 2021; Barucq et al., 2022; Lyche et al., 2022) and
has been considered for use in the finite element method and isogeometric
analysis (Cao et al., 2019; Wang et al., 2022). However, their disadvantage
appears to be a rather complicated underlying configuration of knot lines,
which makes the assembly and (numerical) integration a cumbersome problem when
setting up the linear systems.

In the recent literature, several spline constructions have also been proposed
over unstructured quadrilateral partitions consisting of large structured
parts. In that context it is a common practice to use standard (tensor-
product) B-spline basis functions supported over the structured parts of such
a partition and specifically tailored basis functions (with often lower
smoothness) in the vicinity of the interfaces between the structured parts;
see, e.g., Karčiauskas et al. (2016); Toshniwal et al. (2017); Wei et al.
(2018); Chan et al. (2018); Kapl et al. (2019); Weinmüller and Takacs (2022).
The latter basis functions are usually not without shortcomings, e.g., they
might require locally higher degrees, might not be nonnegative, might not be
$C^{1}$ smooth, and/or might negatively affect the approximation power of the
resulting space. In this paper, we reverse this approach by starting with
basis functions that are designed for general unstructured partitions and have
none of these shortcomings, and then we reduce the basis functions over the
structured parts of the partition to achieve locally higher smoothness. This
process is described by extraction matrices allowing a fast transition from
one approximation space to the other.

More particularly, we begin with the full space spanned by the $C^{1}$ cubic
Powell–Sabin B-splines introduced by Grošelj and Speleers (2017); these are
associated with the vertices and edges of the triangulation. Then, firstly, we
reduce the B-splines associated with the edges, which results in the spline
space introduced by Speleers (2015) and asymptotically (with respect to the
refinement of the initial triangulation) diminishes the dimension of the full
space by a factor $5/3$. Secondly, we combine the previously reduced edge
B-splines over all triangles that possess a specific geometric symmetry. The
resulting space is locally the same as the one introduced by Grošelj and
Speleers (2021) and asymptotically diminishes the dimension of the full space
by a factor $3$.

All three considered spaces reproduce cubic polynomials and have optimal
approximation power, which we illustrate with a selection of numerical
examples in the context of least squares approximation and finite element
approximation for second and fourth order boundary value problems. Through
these examples we also demonstrate that locally increased smoothness improves
accuracy of the approximations with respect to the number of degrees of
freedom, i.e., the reduced spaces prove to be computationally more efficient.

The remainder of the paper is organized as follows. In Section 2, we present
$C^{1}$ cubic Powell–Sabin B-splines and explain the extraction process. In
Section 3, we derive dimension formulas for the considered spaces and
investigate their asymptotic behavior. In Section 4, we describe the
application of the extraction process to a selection of approximation methods
and provide a few numerical examples. We conclude with some remarks in Section
5.

##  2 Basis functions

Suppose $\triangle$ is an arbitrary triangulation of a planar bounded domain
$\Omega$, and $\triangle^{\ell}$, $\ell\in\mathbb{N}$, a refinement of
$\triangle$ obtained by splitting each triangle of $\triangle$ uniformly such
that each edge of $\triangle$ is split into $\ell$ intervals as shown in
Figure 1. In this section, we will consider different sets of $C^{1}$ cubic
B-spline basis functions over the triangulation $\triangle^{\ell}$.

###  2.1 Powell–Sabin refinement

Let us represent the triangulation $\triangle^{\ell}$ by the triplet
$(V^{\ell},E^{\ell},T^{\ell})$, where
$V^{\ell}=\\{v_{1},v_{2},\ldots,v_{n_{v}^{\ell}}\\}$ is the set of vertices,
$E^{\ell}=\\{e_{1},e_{2},\ldots,e_{n_{e}^{\ell}}\\}$ the set of edges, and
$T^{\ell}=\\{t_{1},t_{2},\ldots,t_{n_{t}^{\ell}}\\}$ the set of triangles.
Note that $n_{v}^{\ell}$, $n_{e}^{\ell}$, and $n_{t}^{\ell}$ denote the number
of elements in $V^{\ell}$, $E^{\ell}$, and $T^{\ell}$, respectively. We
further divide the set of triangles $T^{\ell}$ into three subsets
$T^{\ell}_{r}$, $r=0,1,2$. A triangle $t_{k}\in T^{\ell}$ belongs to
$T^{\ell}_{r}$ if $r$ is the maximum number for which there exists an edge of
$\triangle$ that contains $r$ vertices of $t_{k}$; see Figure 1 for an
illustration.

Figure 1: A part of a triangulation $\triangle$ (depicted by thicker lines)
that is refined to the triangulation $\triangle^{\ell}$ with $\ell=8$. The
color of a triangle in $\triangle^{\ell}$ determines the set $T^{\ell}_{r}$,
$r\in\\{0,1,2\\}$, to which the triangle belongs.

To introduce basis functions over $\triangle^{\ell}$, we refine
$\triangle^{\ell}$ by applying the Powell–Sabin 6-split (Powell and Sabin,
1977). For each triangle $t_{k}\in T^{\ell}$, we choose a split point
$v_{k}^{t}$ inside $t_{k}$ such that for any neighboring triangle
$t_{k^{\prime}}\in T^{\ell}$ the line between $v_{k}^{t}$ and
$v_{k^{\prime}}^{t}$ intersects the interior of the edge $e_{j}\in E^{\ell}$
that is common to $t_{k}$ and $t_{k^{\prime}}$. The intersection point is
denoted by $v_{j}^{e}$. If $e_{j}$ is a boundary edge, the split point
$v_{j}^{e}$ is the midpoint of the edge.

In our setting, a preferred position of $v_{k}^{t}$ is the barycenter of
$t_{k}$, which is always a valid choice when $t_{k}\in T^{\ell}_{0}\cup
T^{\ell}_{1}$ as any two neighboring triangles in a uniform subtriangulation
form a strictly convex quadrilateral. In general, however, choosing the
barycenter of $t_{k}\in T^{\ell}_{2}$ could result in a degenerate or ill-
poised refinement setting. On the other hand, choosing $v_{k}^{t}$ to be the
incenter of $t_{k}$ guarantees the intersection property.

Each $t_{k}$ is split into six smaller triangles by connecting $v_{k}^{t}$ to
the vertices of $t_{k}$ and to the split points on the edges of $t_{k}$. We
denote the induced refinement of $\triangle^{\ell}$ by
$\triangle_{\mathrm{PS}}^{\ell}$. An example of such a refinement is shown in
Figure 2.

Suppose a triangle $t_{k}$ of $\triangle^{\ell}$ lies inside a uniform
subtriangulation, i.e., if $v_{i_{1}}$, $v_{i_{2}}$, $v_{i_{3}}$ are the
vertices of $t_{k}$, then its three neighboring triangles in
$\triangle^{\ell}$ are determined by the vertices

$v_{i_{1}},v_{i_{2}},v_{i_{1}}+v_{i_{2}}-v_{i_{3}},\qquad
v_{i_{2}},v_{i_{3}},v_{i_{2}}+v_{i_{3}}-v_{i_{1}},\qquad
v_{i_{3}},v_{i_{1}},v_{i_{3}}+v_{i_{1}}-v_{i_{2}}.$

Suppose also that the split points of these four triangles are barycenters. We
then say that $t_{k}$ has symmetric configuration. The set of all triangles in
$\triangle^{\ell}$ with symmetric configuration is denoted by
$T_{\mathrm{sym}}^{\ell}$; see Figure 2 for an illustration.

Figure 2: An excerpt from the triangulation $\triangle^{\ell}$ shown in Figure
1 and its Powell–Sabin refinement $\triangle_{\mathrm{PS}}^{\ell}$. Each
triangle of $\triangle^{\ell}$ (depicted by dark gray lines) is split into six
smaller triangles in $\triangle_{\mathrm{PS}}^{\ell}$ (depicted by light gray
lines). The split points are colored gray. The ones inside triangles of
$\triangle^{\ell}$ are either barycenters or incenters, which are represented
by circles or stars, respectively. The split points on the edges are marked by
squares. The triangles that belong to $T_{\mathrm{sym}}^{\ell}$ are red
shaded. The vertices of $\triangle^{\ell}$ are represented by black points.
For each $v_{i}\in V^{\ell}$, there is a triangle $q_{i}^{v}$ depicted in
blue.

###  2.2 Full space

Over the triangulation $\triangle_{\mathrm{PS}}^{\ell}$ of the domain
$\Omega$, we consider the space
$\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell})$ consisting of $C^{1}$
functions such that the restriction to any triangle of
$\triangle_{\mathrm{PS}}^{\ell}$ is a bivariate polynomial of total degree
less than or equal to $3$. The dimension of this space is
$3n_{v}^{\ell}+4n_{e}^{\ell}$; we refer the reader to Grošelj and Speleers
(2017) for more details.

The space $\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell})$ can be represented
by basis functions $B_{1},B_{2},\ldots,B_{m}$, $m=3n_{v}^{\ell}$, associated
with the vertices of $\triangle^{\ell}$ and basis functions
$B_{1}^{0},B_{2}^{0},\ldots,B_{m_{0}}^{0}$, $m_{0}=4n_{e}^{\ell}$, associated
with the edges of $\triangle^{\ell}$. We follow the approach from Grošelj and
Speleers (2017) for their construction.

To each vertex $v_{i}\in V^{\ell}$, we assign a triangle $q_{i}^{v}$ that
contains $v_{i}$ and the points $\frac{2}{3}v_{i}+\frac{1}{3}v_{j}^{e}$ and
$\frac{2}{3}v_{i}+\frac{1}{3}v_{k}^{t}$ for each edge $e_{j}\in E^{\ell}$ and
each triangle $t_{k}\in T^{\ell}$ attached to $v_{i}$; an illustration of such
triangles can be found in Figure 2. The vertices of $q_{i}^{v}$ happen to be
Greville points (see Grošelj and Speleers, 2017, for details). To describe the
other parameters in the construction, let us consider the edge $e_{j}\in
E^{\ell}$ with the endpoints $v_{i},v_{i^{\prime}}\in V^{\ell}$ and the split
point $v_{j}^{e}$. For each triangle of $\triangle_{\mathrm{PS}}$ determined
by the vertices $v_{i}$, $v_{j}^{e}$, $v_{k}^{t}$, we can choose a parameter
$\sigma_{i,j,k}>0$, and for a possible boundary edge of
$\triangle_{\mathrm{PS}}$ determined by the endpoints $v_{i}$, $v_{j}^{e}$, we
can choose a parameter $\sigma_{i,j}>0$. These yield the points
$\frac{2}{3}v_{i}+\frac{1}{3}v_{k}^{t}+\frac{2}{3}\sigma_{i,j,k}(v_{i^{\prime}}-v_{i})$
and
$\frac{2}{3}v_{i}+\frac{1}{3}v_{j}^{e}+\frac{2}{3}\sigma_{i,j}(v_{i^{\prime}}-v_{i})$,
which play again the role of Greville points. An admissible choice that is
assumed in the remainder of the paper is
$\sigma_{i,j,k}=\sigma_{i,j}=\frac{1}{2}$.

Any triangle $q_{i}^{v}$ associated with $v_{i}\in V^{\ell}$ specifies three
basis functions. Examples of such triplets of functions are shown in Figure 3
(left). These basis functions $B_{1},B_{2},\ldots,B_{m}$ possess the following
super-smoothness properties; see Grošelj and Speleers (2017, 2021).

######  Theorem 1.

The functions $B_{1},B_{2},\ldots,B_{m}$ possess the following super-
smoothness properties.

  1. 1.

For each $t_{k}\in T^{\ell}$, the functions are $C^{2}$ smooth at the split
point $v_{k}^{t}$.

  2. 2.

For each $t_{k}\in T^{\ell}$ and each edge $e_{j}$ of $t_{k}$, the functions
are $C^{2}$ smooth across the edge connecting $v_{j}^{e}$ and $v_{k}^{t}$.

  3. 3.

For each $t_{k}\in T_{\mathrm{sym}}^{\ell}$, the functions are $C^{2}$ smooth
over the interior of the triangle $t_{k}$.

Figure 3: Examples of basis functions $B_{1},B_{2},\ldots,B_{m}$ (left) and
$B_{1}^{0},B_{2}^{0},\ldots,B_{m_{0}}^{0}$ (right) over
$\triangle_{\mathrm{PS}}^{\ell}$ shown in Figure 2. The left picture shows the
triplets of basis functions associated with three different types of vertices
of $\triangle^{\ell}$: a vertex encircled by triangles from
$T_{\mathrm{sym}}^{\ell}$, a vertex on an inner edge of $\triangle$, and a
vertex on the boundary of the domain. The right picture shows the quartets of
basis functions associated with three different types of edges of
$\triangle^{\ell}$: an edge shared by triangles from
$T_{\mathrm{sym}}^{\ell}$, an edge lying on an inner edge of $\triangle$, and
an edge on the boundary of the domain.

Additionally, there are four basis functions associated with each edge of
$\triangle^{\ell}$. Examples of them are shown in Figure 3 (right). To
distinguish between these basis functions
$B_{1}^{0},B_{2}^{0},\ldots,B_{m_{0}}^{0}$, we assign them to triangles and
boundary edges of $\triangle_{\mathrm{PS}}^{\ell}$. Let $\Xi_{0}$ be a subset
of $V^{\ell}\times E^{\ell}\times(T^{\ell}\cup\\{\emptyset\\})$ consisting of
$m_{0}$ elements. For each edge $e_{j}$ of $\triangle^{\ell}$, it contains
four triplets $(v_{i},e_{j},t_{k})$, $(v_{i^{\prime}},e_{j},t_{k})$,
$(v_{i},e_{j},t_{k^{\prime}})$, $(v_{i^{\prime}},e_{j},t_{k^{\prime}})$ that
identify

  * •

four triangles of $\triangle_{\mathrm{PS}}^{\ell}$ attached to $e_{j}$ if
$e_{j}$ is not a boundary edge, or

  * •

two triangles and two boundary edges of $\triangle_{\mathrm{PS}}^{\ell}$ if
$e_{j}$ is a boundary edge ($t_{k^{\prime}}=\emptyset$).

We denote by $\xi_{0}$ the bijection from $\Xi_{0}$ to
$\\{1,2,\ldots,m_{0}\\}$.

Contrary to $B_{1},B_{2},\ldots,B_{m}$, the functions
$B_{1}^{0},B_{2}^{0},\ldots,B_{m_{0}}^{0}$ in general do not possess any
additional super-smoothness properties. However, they complement the basis
functions associated with the vertices in the following sense; see Grošelj and
Speleers (2017).

######  Theorem 2.

The functions $B_{1},B_{2},\ldots,B_{m}$ and the functions
$B_{1}^{0},B_{2}^{0},\ldots,B_{m_{0}}^{0}$ form a convex partition of unity
and a locally supported basis of
$\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell})$.

In the next sections, we consider two space reductions with additional local
super-smoothness. The reduced basis functions are obtained through an
extraction process, i.e., they are expressed in terms of less smooth basis
functions.

###  2.3 First space reduction

Let us consider a subset $\Xi_{1}$ of
$E^{\ell}\times(T^{\ell}\cup\\{\emptyset\\})$ with $m_{1}=2n_{e}^{\ell}$
elements. For each edge $e_{j}\in E^{\ell}$ it contains two pairs
$(e_{j},t_{k})$, $(e_{j},t_{k^{\prime}})$ where $t_{k}$ and $t_{k^{\prime}}$
are triangles of $\triangle^{\ell}$ attached to $e_{j}$. If $e_{j}$ is a
boundary edge, then $t_{k^{\prime}}=\emptyset$. Suppose $\xi_{1}$ is a
bijection from $\Xi_{1}$ to $\\{1,2,\ldots,m_{1}\\}$.

Let $\mathcal{B}_{0}:\Omega\rightarrow\mathbb{R}^{m_{0}}$ be a mapping that
assigns to the point $p\in\Omega$ the column vector
$[B_{1}^{0}(p),B_{2}^{0}(p),\ldots,B_{m_{0}}^{0}(p)]^{T}$. Furthermore, let
$H_{1}\in\mathbb{R}^{m_{1}\times m_{0}}$ be the matrix with entries equal to
$1$ at positions $(\xi_{1}(e_{j},t_{k}),\xi_{0}(v_{i},e_{j},t_{k}))$ and
$(\xi_{1}(e_{j},t_{k}),\xi_{0}(v_{i^{\prime}},e_{j},t_{k}))$ for each
$(e_{j},t_{k})\in\Xi_{1}$, where $v_{i},v_{i^{\prime}}\in V^{\ell}$ are the
endpoints of $e_{j}$. All other entries of $H_{1}$ are equal to $0$. The
mapping $\mathcal{B}_{1}:\Omega\rightarrow\mathbb{R}^{m_{1}}$ defined by

$\mathcal{B}_{1}(p)=H_{1}\cdot\mathcal{B}_{0}(p)$

determines a new set of functions, which we denote by
$B_{1}^{1},B_{2}^{1},\ldots,B_{m_{1}}^{1}$. Examples of them are shown in
Figure 4 (left). This procedure results in the reduced set of basis functions
developed by Speleers (2015). We refer the reader to Grošelj and Speleers
(2017) for more details on the relation between the two sets of basis
functions $\mathcal{B}_{0}$ and $\mathcal{B}_{1}$.

Figure 4: Examples of basis functions
$B_{1}^{1},B_{2}^{1},\ldots,B_{m_{1}}^{1}$ (left) and
$B_{1}^{2},B_{2}^{2},\ldots,B_{m_{2}}^{2}$ (right) over
$\triangle_{\mathrm{PS}}^{\ell}$ shown in Figure 2. The left picture shows the
pairs of basis functions associated with three different types of edges of
$\triangle^{\ell}$: an edge shared by triangles from
$T_{\mathrm{sym}}^{\ell}$, an edge lying on an inner edge of $\triangle$, and
an edge on the boundary of the domain. The right picture shows six functions:
two basis functions spanning over four triangles of $\triangle^{\ell}$
associated with two triangles from $T_{\mathrm{sym}}^{\ell}$, a basis function
spanning over two triangles of $\triangle^{\ell}$, one of which is in
$T_{\mathrm{sym}}^{\ell}$ and the other is not, and three basis functions that
remain unchanged after the second reduction (two of which are already shown on
the left). For interpretation of the indices in the legends, see Figures 5 and
6.

######  Theorem 3.

The functions $B_{1}^{1},B_{2}^{1},\ldots,B_{m_{1}}^{1}$ possess the following
super-smoothness properties.

  1. 1.

For each $t_{k}\in T^{\ell}$, the functions are $C^{2}$ smooth at the split
point $v_{k}^{t}$.

  2. 2.

For each $t_{k}\in T^{\ell}$ and each edge $e_{j}$ of $t_{k}$, the functions
are $C^{2}$ smooth across the edge connecting $v_{j}^{e}$ and $v_{k}^{t}$.

Let $\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})$ be the span of the
functions $B_{1},B_{2},\ldots,B_{m}$ and the functions
$B_{1}^{1},B_{2}^{1},\ldots,B_{m_{1}}^{1}$. The dimension of this space is
$3n_{v}^{\ell}+2n_{e}^{\ell}$. From Speleers (2015) we also know the following
properties.

######  Theorem 4.

The functions $B_{1},B_{2},\ldots,B_{m}$ and the functions
$B_{1}^{1},B_{2}^{1},\ldots,B_{m_{1}}^{1}$ form a convex partition of unity
and a locally supported basis of
$\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})$.

######  Example 5.

Figure 5 shows a few rows and columns of the matrix $H_{1}$ used for the
reduction of the three quartets of basis functions depicted in Figure 3
(right) to the three pairs of basis functions depicted in Figure 4 (left).
Thanks to the properties of the original basis functions (see Theorem 2), this
matrix clearly confirms that the reduced basis functions form a convex
partition of unity. Indeed, all its entries are nonnegative and the entries of
each column sum to one. The local supports of the reduced basis functions
follow from the sparsity of the matrix.

######  Remark 6.

A matrix that describes a set of basis functions in terms of another set of
(less smooth) basis functions is typically called an extraction matrix. This
terminology was introduced by Borden et al. (2011) and Scott et al. (2011) to
denote extraction of spline functions in terms of Bernstein polynomials. More
recent examples of extraction matrices can be found in Speleers (2019) and
Speleers and Toshniwal (2021).

Figure 5: An example of the first reduction. The basis functions shown in the
left picture of Figure 4 are obtained from the basis functions shown in the
right picture of Figure 3. An excerpt from the matrix $H_{1}$ contains the six
rows that describe this reduction.

###  2.4 Second space reduction

Let
$\Xi_{2}\subset(E^{\ell}\cup\\{\emptyset\\})\times(T^{\ell}\cup\\{\emptyset\\})$
be the union of $\\{\emptyset\\}\times T_{\mathrm{sym}}^{\ell}$ and
$\Xi_{1}\cap(E^{\ell}\times(T^{\ell}\setminus
T_{\mathrm{sym}}^{\ell}\cup\\{\emptyset\\}))$, and let $m_{2}$ denote the
number of elements of this set. Suppose $\xi_{2}$ is a bijection from
$\Xi_{2}$ to $\\{1,2,\ldots,m_{2}\\}$.

Let $H_{2}\in\mathbb{R}^{m_{2}\times m_{1}}$ be the matrix with entries
defined as follows. For each $t_{k}\in T_{\mathrm{sym}}^{\ell}$ and each edge
$e_{j}$ of $t_{k}$, the entry at position
$(\xi_{2}(\emptyset,t_{k}),\xi_{1}(e_{j},t_{k}))$ is equal to $\frac{2}{3}$.
For the second triangle $t_{k^{\prime}}\in T^{\ell}$ attached to $e_{j}$, the
entry at position $(\xi_{2}(\emptyset,t_{k}),\xi_{1}(e_{j},t_{k^{\prime}}))$
is equal to $\frac{1}{3}$. Moreover, if $t_{k^{\prime}}\notin
T_{\mathrm{sym}}^{\ell}$, the entry at position
$(\xi_{2}(e_{j},t_{k^{\prime}}),\xi_{1}(e_{j},t_{k^{\prime}}))$ is equal to
$\frac{2}{3}$ and the entry at position
$(\xi_{2}(e_{j},t_{k^{\prime}}),\xi_{1}(e_{j},t_{k}))$ is equal to
$\frac{1}{3}$. Finally, for each pair
$(e_{j},t_{k}),(e_{j},t_{k^{\prime}})\in\Xi_{1}$ such that both $t_{k}$ and
$t_{k^{\prime}}$ are in $T^{\ell}\setminus
T_{\mathrm{sym}}^{\ell}\cup\\{\emptyset\\}$, the entries at positions
$(\xi_{2}(e_{j},t_{k}),\xi_{1}(e_{j},t_{k}))$ and
$(\xi_{2}(e_{j},t_{k^{\prime}}),\xi_{1}(e_{j},t_{k^{\prime}}))$ are equal to
$1$. All other entries of $H_{2}$ are equal to $0$. Note that if
$T_{\mathrm{sym}}^{\ell}$ is empty, $H_{2}$ is a permutation of the identity
matrix.

The mapping $\mathcal{B}_{2}:\Omega\rightarrow\mathbb{R}^{m_{2}}$ defined by

$\mathcal{B}_{2}(p)=H_{2}\cdot\mathcal{B}_{1}(p)$

determines a new set of functions, which we denote by
$B_{1}^{2},B_{2}^{2},\ldots,B_{m_{2}}^{2}$. Examples of them are shown in
Figure 4 (right). For all $t_{k}\in T_{\mathrm{sym}}^{\ell}$, the new
functions $B_{\xi_{2}(\emptyset,t_{k})}^{2}$ coincide with the basis functions
described by Grošelj and Speleers (2021).

######  Theorem 7.

The functions $B_{1}^{2},B_{2}^{2},\ldots,B_{m_{2}}^{2}$ possess the following
super-smoothness properties.

  1. 1.

For each $t_{k}\in T^{\ell}$, the functions are $C^{2}$ smooth at the split
point $v_{k}^{t}$.

  2. 2.

For each $t_{k}\in T^{\ell}$ and each edge $e_{j}$ of $t_{k}$, the functions
are $C^{2}$ smooth across the edge connecting $v_{j}^{e}$ and $v_{k}^{t}$.

  3. 3.

For each $t_{k}\in T_{\mathrm{sym}}^{\ell}$, the functions are $C^{2}$ smooth
over the interior of the triangle $t_{k}$.

Let $\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$ be the span of the
functions $B_{1},B_{2},\ldots,B_{m}$ and the functions
$B_{1}^{2},B_{2}^{2},\ldots,B_{m_{2}}^{2}$. The following properties can be
deduced by combining results from Grošelj and Speleers (2021) and Speleers
(2015).

######  Theorem 8.

The functions $B_{1},B_{2},\ldots,B_{m}$ and the functions
$B_{1}^{2},B_{2}^{2},\ldots,B_{m_{2}}^{2}$ form a convex partition of unity
and a locally supported basis of
$\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$.

######  Example 9.

Figure 6 shows a few rows and columns of the matrix $H_{2}$ which determine
the 6 reduced basis functions depicted in Figure 4 (right). Thanks to the
properties of the first reduced basis functions (see Theorem 4), this matrix
clearly confirms that the second reduced basis functions form a convex
partition of unity. Indeed, all its entries are nonnegative and the entries of
each column sum to one. The local supports of the second reduced basis
functions follow from the sparsity of the matrix.

Figure 6: An example of the second reduction. The basis functions shown in the
right picture of Figure 4 are produced from the basis functions shown in the
left picture of Figure 4 and a number of additional basis functions from the
first reduction. An excerpt from the matrix $H_{2}$ contains the six rows that
describe this reduction.

In Appendix A we show that the spline space
$\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$ contains all bivariate
polynomials of total degree less than or equal to $3$. This implies that the
space reduction presented here maintains the same optimal approximation power
as the other spline spaces
$\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})\subset\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell})$.

##  3 Dimension formulas

The three spline spaces considered in the previous section differ in local
super-smoothness, but all reproduce cubic polynomials and possess an optimal
approximation order of four. In this section, we will analyze and compare
their dimensions.

Let $\triangle$ be an arbitrary triangulation consisting of $n_{v}$ vertices,
$n_{e}$ edges, $n_{\partial e}$ boundary edges, and $n_{t}$ triangles. Recall
that $\triangle^{\ell}$, $\ell\in\mathbb{N}$, is obtained from $\triangle$ by
splitting each triangle of $\triangle$ uniformly such that each edge of
$\triangle$ is split into $\ell$ intervals. This refined triangulation has
$n_{v}^{\ell}$ vertices, $n_{e}^{\ell}$ edges, $n_{\partial e}^{\ell}$
boundary edges, and $n_{t}^{\ell}$ triangles, where

$\displaystyle n_{v}^{\ell}$
$\displaystyle=n_{v}+(\ell-1)n_{e}+\frac{1}{2}(\ell-2)(\ell-1)n_{t},$
$\displaystyle n_{e}^{\ell}$
$\displaystyle=\ell\,n_{e}+\frac{3}{2}(\ell-1)\ell\,n_{t},$ $\displaystyle
n_{\partial e}^{\ell}$ $\displaystyle=\ell n_{\partial e}=2\ell n_{e}-3\ell
n_{t},$ $\displaystyle n_{t}^{\ell}$ $\displaystyle=\ell^{2}\,n_{t}.$

Moreover, let $n_{\mathrm{sym}}^{\ell}$ be the number of triangles in
$T_{\mathrm{sym}}^{\ell}$. For simplicity of exposition, we assume that the
triangles in $\triangle$ do not form locally uniform subtriangulations.

The dimension of the full spline space
$\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell})$ is given by

$\dim(\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell}))=3n_{v}^{\ell}+4n_{e}^{\ell}=3n_{v}+(7\ell-3)n_{e}+\frac{3}{2}(5\ell-2)(\ell-1)n_{t}.$

Similarly, we deduce that the spline space after the first reduction
$\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})$ has dimension

$\dim(\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell}))=3n_{v}^{\ell}+2n_{e}^{\ell}=3n_{v}+(5\ell-3)n_{e}+\frac{3}{2}(3\ell-2)(\ell-1)n_{t}.$

Finally, the dimension of the second reduced spline space
$\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$ equals

$\displaystyle\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))$
$\displaystyle=3n_{v}^{\ell}$ (corresponding to vertices)
$\displaystyle\quad+n_{\mathrm{sym}}^{\ell}$ (corresponding to symmetric
triangles)
$\displaystyle\quad+3(n_{t}^{\ell}-n_{\mathrm{sym}}^{\ell})+n_{\partial
e}^{\ell},$ (corresponding to asymmetric triangles)

where the number $n_{\mathrm{sym}}^{\ell}$ depends on the geometry of
$\triangle_{\mathrm{PS}}^{\ell}$. In Example 10 we illustrate two extreme
cases of $n_{\mathrm{sym}}^{\ell}$.

Figure 7: Illustration of two extreme cases of $T_{\mathrm{sym}}^{\ell}$. The
marks and color codes are the same as in Figure 2.

######  Example 10.

If all the triangles in $\triangle$ are refined as shown in the left picture
of Figure 7, then $T_{\mathrm{sym}}^{\ell}=T_{0}^{\ell}\cup T_{1}^{\ell}$ and
the number of triangles in $T_{\mathrm{sym}}^{\ell}$ is maximal, namely

$n_{\mathrm{sym}}^{\ell}=[\ell^{2}-3(\ell-1)]n_{t},\quad\ell\geq 2.$

On the other hand, if all the triangles in $\triangle$ are refined as shown in
the right picture of Figure 7, then $T_{\mathrm{sym}}^{\ell}=T_{0}^{\ell}$ and
the number of triangles in $T_{\mathrm{sym}}^{\ell}$ is minimal, namely

$n_{\mathrm{sym}}^{\ell}=(\ell-3)^{2}n_{t},\quad\ell\geq 3.$

In general, we have

$(\ell-3)^{2}n_{t}\leq
n_{\mathrm{sym}}^{\ell}\leq[\ell^{2}-3(\ell-1)]n_{t},\quad\ell\geq 3.$

An illustration of a general configuration can be found in Figure 2.

In the following we fix $\ell\geq 3$. We can rewrite the above dimension
formula of $\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$ as

$\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))=3n_{v}+(5\ell-3)n_{e}+\frac{3}{2}(3\ell^{2}-5\ell+2)n_{t}-2n_{\mathrm{sym}}^{\ell},$

and from the bounds for $n_{\mathrm{sym}}^{\ell}$ provided in Example 10 we
obtain

$\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))\geq
3n_{v}+(5\ell-3)n_{e}+\frac{1}{2}(5\ell^{2}-3\ell-6)n_{t}$

and

$\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))\leq
3n_{v}+(5\ell-3)n_{e}+\frac{1}{2}(5\ell^{2}+9\ell-30)n_{t}.$

By comparing the different dimension formulas we deduce

$\displaystyle\dim(\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell}))-\dim(\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell}))$
$\displaystyle=2\ell\,n_{e}+3(\ell-1)\ell\,n_{t}>0,$
$\displaystyle\dim(\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell}))-\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))$
$\displaystyle\geq 2\ell\,n_{e}+(5\ell^{2}-15\ell+18)n_{t}>0.$

We also observe

$\displaystyle\lim_{\ell\rightarrow\infty}\frac{\dim(\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell}))}{\dim(\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell}))}$
$\displaystyle=\frac{15/2}{9/2}=\frac{5}{3},$
$\displaystyle\lim_{\ell\rightarrow\infty}\frac{\dim(\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell}))}{\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))}$
$\displaystyle=\frac{15/2}{5/2}=3.$

Thus, asymptotically speaking, the dimension drops by a factor $5/3$ and $3$,
respectively, for the two reduced spaces without loss of approximation order.

######  Remark 11.

The space of $C^{1}$ cubic Clough–Tocher splines is a well-known space over
triangulations; see Clough and Tocher (1965); Lai and Schumaker (2007). In
this case, the triangulation $\triangle^{\ell}$ is further refined to a
triangulation, denoted by $\triangle_{\mathrm{CT}}^{\ell}$, where each
triangle $t_{k}\in T^{\ell}$ is split into three smaller triangles by
connecting a split point $v_{k}^{t}$ inside $t_{k}$ to the vertices of
$t_{k}$. Let us denote the corresponding $C^{1}$ cubic spline space by
$\mathbb{S}_{0}(\triangle_{\mathrm{CT}}^{\ell})$. When the partitions
$\triangle_{\mathrm{CT}}^{\ell}$ and $\triangle_{\mathrm{PS}}^{\ell}$ are
compatible, it has been shown by Speleers (2015) that
$\mathbb{S}_{0}(\triangle_{\mathrm{CT}}^{\ell})$ is a subspace of
$\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})$ and, as a consequence, also
of $\mathbb{S}_{0}(\triangle_{\mathrm{PS}}^{\ell})$. The dimension of
$\mathbb{S}_{0}(\triangle_{\mathrm{CT}}^{\ell})$ equals

$\dim(\mathbb{S}_{0}(\triangle_{\mathrm{CT}}^{\ell}))=3n_{v}^{\ell}+n_{e}^{\ell}=3n_{v}+(4\ell-3)n_{e}+3(\ell-1)^{2}n_{t}.$

It is asymptotically slightly larger than the dimension of
$\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$. More precisely,

$\lim_{\ell\rightarrow\infty}\frac{\dim(\mathbb{S}_{0}(\triangle_{\mathrm{CT}}^{\ell}))}{\dim(\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell}))}=\frac{3}{5/2}=\frac{6}{5}.$

Moreover, there is no B-spline basis available for the space
$\mathbb{S}_{0}(\triangle_{\mathrm{CT}}^{\ell})$. Partial B-spline results can
be found in Speleers (2010b) and Lyche and Merrien (2018).

##  4 Applications

In this section, we provide numerical examples demonstrating the performance
of the considered spline spaces in the context of least squares approximation
and finite element approximation for second and fourth order boundary value
problems. We start off by explaining how to take advantage of the reductions
presented in Sections 2.3 and 2.4.

###  4.1 Practical aspects of basis reduction

Analogous to the mapping $\mathcal{B}_{0}$ defined in Section 2.2, let us
denote by $\mathcal{B}:\Omega\rightarrow\mathbb{R}^{m}$ the mapping that
assigns to the point $p\in\Omega$ the column vector
$[B_{1}(p),B_{2}(p),\ldots,B_{m}(p)]^{T}$. Suppose
$\mathcal{S}_{r}:\Omega\rightarrow\mathbb{R}^{m+m_{r}}$, $r\in\\{0,1,2\\}$, is
the mapping that assigns to $p\in\Omega$ the column vector obtained by
combining the vectors $\mathcal{B}(p)$ and $\mathcal{B}_{r}(p)$. Then,
$\mathcal{S}_{r}$ determines a basis of
$\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$, and for each
$s_{r}\in\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$, there exists a
column vector $c_{r}\in\mathbb{R}^{m+m_{r}}$ such that

$s_{r}(p)=\mathcal{S}_{r}(p)^{T}\cdot c_{r}.$ (1)

The mappings $\mathcal{S}_{r}$ can be related by extended extraction matrices
as

$\mathcal{S}_{1}(p)=\widetilde{H}_{1}\cdot\mathcal{S}_{0}(p),\quad\mathcal{S}_{2}(p)=\widetilde{H}_{2}\cdot\mathcal{S}_{1}(p),$

where

$\widetilde{H}_{1}=\begin{bmatrix}I&\\\
&H_{1}\end{bmatrix},\quad\widetilde{H}_{2}=\begin{bmatrix}I&\\\
&H_{2}\end{bmatrix},$

and $I\in\mathbb{R}^{m\times m}$ denotes the identity matrix. Moreover, for
$\alpha,\beta,\gamma\in\mathbb{N}_{0}$ such that $\alpha,\beta\leq\gamma$, let
$\mathcal{M}_{r}^{(\alpha,\beta,\gamma)}:\Omega\rightarrow\mathbb{R}^{(m+m_{r})\times(m+m_{r})}$
be the mappings defined by

$\mathcal{M}_{r}^{(\alpha,\beta,\gamma)}(p)=\left(D_{x}^{\alpha}D_{y}^{\gamma-\alpha}\mathcal{S}_{r}(p)\right)\cdot\left(D_{x}^{\beta}D_{y}^{\gamma-\beta}\mathcal{S}_{r}(p)\right)^{T}.$
(2)

Since

$\displaystyle\mathcal{M}_{1}^{(\alpha,\beta,\gamma)}(p)$
$\displaystyle=\widetilde{H}_{1}\cdot\mathcal{M}_{0}^{(\alpha,\beta,\gamma)}(p)\cdot\widetilde{H}_{1}^{T},$
$\displaystyle\mathcal{M}_{2}^{(\alpha,\beta,\gamma)}(p)$
$\displaystyle=\widetilde{H}_{2}\cdot\mathcal{M}_{1}^{(\alpha,\beta,\gamma)}(p)\cdot\widetilde{H}_{2}^{T},$

it is sufficient to compute $\mathcal{M}_{0}^{(\alpha,\beta,\gamma)}(p)$.
Then, $\mathcal{M}_{1}^{(\alpha,\beta,\gamma)}(p)$ and
$\mathcal{M}_{2}^{(\alpha,\beta,\gamma)}(p)$ are easy to obtain by matrix
multiplications. The considered mappings can be used to express matrices that
typically arise in approximation methods. For example, as we demonstrate in
the next sections, the mass and stiffness matrices can be acquired by
elementwise integration of these mappings.

###  4.2 Least squares approximation

For a function $f\in L^{2}(\Omega)$, the least squares approximation $s_{r}$
in $\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$ is given by the
orthogonality condition

$\int_{\Omega}\mathcal{S}_{r}(p)\cdot\left(f(p)-s_{r}(p)\right)\,\mathrm{d}p=0,$

where the integral is considered elementwise. By (1) and (2) this is the same
as

$\left(\int_{\Omega}\mathcal{M}_{r}^{(0,0,0)}(p)\,\mathrm{d}p\right)\cdot
c_{r}=\int_{\Omega}\mathcal{S}_{r}(p)\cdot f(p)\,\mathrm{d}p.$

######  Example 12.

We approximate the function

$f(x,y)=\sin\left(7\pi(1-x)(1-y)\right)$ (3)

on $\Omega=[0,1]\times[0,1]$; see Figure 8 (left). The domain is initially
partitioned by the triangulation $\triangle$ shown in Figure 9 and
subsequently refined by the triangulations $\triangle^{\ell}$ for
$\ell=2^{k}$, $k=1,2,3$. Figure 10 shows the $L^{2}$, $H^{1}$, and $H^{2}$
errors of the least squares approximations from the spaces
$\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$, $r=0,1,2$, with respect to
the square root of the numbers of degrees of freedom (NDOF) associated with
the spaces. As the square root of NDOF is inversely proportional to the length
of the longest edge in the triangulation, we observe that the convergence
rates are optimal in all three cases with respect to all three errors. We also
see that the solutions from the reduced spaces outperform the solutions from
the full space in terms of accuracy with respect to NDOF.

Figure 8: Graphs of the functions (3), (6), (8) (from left to right) used in
the numerical examples.

Figure 9: An unstructured triangulation $\triangle=\triangle^{1}$ of the unit
square (left) and three subsequent uniform refinements $\triangle^{\ell}$ for
$\ell=2^{k}$, $k=1,2,3$ (from left to right). The Powell–Sabin refinements of
$\triangle^{\ell}$ are depicted in light gray.

Figure 10: $L^{2}$, $H^{1}$, and $H^{2}$ errors (from left to right) with
respect to the square root of NDOF in the least squares approximation problem
of Example 12.

###  4.3 Approximate solution to a second order boundary value problem

For a function $f\in L^{2}(\Omega)$, an approximation $s_{r}$ from
$\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$ to the solution $u$ of the
second order boundary value problem

$\left\\{\begin{aligned} -\nabla^{2}u(p)&=f(p)&\text{if $p$}&\in\Omega,\\\
u(p)&=0&\text{if $p$}&\in\partial\Omega,\end{aligned}\right.$ (4)

is commonly obtained by considering the problem in the variational formulation

$\int_{\Omega}\nabla u(p)\cdot\nabla
v(p)\,\mathrm{d}p=\int_{\Omega}f(p)v(p)\,\mathrm{d}p,\quad v\in
H_{0}^{1}(\Omega).$ (5)

Suppose $\mathring{\mathcal{S}}_{r}$ is the column vector consisting of the
basis functions in $\mathcal{S}_{r}$ that are zero on the boundary of
$\Omega$. When discretizing $H_{0}^{1}(\Omega)$ with the space spanned by the
components of $\mathring{\mathcal{S}}_{r}$, the conditions in (5) can be
expressed as the system of linear equations

$\left(\int_{\Omega}\left(\mathring{\mathcal{M}}_{r}^{(1,1,1)}(p)+\mathring{\mathcal{M}}_{r}^{(0,0,1)}(p)\right)\,\mathrm{d}p\right)\cdot\mathring{c}_{r}=\int_{\Omega}\mathring{\mathcal{S}}_{r}(p)\cdot
f(p)\,\mathrm{d}p,$

where $\mathring{\mathcal{M}}_{r}^{(1,1,1)}$ and
$\mathring{\mathcal{M}}_{r}^{(0,0,1)}$ are defined in the same way as
$\mathcal{M}_{r}^{(1,1,1)}$ and $\mathcal{M}_{r}^{(0,0,1)}$ in (2) but with
$\mathcal{S}_{r}$ replaced by $\mathring{\mathcal{S}}_{r}$. The solution
$\mathring{c}_{r}$ of the system determines the approximation
$s_{r}=\mathring{\mathcal{S}}_{r}^{T}\cdot\mathring{c}_{r}$.

######  Example 13.

We approximate the manufactured solution

$u(x,y)=16x(1-x)y(1-y)\cos\left(16\pi\left((x-\tfrac{1}{2})^{2}+(y-\tfrac{1}{2})^{2}\right)\right),$
(6)

shown in Figure 8 (middle), to the problem (4) on the square domain
partitioned as in Example 12 (Figure 9). The $L^{2}$, $H^{1}$, and $H^{2}$
errors of the spline approximations from
$\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$ are provided in Figure 11. We
see again that the convergence rates are optimal and that the solutions from
the reduced spaces give better accuracy with respect to NDOF.

Figure 11: $L^{2}$, $H^{1}$, and $H^{2}$ errors (from left to right) with
respect to the square root of NDOF in the second order boundary value problem
of Example 13.

###  4.4 Approximate solution to a fourth order boundary value problem

Finally, for $f\in L^{2}(\Omega)$, we consider the fourth order boundary value
problem

$\left\\{\begin{aligned} \nabla^{4}u(p)&=f(p)&\text{if $p$}&\in\Omega,\\\
u(p)&=0&\text{if $p$}&\in\partial\Omega,\\\ \nabla u(p)\cdot\nu&=0&\text{if
$p$}&\in\partial\Omega,\end{aligned}\right.$ (7)

where $\nu$ is the unit outer normal to $\partial\Omega$. The solution $u$ can
be described through the variational formulation

$\int_{\Omega}\nabla^{2}u(p)\,\nabla^{2}v(p)\,\mathrm{d}p=\int_{\Omega}f(p)v(p)\,\mathrm{d}p,\quad
v\in H_{0}^{2}(\Omega).$

Let $\mathring{\mathcal{S}}_{r}$ be now the column vector consisting of the
basis functions in $\mathcal{S}_{r}$ that are zero on the boundary of $\Omega$
and whose normal derivatives are zero on the boundary of $\Omega$. After
discretizing $H_{0}^{2}(\Omega)$ with the space spanned by the components of
$\mathring{\mathcal{S}}_{r}$, we arrive at the system of linear equations

$\left(\int_{\Omega}\left(\mathring{\mathcal{M}}_{r}^{(2,2,2)}(p)+\mathring{\mathcal{M}}_{r}^{(2,0,2)}(p)+\mathring{\mathcal{M}}_{r}^{(0,2,2)}(p)+\mathring{\mathcal{M}}_{r}^{(0,0,2)}(p)\right)\,\mathrm{d}p\right)\cdot\mathring{c}_{r}=\int_{\Omega}\mathring{\mathcal{S}}_{r}(p)\cdot
f(p)\,\mathrm{d}p,$

where the mappings appearing in the integral on the left side are defined as
their analogous counterparts in (2) with $\mathcal{S}_{r}$ replaced by
$\mathring{\mathcal{S}}_{r}$. The solution $\mathring{c}_{r}$ of the system
specifies an approximation
$s_{r}=\mathring{\mathcal{S}}_{r}^{T}\cdot\mathring{c}_{r}$ to $u$ from
$\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$.

######  Example 14.

As in Examples 12 and 13, we consider the square domain partitioned by a
sequence of triangulations shown in Figure 9. Figure 12 shows the $L^{2}$,
$H^{1}$, and $H^{2}$ errors of the approximate solutions to the problem (7)
with the manufactured solution

$u(x,y)=\sin\left(2\pi(2x-y)\right)\sin(\pi x)^{4}\sin(\pi y)^{4}$ (8)

plotted in Figure 8 (right). The results indicate the same behavior of the
approximations from the spaces
$\mathbb{S}_{r}(\triangle_{\mathrm{PS}}^{\ell})$ as observed in the previous
examples. Indeed, they all exhibit optimal convergence rates, and the
reductions prove to be beneficial in achieving better accuracy with respect to
NDOF.

Figure 12: $L^{2}$, $H^{1}$, and $H^{2}$ errors (from left to right) with
respect to the square root of NDOF in the fourth order boundary value problem
of Example 14.

######  Remark 15.

The basis functions in $\mathring{\mathcal{S}}_{r}$ used for the fourth order
problem here are different from the ones used for the second order problem in
Section 4.3. Let us elaborate this in more detail for $r=2$. For each boundary
edge $e_{j}\in E^{\ell}$ and the corresponding triangle $t_{k}\in T^{\ell}$
that contains $e_{j}$, the basis function $B_{\xi_{2}(e_{j},t_{k})}^{2}$ has
zero values but not zero normal derivatives on the boundary. An example of
such a basis function is depicted in red in Figure 4 (right). Moreover, for
each boundary vertex $v_{i}\in V^{\ell}$ connecting two collinear boundary
edges and a proper choice of the triangle $q_{i}^{v}$ (with one side on the
boundary), one of the three corresponding basis functions also has zero values
but not zero normal derivatives on the boundary. An example of such a basis
function is depicted in blue in Figure 3 (left).

##  5 Conclusion

It is very common that the domain over which spline methods are used is in
large part partitioned regularly, i.e., only small regions are subjected to
irregular meshing. Even though often negligible in area, the irregular parts
are far from insignificant in spline constructions and may be, when not
properly handled, the cause of decreased accuracy and shape imperfections. In
this paper, we have chosen to address such a problem by applying an
unstructured spline technology and by reducing the number of degrees of
freedom associated with regularly refined regions through the imposition of
extra smoothness conditions.

The main contribution is an extraction process based on cubic Powell–Sabin
B-splines. Thanks to the linear independence, local support, $C^{1}$
continuity, non-negativity, partition of unity, and cubic precision, these
spline functions provide a sound and reliable approximation framework over any
given triangulation. The extraction of smoother functions from the initial set
of basis functions does not compromise the optimal approximation power, which
in numerical examples manifests in better accuracy with respect to the number
of degrees of freedom.

The extraction is presented as a two-step process. Firstly, a partial $C^{2}$
smoothness is gained inside triangles, i.e., over certain edges and split
points of the Powell–Sabin refinement. Secondly, a complete $C^{2}$ smoothness
is obtained over the interior of each triangle with symmetric configuration.
The transition between the full set of basis functions and the two reduced
alternatives is simple and is expressed in terms of sparse extraction
matrices. These can be directly applied in approximation methods, which is in
the present paper demonstrated in the context of the least squares
approximation method and the finite element method for solving second and
fourth order boundary value problems.

Although the paper considers the extraction process over an initial
unstructured triangulation that is subsequently uniformly refined, the
presented techniques are applicable whenever one can identify the set of
triangles with symmetric configuration. Moreover, it is likely that the notion
of symmetric configuration can be loosened to a more general geometric
symmetry when the edges of the Powell–Sabin refinement meet with only three
different slopes at the triangle split point. This extension, which requires a
different recombination of basis functions in order to preserve cubic
precision, is a subject for future research.

## Acknowledgements

J. Grošelj was partially supported by the research programme P1-0294 of Javna
agencija za raziskovalno dejavnost Republike Slovenije (ARRS). H. Speleers was
supported in part by a GNCS 2022 project (CUP E55F22000270001) of Gruppo
Nazionale per il Calcolo Scientifico – Istituto Nazionale di Alta Matematica
(GNCS – INdAM) and by the MUR Excellence Department Project MatMod@TOV (CUP
E83C23000330006) awarded to the Department of Mathematics of the University of
Rome Tor Vergata.

##  Appendix A Cubic polynomial reproduction

In this appendix, we show that the second reduced spline space
$\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$ contains the space
$\mathbb{P}_{3}$ of bivariate polynomials of total degree less than or equal
to $3$. The characterization of the first reduced spline space
$\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})$ given by Speleers (2015)
tells us that any $g\in\mathbb{P}_{3}$ belongs to
$\mathbb{S}_{1}(\triangle_{\mathrm{PS}}^{\ell})$. Considering (1), it can be
uniquely represented as

$g(p)=\mathcal{S}_{1}(p)^{T}\cdot c_{1}$

for some column vector $c_{1}\in\mathbb{R}^{m+m_{1}}$. Following the reduction
procedure from Section 2.4, we see that the basis functions
$B_{1},B_{2},\ldots,B_{m}$ remain unchanged in the basis $\mathcal{S}_{2}$.
The same holds true for the basis functions $B_{\xi_{1}(e_{j},t_{k})}^{1}$,
$B_{\xi_{1}(e_{j},t_{k^{\prime}})}^{1}$ for each pair
$(e_{j},t_{k}),(e_{j},t_{k^{\prime}})\in\Xi_{1}$ such that both $t_{k}$ and
$t_{k^{\prime}}$ are in $T^{\ell}\setminus
T_{\mathrm{sym}}^{\ell}\cup\\{\emptyset\\}$. Therefore, in order to verify
that $g(p)$ can be represented in terms of $\mathcal{S}_{2}(p)$, we only need
to focus on the remaining set of basis functions.

Let $t_{k}\in T_{\mathrm{sym}}^{\ell}$ be a triangle with vertices
$v_{i_{1}}$, $v_{i_{2}}$, and $v_{i_{3}}$. Moreover, let $t_{k^{\prime}}\in
T^{\ell}$ be the neighboring triangle with vertices $v_{i_{1}}$, $v_{i_{2}}$,
and $v_{i_{3}^{\prime}}$. Finally, let $e_{j}\in E^{\ell}$ be the edge between
the vertices $v_{i_{1}}$ and $v_{i_{2}}$. Since $t_{k}$ has symmetric
configuration, we know that

$\displaystyle v_{k}^{t}$
$\displaystyle=\frac{v_{i_{1}}+v_{i_{2}}+v_{i_{3}}}{3}=\frac{2v_{i_{3}}+v_{i_{3}^{\prime}}}{3},$
$\displaystyle v_{k^{\prime}}^{t}$
$\displaystyle=\frac{v_{i_{1}}+v_{i_{2}}+v_{i_{3}^{\prime}}}{3}=\frac{v_{i_{3}}+2v_{i_{3}^{\prime}}}{3}.$

Let $\mathfrak{B}[g]$ denote the blossom of the polynomial
$g\in\mathbb{P}_{3}$; see Seidel (1993); Speleers (2011) for more details
about blossoming. From Grošelj and Speleers (2017) we know that $g(p)$ can be
represented in terms of $\mathcal{S}_{1}(p)$ as

$g(p)=\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{k}^{t})B_{\xi_{1}(e_{j},t_{k})}^{1}(p)+\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{k^{\prime}}^{t})B_{\xi_{1}(e_{j},t_{k^{\prime}})}^{1}(p)+R_{1}(p),$

where $R_{1}(p)$ collects the contributions of the remaining basis functions
in $\mathcal{S}_{1}(p)$. By exploiting the multi-affinity property of
blossoms, we obtain

$\displaystyle g(p)$
$\displaystyle=\mathfrak{B}[g]\Bigl{(}v_{i_{1}},v_{i_{2}},\frac{2v_{i_{3}}+v_{i_{3}^{\prime}}}{3}\Bigr{)}B_{\xi_{1}(e_{j},t_{k})}^{1}(p)+\mathfrak{B}[g]\Bigl{(}v_{i_{1}},v_{i_{2}},\frac{v_{i_{3}}+2v_{i_{3}^{\prime}}}{3}\Bigr{)}B_{\xi_{1}(e_{j},t_{k^{\prime}})}^{1}(p)+R_{1}(p)$
$\displaystyle=\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{i_{3}})\Bigl{[}\frac{2}{3}B_{\xi_{1}(e_{j},t_{k})}^{1}(p)+\frac{1}{3}B_{\xi_{1}(e_{j},t_{k^{\prime}})}^{1}(p)\Bigr{]}$
$\displaystyle\quad+\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{i_{3}^{\prime}})\Bigl{[}\frac{1}{3}B_{\xi_{1}(e_{j},t_{k})}^{1}(p)+\frac{2}{3}B_{\xi_{1}(e_{j},t_{k^{\prime}})}^{1}(p)\Bigr{]}+R_{1}(p).$

After repeating the above argument for all pairs of triangles with at least
one triangle in $T_{\mathrm{sym}}^{\ell}$, we arrive at

$g(p)=\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{i_{3}})B_{\xi_{2}(\emptyset,t_{k})}^{2}(p)+\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{i_{3}^{\prime}})B_{\xi_{2}(\emptyset,t_{k^{\prime}})}^{2}(p)+R_{2}(p)$

if $t_{k^{\prime}}\in T_{\mathrm{sym}}$, or

$g(p)=\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{i_{3}})B_{\xi_{2}(\emptyset,t_{k})}^{2}(p)+\mathfrak{B}[g](v_{i_{1}},v_{i_{2}},v_{i_{3}^{\prime}})B_{\xi_{2}(e_{j},t_{k^{\prime}})}^{2}(p)+R_{2}(p)$

if $t_{k^{\prime}}\in T^{\ell}\setminus T_{\mathrm{sym}}^{\ell}$, where in
both cases $R_{2}(p)$ collects the contributions of the remaining basis
functions in $\mathcal{S}_{2}(p)$. This means that every polynomial
$g\in\mathbb{P}_{3}$ can be represented in the basis $\mathcal{S}_{2}$ and so
belongs to the space $\mathbb{S}_{2}(\triangle_{\mathrm{PS}}^{\ell})$.

## References

  * Barucq et al. (2022) Barucq, H., Calandra, H., Diaz, J., Frambati, S., 2022\.  Polynomial-reproducing spline spaces from fine zonotopal tilings.  J. Comput. Appl. Math. 402, 113812\.
  * Beirão da Veiga et al. (2015) Beirão da Veiga, L., Hughes, T.J.R., Kiendl, J., Lovadina, C., Niiranen, J., Reali, A., Speleers, H., 2015.  A locking-free model for Reissner–Mindlin plates: Analysis and isogeometric implementation via NURBS and triangular NURPS.  Math. Models Methods Appl. Sci. 25, 1519–1551.
  * Borden et al. (2011) Borden, M.J., Scott, M.A., Evans, J.A., Hughes, T.J.R., 2011\.  Isogeometric finite element data structures based on Bézier extraction of NURBS.  Int. J. Numer. Methods Eng. 87, 15–47.
  * Bressan and Sande (2019) Bressan, A., Sande, E., 2019\.  Approximation in FEM, DG and IGA: a theoretical comparison.  Numer. Math. 143, 923–942.
  * Cao et al. (2019) Cao, J., Chen, Z., Wei, X., Zhang, Y.J., 2019.  A finite element framework based on bivariate simplex splines on triangle configurations.  Comput. Methods Appl. Mech. Engrg. 357, 112598.
  * Chan et al. (2018) Chan, C.L., Anitescu, C., Rabczuk, T., 2018.  Isogeometric analysis with strong multipatch $C^{1}$-coupling.  Comput. Aided Geom. Design 62, 294–310.
  * Clough and Tocher (1965) Clough, R.W., Tocher, J.L., 1965\.  Finite element stiffness matrices for analysis of plates in bending, in: Conf. on Matrix Methods in Structural Mechanics, Wright–Patterson Air Force Base, Ohio. pp. 515–545.
  * Cottrell et al. (2009) Cottrell, J.A., Hughes, T.J.R., Bazilevs, Y., 2009.  Isogeometric Analysis: Toward Integration of CAD and FEA.  John Wiley & Sons.
  * Dierckx (1997) Dierckx, P., 1997.  On calculating normalized Powell–Sabin B-splines.  Comput. Aided Geom. Design 15, 61–78.
  * Evans et al. (2009) Evans, J.A., Bazilevs, Y., Babuska, I., Hughes, T.J.R., 2009\.  $n$-Widths, sup-infs, and optimality ratios for the $k$-version of the isogeometric finite element method.  Comput. Methods Appl. Mech. Engrg. 198, 1726–1741.
  * Grošelj and Knez (2021) Grošelj, J., Knez, M., 2021\.  A construction of edge B-spline functions for a $C^{1}$ polynomial spline on two triangles and its application to Argyris type splines.  Comput. Math. Appl. 99, 329–344.
  * Grošelj and Knez (2022) Grošelj, J., Knez, M., 2022\.  Generalized $C^{1}$ Clough–Tocher splines for CAGD and FEM.  Comput. Methods Appl. Mech. Engrg. 395, 114983.
  * Grošelj and Speleers (2017) Grošelj, J., Speleers, H., 2017\.  Construction and analysis of cubic Powell–Sabin B-splines.  Comput. Aided Geom. Design 57, 1–22.
  * Grošelj and Speleers (2021) Grošelj, J., Speleers, H., 2021\.  Super-smooth cubic Powell–Sabin splines on three-directional triangulations: B-spline representation and subdivision.  J. Comput. Appl. Math. 386, 113245\.
  * Hughes et al. (2008) Hughes, T.J.R., Reali, A., Sangalli, G., 2008.  Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of $p$-method finite elements with $k$-method NURBS.  Comput. Methods Appl. Mech. Engrg. 197, 4104–4124.
  * Jaxon and Qian (2014) Jaxon, N., Qian, X., 2014.  Isogeometric analysis on triangulations.  Comput. Aided Design 46, 45–57.
  * Kapl et al. (2019) Kapl, M., Sangalli, G., Takacs, T., 2019.  An isogeometric $C^{1}$ subspace on unstructured multi-patch planar domains.  Comput. Aided Geom. Design 69, 55–75.
  * Karčiauskas et al. (2016) Karčiauskas, K., Nguyen, T., Peters, J., 2016.  Generalizing bicubic splines for modeling and IGA with irregular layout.  Comput. Aided Design 70, 23–35.
  * Lai and Schumaker (2007) Lai, M.J., Schumaker, L.L., 2007\.  Spline Functions on Triangulations.  Cambridge University Press.
  * Liu and Jeffers (2018) Liu, N., Jeffers, A.E., 2018\.  A geometrically exact isogeometric Kirchhoff plate: Feature-preserving automatic meshing and $C^{1}$ rational triangular Bézier spline discretizations.  Int. J. Numer. Methods Eng. 115, 395–409.
  * Liu and Snoeyink (2007) Liu, Y., Snoeyink, J., 2007\.  Quadratic and cubic B-splines by generalizing higher-order Voronoi diagrams, in: Proc. of the 23rd Annual Symposium on Computational Geometry, ACM. pp. 150–157.
  * Lyche et al. (2022) Lyche, T., Manni, C., Speleers, H., 2022.  Construction of $C^{2}$ cubic splines on arbitrary triangulations.  Found. Comput. Math. 22, 1309–1350.
  * Lyche and Merrien (2018) Lyche, T., Merrien, J.L., 2018\.  Simplex-splines on the Clough–Tocher element.  Comput. Aided Geom. Design 65, 76–92.
  * Manni et al. (2022) Manni, C., Sande, E., Speleers, H., 2022.  Application of optimal spline subspaces for the removal of spurious outliers in isogeometric discretizations.  Comput. Methods Appl. Mech. Engrg. 389, 114260.
  * Neamtu (2007) Neamtu, M., 2007.  Delaunay configurations and multivariate splines: A generalization of a result of B. N. Delaunay.  Trans. Amer. Math. Soc. 359, 2993–3004.
  * Powell and Sabin (1977) Powell, M.J.D., Sabin, M.A., 1977\.  Piecewise quadratic approximations on triangles.  ACM Trans. Math. Softw. 3, 316–325.
  * Sande et al. (2019) Sande, E., Manni, C., Speleers, H., 2019.  Sharp error estimates for spline approximation: Explicit constants, $n$-widths, and eigenfunction convergence.  Math. Models Methods Appl. Sci. 29, 1175–1205.
  * Schmitt (2021) Schmitt, D., 2021.  Bivariate B-splines from convex configurations.  J. Comput. Syst. Sci. 120, 42–61.
  * Scott et al. (2011) Scott, M.A., Borden, M.J., Verhoosel, C.V., Sederberg, T.W., Hughes, T.J.R., 2011.  Isogeometric finite element data structures based on Bézier extraction of T-splines.  Int. J. Numer. Methods Eng. 88, 126–156.
  * Seidel (1993) Seidel, H.P., 1993.  An introduction to polar forms.  IEEE Comp. Graph. Appl. 13, 38–46.
  * Speleers (2010a) Speleers, H., 2010a.  A normalized basis for quintic Powell–Sabin splines.  Comput. Aided Geom. Design 27, 438–457.
  * Speleers (2010b) Speleers, H., 2010b.  A normalized basis for reduced Clough–Tocher splines.  Comput. Aided Geom. Design 27, 700–712.
  * Speleers (2011) Speleers, H., 2011.  On multivariate polynomials in Bernstein–Bézier form and tensor algebra.  J. Comput. Appl. Math. 236, 589–599.
  * Speleers (2013) Speleers, H., 2013.  Construction of normalized B-splines for a family of smooth spline spaces over Powell–Sabin triangulations.  Constr. Approx. 37, 41–72.
  * Speleers (2015) Speleers, H., 2015.  A new B-spline representation for cubic splines over Powell–Sabin triangulations.  Comput. Aided Geom. Design 37, 42–56.
  * Speleers (2019) Speleers, H., 2019.  Algorithm 999: Computation of multi-degree B-splines.  ACM Trans. Math. Softw. 45, 43\.
  * Speleers and Manni (2015) Speleers, H., Manni, C., 2015\.  Optimizing domain parameterization in isogeometric analysis based on Powell–Sabin splines.  J. Comput. Appl. Math. 289, 68–86.
  * Speleers et al. (2012) Speleers, H., Manni, C., Pelosi, F., Sampoli, M.L., 2012\.  Isogeometric analysis with Powell–Sabin splines for advection-diffusion-reaction problems.  Comput. Methods Appl. Mech. Engrg. 221–222, 132–148.
  * Speleers and Toshniwal (2021) Speleers, H., Toshniwal, D., 2021\.  A general class of $C^{1}$ smooth rational splines: Application to construction of exact ellipses and ellipsoids.  Comput. Aided Design 132, 102982\.
  * Toshniwal et al. (2017) Toshniwal, D., Speleers, H., Hughes, T.J.R., 2017.  Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations.  Comput. Methods Appl. Mech. Engrg. 327, 411–458.
  * Wang et al. (2022) Wang, Z., Cao, J., Wei, X., Zhang, Y.J., 2022.  TCB-spline-based isogeometric analysis method with high-quality parameterizations.  Comput. Methods Appl. Mech. Engrg. 393, 114771.
  * Wei et al. (2018) Wei, X., Zhang, Y.J., Toshniwal, D., Speleers, H., Li, X., Manni, C., Evans, J.A., Hughes, T.J.R., 2018\.  Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis.  Comput. Methods Appl. Mech. Engrg. 341, 609–639.
  * Weinmüller and Takacs (2022) Weinmüller, P., Takacs, T., 2022\.  An approximate $C^{1}$ multi-patch space for isogeometric analysis with a comparison to Nitsche’s method.  Comput. Methods Appl. Mech. Engrg. 401, 115592.
  * Zareh and Qian (2019) Zareh, M., Qian, X., 2019.  Kirchhoff–Love shell formulation based on triangular isogeometric analysis.  Comput. Methods Appl. Mech. Engrg. 347, 853–873.