# Tchebycheffian B-splines in isogeometric Galerkin methods

Krunal Raval<EMAIL_ADDRESS>Carla Manni<EMAIL_ADDRESS>Hendrik
Speleers<EMAIL_ADDRESS>Department of Mathematics, University of
Rome Tor Vergata, 00133 Rome, Italy

###### Abstract

Tchebycheffian splines are smooth piecewise functions whose pieces are drawn
from (possibly different) Tchebycheff spaces, a natural generalization of
algebraic polynomial spaces. They enjoy most of the properties known in the
polynomial spline case. In particular, under suitable assumptions,
Tchebycheffian splines admit a representation in terms of basis functions,
called Tchebycheffian B-splines (TB-splines), completely analogous to
polynomial B-splines. A particularly interesting subclass consists of
Tchebycheffian splines with pieces belonging to null-spaces of constant-
coefficient linear differential operators. They grant the freedom of combining
polynomials with exponential and trigonometric functions with any number of
individual shape parameters. Moreover, they have been recently equipped with
efficient evaluation and manipulation procedures. In this paper we consider
the use of TB-splines with pieces belonging to null-spaces of constant-
coefficient linear differential operators as an attractive substitute for
standard polynomial B-splines and rational NURBS in isogeometric Galerkin
methods. We discuss how to exploit the large flexibility of the geometrical
and analytical features of the underlying Tchebycheff spaces according to
problem-driven selection strategies. TB-splines offer a wide and robust
environment for the isogeometric paradigm beyond the limits of the rational
NURBS model.

###### keywords:

Isogeometric analysis, Tchebycheffian B-splines, B-splines with shape
parameters

††journal: Computer Methods in Applied Mechanics and Engineering

##  1 Introduction

Spline functions are ubiquitous in numerical methods. Besides their
theoretical interest, they have application in several branches of the
sciences including geometric modeling, signal processing, data analysis,
visualization and numerical simulation just to mention a few. Splines, in the
broad sense of the term, are functions consisting of pieces of smooth
functions glued together in a certain smooth way. There is a large variety of
spline species, often referred to as the zoo of splines. The most popular
species is the one where the pieces are algebraic polynomials of a given
degree $p$ and inter-smoothness is imposed by means of equality of derivatives
up to a certain order. Their popularity can be mainly attributed to their
representation in terms of the so-called B-splines. B-splines enjoy properties
such as local linear independence, minimal support, non-negativity, and
partition of unity; they can be computed through a stable recurrence relation
and can even be seen as the geometrically optimal basis for piecewise
polynomial spaces.

Tchebycheffian splines are smooth piecewise functions whose pieces are drawn
from (possibly different) extended Tchebycheff (ET-) spaces which are natural
generalizations of algebraic polynomial spaces [23, 36, 41]. Any non-trivial
element of an ET-space of dimension $p+1$ has at most $p$ zeros counting
multiplicity. We will refer to $p$ as the degree, in analogy with the
polynomial case. Extended complete Tchebycheff (ECT-) spaces are an important
subclass that can be generated through a set of positive weight functions [35,
41] and are spanned by generalized power functions [26]; the latter are the
natural extension of the monomial basis functions for algebraic polynomials;
see Section 2.1. Relevant examples are null-spaces of linear differential
operators on suitable intervals [41]; see Section 2.2. On bounded and closed
intervals the concepts of ET-space and ECT-space coincide [34], so from now on
we will focus on ECT-spaces.

Most of the results known for splines in the polynomial case extend in a
natural way to the Tchebycheffian setting. In particular, under suitable
assumptions on the involved ECT-spaces, Tchebycheffian splines admit a
representation in terms of basis functions, called Tchebycheffian B-splines
(TB-splines), with similar properties to polynomial B-splines. TB-splines were
introduced in 1968 by Karlin [22] using generalized divided differences. There
are several other ways to define them, including Hermite interpolation [11,
38], de Boor-like recurrence relations [17, 24], integral recurrence relations
[6], and blossoming [36]; see also the historical notes in [41, Chapters 9 and
11] for further details. Each of these definitions has advantages according to
the properties to be proved and lead to the same functions, up to a proper
scaling.

Multivariate extensions of Tchebycheffian splines can be easily obtained by
means of (local) tensor-product structures [7, 8].

Due to the richness of ECT-spaces, Tchebycheffian splines can find
applications in several contexts including data approximation/interpolation,
geometric modeling and numerical simulation; see [26] and references therein.
Thanks to their structural similarities, TB-splines are theoretically plug-to-
plug compatible with classical polynomial B-splines, so they can be
potentially easily incorporated in any software library supporting polynomial
B-splines to enrich its capability.

Unfortunately, despite their theoretical interest and applicative potential,
TB-splines have not gained much attention in practice so far. The reason
behind this is that TB-splines are generally difficult to compute. The
classical approaches mentioned above, based on generalized divided
differences, Hermite interpolation, or repeated integration, are
computationally expensive and/or numerically unstable. An important step
forward was recently made in [19] where the authors proposed a strategy that
represents TB-splines as linear combinations of local Tchebycheffian Bernstein
functions through a suitable extraction operator. The local Tchebycheffian
Bernstein functions form a basis of the local ECT-spaces involved in the
definition of the TB-splines. In the polynomial case, these are nothing but
the classical Bernstein polynomial basis functions. Following the approach in
[19], an object-oriented Matlab toolbox has been developed in [43] for the
construction and manipulation of TB-splines whenever they exist. The toolbox
supports TB-splines whose pieces belong to ECT-spaces that are null-spaces of
constant-coefficient linear differential operators, and is publicly available.
Note that both [19] and [43] address the more general setting of multi-degree
TB-splines (MDTB-splines) where the local ECT-spaces of the Tchebycheffian
spline space are not required to be of the same dimension. Having at our
disposal a publicly available implementation for such a large class of TB-
splines paves the path for their effective use in practical applications. An
alternative strategy has been proposed in [5] but no implementation is
available.

Isogeometric analysis (IgA) is a well established methodology for the analysis
of problems governed by partial differential equations (PDEs); see [15, 20].
IgA aims to simplify the interoperability between geometric modeling and
numerical simulation by constructing a fully integrated framework for
computer-aided design (CAD) and finite element analysis (FEA). In its original
formulation, IgA uses B-splines, or their rational generalization (NURBS),
both to represent the computational domain and to approximate the solution of
the PDE that models the problem of interest. The isogeometric paradigm shows
important advantages over classical FEA: complicated geometries can be
represented more accurately and some common profiles, such as conic sections,
are described exactly; the precise geometry description incorporated at the
coarsest mesh level can be maintained during mesh refinement through the
process of knot insertion; and the inherently high smoothness of
B-splines/NURBS allows for a higher accuracy per degree of freedom.
Nevertheless, the above advantages are not a distinguishing property of
B-splines/NURBS, and B-splines/NURBS are not a requisite ingredient in IgA.

An interesting subclass of TB-splines, the so-called generalized polynomial
B-splines (GB-splines), has been proposed as an alternative to classical
polynomial and rational splines in the context of IgA; see [28, 32] and
references therein. GB-splines can be seen as the minimal extension of
polynomial B-splines towards the wide variety of TB-splines: their pieces
belong to ECT-spaces obtained by enriching an algebraic polynomial space with
a pair of functions, typically hyperbolic or trigonometric functions
identified by a single shape parameter. GB-splines overcome two issues of the
rational NURBS model: they allow for an exact (almost) arc-length
parameterization of conic sections and behave with respect to differentiation
and integration as nicely as polynomial B-splines (for instance, the
derivative of a trigonometric generalized spline with a given phase parameter
$\beta$ and degree $p$ is again a trigonometric generalized spline with the
same shape parameter $\beta$ and degree $p-1$). These properties make GB-
splines an interesting tool to face both geometrical and analytical hurdles in
IgA. The effectiveness of hyperbolic or trigonometric GB-splines in
isogeometric Galerkin and collocation methods has been illustrated in a
sequence of papers, where their properties have been exploited to obtain exact
representations of common geometries [1, 28, 30] or to beneficially deal with
advection-dominated problems [29, 31]; see also [32] for the related spectral
properties.

Being a so minimal extension of the polynomial setting, however, GB-splines
are not always flexible enough for practical applications. In particular, in a
tensor-product GB-spline space only two hyperbolic or trigonometric functions
identified by the same shape parameter are added to polynomials along each
parametric direction. Therefore, a given tensor-product GB-spline space does
not allow for:

  * 1.

an exact representation of different arcs of conic sections at opposite sides;

  * 2.

a proper treatment of different analytic features (like sharp layers) along a
given parametric direction;

  * 3.

a simultaneous treatment of geometrical and analytical features along the same
parametric direction.

Tchebycheffian splines with pieces belonging to ECT-spaces that are null-
spaces of constant-coefficient linear differential operators are an extension
of hyperbolic and trigonometric generalized polynomial splines. They enjoy all
the geometrical and analytical features that motivate the interest in GB-
splines without suffering from the above mentioned restrictions. They grant
the freedom of combining polynomials with exponential and trigonometric
functions with any number of individual shape parameters. For this class of
Tchebycheffian splines, when the various pieces are drawn from a single ECT-
space which contains constants, the existence of TB-splines is always ensured,
possibly with some restriction on the partition; see Section 2.4 for more
details. Furthermore, they are supported by the Matlab toolbox available in
[43]. In summary, TB-splines with pieces belonging to ECT-spaces that are
null-spaces of constant-coefficient linear differential operators offer a
suitable balance between the immense variety of ECT-spaces and the practical
needs of a problem-driven space selection and efficient evaluation procedures
for the space elements. They can be seen as a flexible environment for the IgA
paradigm, beyond the limits of the rational NURBS model. In this paper we
present and discuss the use of such TB-splines in Galerkin isogeometric
methods.

The remainder of the paper consists of four sections. In the next section we
introduce ECT-spaces and Tchebycheffian spline spaces. We also briefly discuss
the construction and the properties of the corresponding TB-spline basis.
Section 3 is devoted to the use of TB-splines in isogeometric analysis and to
describe a general road map for the selection of suitable ECT-spaces in a
problem-driven fashion. A bunch of case studies illustrating the performance
of the proposed approach are presented in Section 4. Section 5 collects some
concluding remarks and outlines possible future research lines.

##  2 TB-spline theory

In this section we outline the mathematical framework for the spaces we are
going to use. We recall from the literature the definition of ECT-spaces and
we summarize some of their main properties. We also discuss a particular
subclass of ECT-spaces identified by linear differential operators with real
constant coefficients that will be our target class for isogeometric analysis.
Afterwards, we introduce Tchebycheffian splines – piecewise functions where
the pieces belong to given ECT-spaces glued together with certain smoothness –
and we briefly discuss the construction of TB-splines, which are B-spline like
basis functions for Tchebycheffian splines. We refer the reader to [25, 26,
36, 41] for in-depth information about Tchebycheffian splines and spline
functions in general.

The presentation of this section is somehow technical for the sake of
completeness. However, all the technical details can be skipped at a first
reading, focusing on the features of the (spline) spaces we are going to
describe.

###  2.1 Extended complete Tchebycheff spaces

We consider a class of $(p+1)$-dimensional spaces that are a natural extension
of the space of algebraic polynomials of degree less than or equal to $p$.

######  Definition 1 (Extended Tchebycheff space).

Given an interval $J$ and an integer $p\geq 0$, a space
$\mathbb{T}_{p}(J)\subset C^{p}(J)$ of dimension $p+1$ is an extended
Tchebycheff (ET-) space on $J$, if any Hermite interpolation problem with
$p+1$ data on $J$ has a unique solution in $\mathbb{T}_{p}(J)$. More in
details, let $z_{1},z_{2},\ldots,z_{\bar{m}}$ be distinct points on $J$, with
any positive integer $\bar{m}$ and $d_{1},d_{2},\ldots,d_{\bar{m}}$ be non-
negative integers such that $p+1=\sum_{i=1}^{\bar{m}}(d_{i}+1)$. Then for any
$f_{i,j}\in\mathbb{R}$, there exists a unique $g\in\mathbb{T}_{p}(J)$ such
that

$D^{j}g(z_{i})=f_{i,j},\quad j=1,\ldots,d_{i},\quad i=1,\ldots,\bar{m}.$

As an immediate consequence of the previous definition, we have that any non-
trivial element of an ET-space on $J$ has at most $p$ roots in $J$, counting
multiplicity. This property shows the strong interconnection between ET-spaces
and the space of algebraic polynomials.

We are now going to identify an interesting subclass of ET-spaces that
strengthens the similarity with algebraic polynomials even further. We do this
by generalizing the monomial basis functions, which form the so-called power
basis of the polynomial space. We denote by $\langle
g_{0},\ldots,g_{k}\rangle$ the linear space spanned by
$\\{g_{0},\ldots,g_{k}\\}$.

######  Definition 2 (Extended complete Tchebycheff space).

Given an interval $J$ and an integer $p\geq 0$, a space
$\mathbb{T}_{p}(J)\subset C^{p}(J)$ of dimension $p+1$ is an extended complete
Tchebycheff (ECT-) space on $J$, if there exist functions $g_{0},\ldots,g_{p}$
such that $\mathbb{T}_{p}(J)=\langle g_{0},\ldots,g_{p}\rangle$ and every
subspace $\langle g_{0},\ldots,g_{k}\rangle$ for $k=0,\ldots,p$ is an ET-space
on $J$. The basis $\big{\\{}g_{0},\ldots,g_{p}\big{\\}}$ is called an ECT-
system.

######  Example 3.

The space of algebraic polynomials $\mathbb{P}_{p}=\langle
1,x,\ldots,x^{p}\rangle$ is an ECT-space on the real line. Indeed, with any
fixed point $z\in\mathbb{R}$, it can be accounted as the span of the ECT-
system

$\Bigg{\\{}1,x-z,\frac{(x-z)^{2}}{2},\ldots,\frac{(x-z)^{p}}{p!}\Bigg{\\}}.$
(1)

It is evident from Definition 2 that every ECT-space on the interval $J$ is an
ET-space of the same dimension on the same interval. The contrary case is not
always true. However, any ET-space on a bounded and closed interval is an ECT-
space on the same interval [34]. In the scope of this paper, we are only
interested in bounded closed intervals. Hence, the notions of ET-space and
ECT-space are interchangeable. Throughout the paper, we consider ECT-spaces to
allow for a neat construction of a proper basis of related spline spaces. In
this perspective, we recall some additional properties of ECT-spaces.

A space $\mathbb{T}_{p}(J)\subset C^{p}(J)$ of dimension $p+1$ is an ECT-space
on $J$ if and only if there exist positive functions $w_{j}\in C^{p-j}(J)$,
$j=0,1,\ldots,p$ such that $\mathbb{T}_{p}(J)$ is spanned by the following
functions [26, 34, 41]:

$\begin{cases}g_{0}(x)\coloneqq w_{0}(x),\\\ g_{1}(x)\coloneqq
w_{0}(x)\int_{z}^{x}w_{1}(y_{1}){\rm d}y_{1},\\\ \hskip 31.2982pt\vdots\\\
g_{p}(x)\coloneqq
w_{0}(x)\int_{z}^{x}w_{1}(y_{1})\int_{z}^{y_{1}}\cdots\int_{z}^{y_{p-1}}w_{p}(y_{p}){\rm
d}y_{p}{\rm d}y_{p-1}\cdots{\rm d}y_{1},\end{cases}$ (2)

where $z$ is any point in $J$. The functions $g_{0},\ldots,g_{p}$ are called
generalized power functions, the functions $w_{0},\ldots,w_{p}$ are usually
referred to as weights, and the set $\\{w_{0},\ldots,w_{p}\\}$ is called a
weight system generating $\mathbb{T}_{p}(J)$.

######  Example 4.

A weight system for the space of algebraic polynomials $\mathbb{P}_{p}$ is
given by $\\{w_{j}=1:j=0,\ldots,p\\}$. Through (2) these weights provide the
basis (1).

######  Example 5.

Let $1\leq\ell\leq p$ and $\alpha_{1},\ldots,\alpha_{\ell}\in\mathbb{R}$ with
$\alpha_{i}\neq\alpha_{j}$ for each $i\neq j$. Then, we consider the functions

$w_{0}(x)=\dots=w_{p-\ell}(x)=1,\quad w_{p-\ell+1}(x)=e^{\alpha_{1}x},\quad
w_{p-\ell+k}(x)=e^{(\alpha_{k}-\alpha_{k-1})x},\quad k=2,\ldots,\ell.$ (3)

For any $z\in\mathbb{R}$, the identified generalized power functions in (2)
span the space

$\Big{\langle}{1,x,\ldots,x^{p-\ell},e^{\alpha_{1}x},\ldots,e^{\alpha_{\ell}x}}\Big{\rangle}.$

Since the functions in (3) are positive on $\mathbb{R}$, the above space is a
$(p+1)$-dimensional ECT-space on $\mathbb{R}$ for any $p\geq 1$.

######  Example 6.

For $\beta\in\mathbb{R}$, $\beta>0$, and $x\in J$ with

$J=\left(-\frac{\pi}{2\beta},\frac{\pi}{2\beta}\right),$

we consider the functions

$w_{0}(x)=\dots=w_{p-2}(x)=1,\quad w_{p-1}(x)=\cos(\beta x),\quad
w_{p}(x)=\frac{1}{\cos^{2}(\beta x)},\quad p\geq 1.$ (4)

For any $z\in J$, the identified generalized power functions in (2) span the
space

$\Big{\langle}{1,x,\ldots,x^{p-2},\cos(\beta x),\sin(\beta x)}\Big{\rangle}.$
(5)

Since the functions in (4) are positive on $J$, the above space is a
$(p+1)$-dimensional ECT-space on $J$ for any $p\geq 1$. The space in (5) is
often referred to as cycloidal space. Actually, for $p\geq 2$, it can be shown
that the space in (5) is an ECT-space on any interval of length less than
$2\pi/\beta$; see [14] and also [26, Section 2].

######  Remark 7.

A given ECT-space on an interval $J$ can be identified by different systems of
weights; see [26, Section 2] for details and examples. In particular, it is
clear that the two weight systems

$w_{0},\ldots,w_{p}\quad{\rm and}\quad K_{0}w_{0},\ldots,K_{p}w_{p},$

where $K_{0},\ldots,K_{p}$ are positive constants, identify the same ECT-
space.

######  Remark 8.

Let $\mathbb{T}_{p}$ be a space of dimension $p+1$ on $J$. Assume that the
derivative space of $\mathbb{T}_{p}$, i.e., the space spanned by the
derivatives of functions in $\mathbb{T}_{p}$, is a $p$-dimensional ECT-space
on the interval $J$ generated by the weights $\\{w_{1},\ldots,w_{p}\\}$. Then,
$\mathbb{T}_{p}$ is a $(p+1)$-dimensional ECT-space on $J$ generated by the
weights $\\{1,w_{1},\ldots,w_{p}\\}$.

Many properties of polynomials involve working with derivatives. When dealing
with ECT-spaces generated by the weight system $\\{w_{0},\ldots,w_{p}\\}$ it
is often convenient to replace the usual derivatives by some related
differential operators. We define $D_{0}f\coloneqq f$ and

$D_{j}f\coloneqq D\left(\frac{f}{w_{j-1}}\right),\quad j=1,\ldots,p+1,$

and we set

$L_{j}\coloneqq D_{j}D_{j-1}\cdots D_{0},\quad j=0,\ldots,p+1.$ (6)

The operator $L_{j}$ can be seen as a natural substitute for $D^{j}$; see [41,
Section 9.1]. In particular, we have

$L_{j}g_{k}(x)=\begin{cases}w_{j}(x)\int_{z}^{x}w_{j+1}(y_{j+1})\int_{z}^{y_{j+1}}\cdots\int_{z}^{y_{k-1}}w_{k}(y_{k}){\rm
d}y_{k}{\rm d}y_{k-1}\cdots{\rm d}y_{j+1},&k=j,\ldots,p,\\\
0,&k=0,\ldots,j-1.\end{cases}$

where $g_{k}$ are defined in (2). This confirms that the generalized power
functions are the natural extension of the monomial basis functions for
algebraic polynomials.

###  2.2 ECT-spaces and null-spaces of linear differential operators

A large and interesting class of ECT-spaces is given by the null-spaces of
linear differential operators with real constant coefficients. More precisely,
let us consider the differential operator defined by

$\mathcal{L}_{p}f\coloneqq D^{p+1}f+\sum_{j=0}^{p}a_{j}D^{j}f,\quad f\in
C^{p+1}(\mathbb{R}),\quad a_{j}\in\mathbb{R},\quad j=0,\ldots,p,$ (7)

and let

$\mathfrak{p}_{p}(\omega)\coloneqq\omega^{p+1}+\sum_{j=0}^{p}a_{j}\omega^{j}$
(8)

be its characteristic polynomial and $\mathbb{N}_{p}$ its null-space.

It is well known that the null-space of the differential operator (7) is
easily described by means of the roots of the characteristic polynomial (8).
Let $\omega=\alpha+\mathrm{i}\beta$ be a root of multiplicity $\mu\geq 1$ of
the polynomial in (8) for $\alpha,\beta\in\mathbb{R}$ and
$\mathrm{i}\coloneqq\sqrt{-1}$. Then, this root generates the following
fundamental subspace:

  * 1.

if $\beta=0$, then

$\Big{\langle}{x^{k}e^{\alpha
x}:k=0,\ldots,\mu-1}\Big{\rangle}\subseteq\mathbb{N}_{p};$

  * 2.

if $\beta\neq 0$, then the complex conjugate of $\omega$ is also a root of
multiplicity $\mu$, and

$\Big{\langle}{x^{k}e^{\alpha x}\cos(\beta x),x^{k}e^{\alpha x}\sin(\beta
x):k=0,\ldots,\mu-1}\Big{\rangle}\subseteq\mathbb{N}_{p}.$

The above subspaces for all roots together form the null-space
$\mathbb{N}_{p}$. To ensure that the constants belong to $\mathbb{N}_{p}$, we
have to assume that $\omega=0$ is a root of the characteristic polynomial in
(8).

######  Example 9.

Setting $a_{j}=0$, $j=0,\ldots,p$ in (7) gives the linear differential
operator $\mathcal{L}_{p}f=D^{p+1}f$. In this case, the null-space is the
space of algebraic polynomials $\mathbb{N}_{p}=\mathbb{P}_{p}$; see also
Example 3.

######  Remark 10.

If the polynomial (8) has only real roots, then $\mathbb{N}_{p}$ is an ECT-
space on the whole real line, and in particular on any bounded interval
$[a,b]$. On the other hand, if the characteristic polynomial has also complex
roots, then $\mathbb{N}_{p}$ is an ECT-space on sufficiently small intervals;
see [41]. More precisely, $\mathbb{N}_{p}$ is an ECT-space on any interval
$[a,b]$ such that

$b-a<\mathfrak{l}_{p}.$

The value $\mathfrak{l}_{p}$ is referred to as critical length and it can be
bounded from below as

$\mathfrak{l}_{p}\geq\pi/B>0,$ (9)

where $B$ denotes the maximum imaginary part of all non-real roots of the
characteristic polynomial. For a detailed study about the critical lengths of
such ECT-spaces, the reader is referred to [4] and [13].

In the rest of the paper we will mainly focus on $(p+1)$-dimensional ECT-
spaces that are null-spaces of linear differential operators with real
constant coefficients, as in (7). The reference interval $J$ ensuring the ECT-
property will be specified when necessary. We will always assume that
constants belong to the considered ECT-spaces. Let us now introduce the
notation we are going to use for such spaces and give some examples.

Consider linear differential operators whose characteristic polynomial has
roots $\omega_{k}=\alpha_{k}+\mathrm{i}\beta_{k}$ of multiplicity $\mu_{k}\geq
1$, with $k=0,\ldots,M$, such that $\sum_{k=0}^{M}\mu_{k}=p+1$, and
$\omega_{0}=0$. Under this assumption, the corresponding null-space is
uniquely characterized by the following vector with $p+1-\mu_{0}$ components:

${\scriptstyle{\mathcal{W}}}\coloneqq\Big{(}\underbrace{\omega_{1},\,\ldots,\,\omega_{1}}_{\mu_{1}\text{
times}},\,\underbrace{\omega_{2},\,\ldots,\,\omega_{2}}_{\mu_{2}\text{
times}},\,\ldots\,,\underbrace{\omega_{M},\,\ldots,\,\omega_{M}}_{\mu_{M}\text{
times}}\Big{)}.$ (10)

To highlight the dependence of the null-space on the roots (10), we use the
notation

$\mathbb{P}_{p}^{\scriptstyle{\mathcal{W}}}=\mathbb{P}_{p}^{(\alpha_{1}+\mathrm{i}\beta_{1},\ldots,\alpha_{1}+\mathrm{i}\beta_{1},\alpha_{2}+\mathrm{i}\beta_{2},\ldots,\alpha_{2}+\mathrm{i}\beta_{2},\ldots,\alpha_{M}+\mathrm{i}\beta_{M})}.$
(11)

The ECT-spaces constructed from the null-spaces of linear differential
operators with real constant coefficients offer great flexibility through the
shape parameters to be chosen in (10). A particular selection of these shape
parameters results into familiar subclasses of ECT-spaces. We now list some
popular examples.

  * 1.

The class of algebraic polynomial spaces is the most established subclass (see
Example 3):

$\mathbb{P}_{p}=\Big{\langle}{1,x,\ldots,x^{p}}\Big{\rangle}.$

  * 2.

Generalized polynomial spaces, see [26, 28], enrich algebraic polynomial
spaces by a pair of functions. An important class of generalized polynomial
spaces can be obtained as null-spaces of linear differential operators with
constant coefficients. A first example is the algebraic polynomial space of
degree $p-2$ enriched with two exponential functions. Taking
$\omega_{1}=\alpha$, $\omega_{2}=-\alpha$, and $\mu_{1}=\mu_{2}=1$ gives the
following null-space:

$\mathbb{P}_{p}^{(\alpha,-\alpha)}=\Big{\langle}{1,x,\ldots,x^{p-2},e^{\alpha
x},e^{-\alpha x}}\Big{\rangle}=\Big{\langle}{1,x,\ldots,x^{p-2},\cosh(\alpha
x),\sinh(\alpha x)}\Big{\rangle},\quad p\geq 2;$ (12)

see also Example 5. This is an ECT-space on any interval according to Remark
10. Another instance of generalized polynomial spaces is the combination of
the algebraic polynomial space of degree $p-2$ with two trigonometric
functions sharing the same phase $\beta>0$. Taking
$\omega_{1}=\mathrm{i}\beta$, $\omega_{2}=-\mathrm{i}\beta$, and
$\mu_{1}=\mu_{2}=1$ gives the following null-space:

$\mathbb{P}_{p}^{(\mathrm{i}\beta,-\mathrm{i}\beta)}=\Big{\langle}{1,x,\ldots,x^{p-2},\cos(\beta
x),\sin(\beta x)}\Big{\rangle},\quad p\geq 2,$ (13)

which is a cycloidal space; see also Example 6. Remark 10 provides a lower
bound on its critical length, while it can be bounded from above by

$\mathfrak{l}_{p}\leq\frac{2\pi}{\beta}\left\lfloor\frac{p}{2}\right\rfloor,\quad
p\geq 2;$

see [13]. A detailed and precise study on the critical length of cycloidal
spaces can be found in [14].

  * 3.

An interesting subclass of ECT-spaces obtained as null-spaces of linear
differential operators are the following exponential and trigonometric spaces,
both defined for $p=2q$:

$\displaystyle\mathbb{P}_{p}^{(\alpha,-\alpha,\ldots,q\alpha,-q\alpha)}$
$\displaystyle=\Big{\langle}{1,e^{\alpha x},e^{-\alpha x},e^{2\alpha
x},e^{-2\alpha x},\ldots,e^{q\alpha x},e^{-q\alpha x}}\Big{\rangle}$
$\displaystyle=\Big{\langle}{1,\cosh(\alpha x),\sinh(\alpha x),\cosh(2\alpha
x),\sinh(2\alpha x),\ldots,\cosh(q\alpha x),\sinh(q\alpha x)}\Big{\rangle},$

and

$\mathbb{P}_{p}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\ldots,\mathrm{i}q\beta,-\mathrm{i}q\beta)}=\Big{\langle}{1,\cos(\beta
x),\sin(\beta x),\ldots,\cos(q\beta x),\sin(q\beta x)}\Big{\rangle},$

identified by the vector (10) with $\omega_{k}=\pm k\alpha$ and
$\omega_{k}=\pm\mathrm{i}k\beta$, $k=1,\ldots,q$, respectively. In the former
case, the space is an ECT-space on any interval, while in the latter case, the
space is an ECT-space on any interval of length less than $\frac{\pi}{\beta}$
for $\beta>0$; see [33, 39].

All the above examples are special instances of the following subclass of ECT-
spaces that combines polynomial, exponential, and trigonometric functions. Let
$0\leq\ell\leq p$ and $\alpha_{1},\ldots,\alpha_{\ell}\in\mathbb{R}$ with
$\alpha_{i}\neq\alpha_{j}$ for each $i\neq j$. Moreover, let $0\leq 2q\leq
p-\ell$ and $\beta_{1},\ldots,\beta_{q}\in\mathbb{R}$ with
$\beta_{i}\neq\beta_{j}$ for each $i\neq j$. Then, let us consider the space

$\mathbb{P}_{p}^{(\alpha_{1},\ldots,\alpha_{\ell},\mathrm{i}\beta_{1},-\mathrm{i}\beta_{1},\ldots,\mathrm{i}\beta_{q},-\mathrm{i}\beta_{q})}=\Big{\langle}{1,x,\ldots,x^{p-\ell-2q},e^{\alpha_{1}x},\ldots,e^{\alpha_{\ell}x},\cos(\beta_{1}x),\sin(\beta_{1}x),\ldots,\cos(\beta_{q}x),\sin(\beta_{q}x)}\Big{\rangle},$
(14)

defined on a reference interval ensuring the ECT-property. Such space is
invariant under translation of the interval. We are going to utilize this
comprehensive subclass of ECT-spaces in the case studies presented in Section
4. In the following we will refer to $p$ as the degree of the ECT-space.

######  Remark 11.

The subclass of ECT-spaces in (14) is a proper subclass of null-spaces of
differential operators with real constant coefficients, see (11), but still
grants the freedom of combining polynomials with exponential and trigonometric
functions with any number of individual shape parameters to be selected
according to (automatic) problem-driven strategies. Moreover, derivatives and
integrals of functions belonging to spaces of the form (14) belong to spaces
of the same form, possibly not including constants anymore.

###  2.3 Tchebycheffian spline spaces

In a complete similarity with polynomial splines, Tchebycheffian splines are
piecewise functions with pieces belonging to ECT-spaces which can be glued
together with certain prescribed smoothness. More precisely, let $\Delta$ be a
partition of the interval $[a,b]\subset\mathbb{R}$ specified in terms of a
sequence of breakpoints,

$\Delta\coloneqq\\{a=:x_{0}<x_{1}<\dots<x_{m-1}<x_{m}\coloneqq b\\}.$

We set $J_{i}\coloneqq[x_{i-1},x_{i})$, $i=1,\ldots,m-1$, and
$J_{m}\coloneqq[x_{m-1},x_{m}]$. Furthermore, for $i=1,\ldots,m$ let
$\mathbb{T}_{p,i}$ be an ECT-space of dimension $p+1$ on the closed interval
$[x_{i-1},x_{i}]$. Finally, let us consider a sequence of $m-1$ integers

$\mathbf{r}\coloneqq\Big{\\{}r_{i}\in\mathbb{Z}:-1\leq r_{i}\leq
p-1,\,i=1,\ldots,m-1\Big{\\}}.$

The elements of the set

$\begin{split}\mathbb{S}_{p}^{\mathbf{r}}(\Delta)\coloneqq\Big{\\{}f:[a,b]\rightarrow\mathbb{R}:\,&f|_{J_{i}}\in\mathbb{T}_{p,i},\
i=1,\ldots,m;\\\ &D^{l}_{-}f(x_{i})=D^{l}_{+}f(x_{i}),\ l=0,\ldots,r_{i},\
i=1,\ldots,m-1\Big{\\}}\end{split}$ (15)

are Tchebycheffian splines of degree $p$ and smoothness $\mathbf{r}$ on the
partition $\Delta$ of the interval $[a,b]$. The dimension of the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ is given by

$n\coloneqq p+1+\sum_{i=1}^{m-1}(p-r_{i}).$ (16)

The reader is referred to [26, 41] for more details.

When building the spline space $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$, all the
pieces of the spline functions can be taken from the “same” ECT-space on the
whole interval $[a,b]$, or from different ECT-spaces considered on the
different intervals. The latter case provides a more general framework which
allows us to locally exploit the full richness of ECT-spaces; however, it also
entails theoretical difficulties, mainly concerning the existence of a basis
possessing all the nice properties of classical polynomial B-splines. This
will be discussed in the next subsection.

Being interested in the application of Tchebycheffian splines in IgA, it is
important to understand the approximation power of the space in (15). It turns
out that, in general, smooth functions can be approximated by Tchebycheffian
splines with the same orders of approximation as in the case of classical
polynomial splines. The reader is referred to [41, Section 9.7] for details
and we state a simplified version of [41, Theorem 9.38] below. Let
$\mathbb{W}^{r,q}([a,b])$ be the standard Sobolev space on $[a,b]$ and let
$\|\cdot\|_{q}$ denote the standard $L^{q}$-norm on $[a,b]$ for $1\leq
q\leq\infty$.

######  Theorem 12.

Let $\mathbb{T}_{p}$ be a $(p+1)$-dimensional ECT-space on the interval
$[a,b]$ and let $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ be a Tchebycheffian
spline space where $\mathbb{T}_{p,i}=\mathbb{T}_{p}$ for $i=1,\ldots,m$. For
any $f\in\mathbb{W}^{p+1,q}([a,b])$ there exists a spline
$Qf\in\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ such that

$\|f-Qf\|_{q}\leq C\,(h_{\Delta})^{p+1}\,\|L_{p+1}f\|_{q},$

where $L_{p+1}$ is defined in (6),

$h_{\Delta}\coloneqq\max_{i=1,\ldots,m}(x_{i}-x_{i-1}),$

and $C$ is a constant independent of $f$ and $\Delta$.

###  2.4 Tchebycheffian B-splines

We now introduce a basis of the space $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$
that possesses all the characterizing properties of classical polynomial
B-splines. In particular, the elements of this basis form a non-negative
partition of unity and have minimal support. These basis functions will be
referred to as Tchebycheffian B-splines (TB-splines).

The existence of TB-splines requires some assumptions on the sequence of ECT-
spaces $\mathbb{T}_{p,i}$, $i=1,\ldots,m$ in (15) which can be expressed in
terms of weight systems; see Section 2.1.

######  Definition 13 (Admissible weights).

For $i=1,\ldots,m$, let $\\{w_{j,i}:j=0,\ldots,p\\}$ be a weight system
generating the ECT-space $\mathbb{T}_{p,i}$ on the interval $[x_{i-1},x_{i}]$.
These weight systems are admissible for the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ in (15), if

$w_{0,i}=1,\quad i=1,\ldots,m,$

and for $i=1,\ldots,m-1$ and $j=0,\ldots,r_{i}$,

$D^{l}_{-}w_{j,i}(x_{i})=D^{l}_{+}w_{j,i+1}(x_{i}),\quad l=0,\ldots,r_{i}-j.$

Given a set of admissible weight systems for the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$, we define the global weight system
$\\{w_{j}:j=0,\ldots,p\\}$ as

$w_{j}(x)\coloneqq w_{j,i}(x),\quad x\in J_{i},\quad i=1,\ldots,m.$ (17)

Similar to the standard polynomial B-splines, Tchebycheffian B-splines can be
defined using a vector of non-decreasing (open) knots,

$\bm{\xi}\coloneqq(\xi_{k})_{k=1}^{n+p+1}\coloneqq\Bigl{(}\underbrace{x_{0},\,\ldots,\,x_{0}}_{p+1\text{
times}},\underbrace{x_{1},\,\ldots,\,x_{1}}_{p-r_{1}\text{
times}},\,\ldots\,,\underbrace{x_{m-1},\,\ldots,\,x_{m-1}}_{p-r_{m-1}\text{
times}},\underbrace{x_{m},\,\ldots,\,x_{m}}_{p+1\text{ times}}\Bigr{)},$ (18)

with $n$ being the dimension of the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$; see (16).

######  Definition 14 (TB-splines).

Assume there exist admissible weight systems as in Definition 13 for the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$. Given the knot vector (18) and the
global weight system (17), the TB-splines $\\{N_{k,p}:k=1,\ldots,n\\}$ are
defined recursively as follows [26]: for $x\in[a,b]$, $q=0,\ldots,p$ and
$k=1,\ldots,n+p-q$,

$N_{k,0}(x)\coloneqq\begin{cases}w_{p}(x),\quad&x\in[\xi_{k},\xi_{k+1}),\\\
0,\quad&\text{otherwise},\end{cases}$

and

$N_{k,q}(x)\coloneqq
w_{p-q}(x)\int_{a}^{x}\Bigg{[}\frac{N_{k,q-1}(y)}{d_{k,q-1}}-\frac{N_{k+1,q-1}(y)}{d_{k+1,q-1}}\Bigg{]}\,{\rm
d}y,\quad q>0,$

where

$d_{j,q-1}\coloneqq\int_{a}^{b}N_{j,q-1}(y)\,{\rm d}y,$

and if $d_{j,q}=0$ then

$\int_{a}^{x}\frac{N_{j,q}(y)}{d_{j,q}}\,dy\coloneqq\begin{cases}1,\quad&x\geq\xi_{j+q+1},\\\
0,\quad&\text{otherwise.}\end{cases}$

TB-splines exhibit a complete structural similarity with polynomial B-splines.
In particular, they enjoy the nice properties collected in the following
proposition; see [26, 36, 38] and references therein.

######  Proposition 15.

Let $\\{N_{k,p}:k=1,\ldots,n\\}$ be the functions in Definition 14. Then,

$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)=\langle N_{k,p}:k=1,\ldots,n\rangle,$

and the following properties hold (see LABEL:{fig_tb_spline_properties}).

  * 1.

Non-negativity: $N_{k,p}(x)>0$ for all $x\in(\xi_{k},\xi_{k+p+1})$;

  * 2.

Local support: $N_{k,p}(x)=0$ for all $x\notin[\xi_{k},\xi_{k+p+1}]$;

  * 3.

Partition of unity: $\sum_{k=1}^{n}N_{k,p}(x)=1$ for all $x\in[a,b]$;

  * 4.

Local linear independence: the set $\\{N_{k,p}\\}_{k=j-p}^{j}$ forms a basis
of $\mathbb{T}_{p,i_{j}+1}$ on $J_{i_{j}+1}$ for all $p+1\leq j\leq n$, where
$i_{j}$ denotes the index such that $\xi_{j}=x_{i_{j}}$;

  * 5.

Interpolation at end-points:

$\displaystyle N_{1,p}(a)$ $\displaystyle=1,\quad N_{k,p}(a)=0,\quad
k=2,\ldots,n;$ $\displaystyle N_{n,p}(b)$ $\displaystyle=1,\quad
N_{k,p}(b)=0,\quad k=1,\ldots,n-1.$

Figure 1: Illustration of TB-spline properties.

(a) $\mathbb{P}_{3}=\langle\,1,x,x^{2},x^{3}\,\rangle$

(b) $\mathbb{P}_{3}^{(10,-10)}=\langle\,1,x,e^{10x},e^{-10x}\,\rangle$

(c) $\mathbb{P}_{3}^{(50,-50)}=\langle\,1,x,e^{50x},e^{-50x}\,\rangle$

(d) $\mathbb{P}_{3}^{(\mathrm{i}\pi.-\mathrm{i}\pi)}=\langle\,1,x,\cos(\pi
x),$ $\sin(\pi x)\,\rangle$

(e)
$\mathbb{P}_{3}^{(\mathrm{i}\frac{3}{2}\pi.-\mathrm{i}\frac{3}{2}\pi)}=\big{\langle}\,1,x,\cos\Big{(}\frac{3\pi}{2}x\Big{)},\sin\Big{(}\frac{3\pi}{2}x\Big{)}\,\big{\rangle}$

(f)
$\mathbb{P}_{6}^{(50,-50,\mathrm{i}\frac{3}{2}\pi.-\mathrm{i}\frac{3}{2}\pi)}=\big{\langle}\,1,x,x^{2},e^{50x},e^{-50x},$
$\cos\Big{(}\frac{3\pi}{2}x\Big{)},\sin\Big{(}\frac{3\pi}{2}x\Big{)}\,\big{\rangle}$

Figure 2: Tchebycheffian Bernstein basis on the interval $[0,1]$ for different
ECT-spaces.

(a) $\mathbb{P}_{3}=\langle\,1,x,x^{2},x^{3}\,\rangle$

(b) $\mathbb{P}_{3}^{(50,-50)}=\langle\,1,x,e^{50x},e^{-50x}\,\rangle$

(c)
$\mathbb{P}_{3}^{(\mathrm{i}\frac{3}{2}\pi.-\mathrm{i}\frac{3}{2}\pi)}=\big{\langle}\,1,x,\cos\Big{(}\frac{3\pi}{2}x\Big{)},$
$\sin\Big{(}\frac{3\pi}{2}x\Big{)}\,\big{\rangle}$

(d)
$\mathbb{P}_{6}^{(50,-50,\mathrm{i}\frac{3}{2}\pi.-\mathrm{i}\frac{3}{2}\pi)}=\big{\langle}\,1,x,x^{2},e^{50x},e^{-50x},\cos\Big{(}\frac{3\pi}{2}x\Big{)},$
$\sin\Big{(}\frac{3\pi}{2}x\Big{)}\,\big{\rangle}$

Figure 3: TB-splines on the partition $\Delta=\\{0,1/3,2/3,1\\}$ with
$\mathbf{r}=\\{2,2\\}$ for different local ECT-spaces.

When $m=1$, so that the partition $\Delta$ has no interior breakpoints, the
space $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ reduces to the ECT-space
$\mathbb{T}_{p}=\mathbb{T}_{p,1}$. The corresponding TB-splines are called
Tchebycheffian Bernstein functions, and play the Tchebycheffian counterpart of
Bernstein polynomials in case of classical algebraic polynomials.
Tchebycheffian Bernstein functions provide an interesting alternative basis to
the generalized power functions of the space $\mathbb{T}_{p}$.

Some examples of the Tchebycheffian Bernstein basis for different ECT-spaces
are depicted in Figure 2, and considering a partition consisting of three
intervals, the corresponding $C^{2}$ TB-splines based on the same ECT-spaces
are given in Figure 3.

Definition 14 allows for the use of different ECT-spaces on different
intervals. However, the requirement of admissible weight systems states that
the various ECT-spaces we are dealing with should be identified by a sequence
of weights which connect smoothly across the different segments. In general,
the construction of such weights is complicated; see [3, 37]. Some comments
about the existence of a TB-spline basis of the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ are in order.

  * 1.

Assume that $\mathbb{T}_{p}$ is a $(p+1)$-dimensional ECT-space on the
interval $[a,b]$ and it contains constants. Moreover, assume that the
derivative space of $\mathbb{T}_{p}$ is a $p$-dimensional ECT-space on $[a,b]$
generated by the weights $\\{w_{1},\ldots,w_{p}\\}$. If all the spaces
$\mathbb{T}_{p,i}$, $i=1,\ldots,m$ in (15) are selected as the restriction of
$\mathbb{T}_{p}$, then from Remark 8 it follows that the weights
$\\{1,w_{1},\ldots,w_{p}\\}$ restricted to the intervals $[x_{i-1},x_{i}]$,
$i=1,\ldots,m$ form admissible weight systems for the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$. In other words, if all the pieces of
the Tchebycheffian spline functions are taken from a single space
$\mathbb{T}_{p}$ which is an ECT-space on $[a,b]$ containing constants, then
the space $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ admits a TB-spline basis
provided that the derivative space of $\mathbb{T}_{p}$ is also an ECT-space on
$[a,b]$; see [26, 33, 41] and references therein.

  * 2.

Assume now that $\mathbb{T}_{p}$ is a $(p+1)$-dimensional ECT-space on each
interval $[x_{i-1},x_{i}]$ separately, but not necessarily on the entire
interval $[a,b]$. Moreover, assume that $\mathbb{T}_{p}$ is the null-space of
a linear differential operator as in (7) and it contains constants. If all the
pieces of the Tchebycheffian spline functions are taken from such space
$\mathbb{T}_{p}$, then the space $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ admits
a TB-spline basis on condition that
$h_{\Delta}\coloneqq\max_{i=1,\ldots,m}(x_{i}-x_{i-1})$ is sufficiently small.
In particular, it admits a TB-spline basis when the derivative space of
$\mathbb{T}_{p}$ is an ECT-space on each interval $[\xi_{k+1},\xi_{k+p}]$,
where $\xi_{k+1}<\xi_{k+p}$ and $1\leq k\leq n$. We refer the reader to [9]
for details (and also [27]). Note that these intervals are closely related to
the supports of the TB-splines; see Proposition 15.

  * 3.

Admissible weight systems can be easily obtained for the interesting subclass
of generalized polynomial splines, i.e., splines where the pieces are taken
from generalized polynomial spaces (such as (12) and (13)), even considering
different spaces on different intervals. This ensures the existence of TB-
splines, called generalized polynomial B-splines (GB-splines) in this setting,
under very mild restrictions. In particular, the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ admits a TB-spline basis when the
$(p-1)$-th order derivative space of $\mathbb{T}_{p,i}$ in (15) is an ECT-
space on $[x_{i-1},x_{i}]$, $i=1,\ldots,m$. We refer the reader to [26,
Section 4] for details.

######  Example 16.

Let $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ be a spline space as in (15) where

$\mathbb{T}_{p,i}=\mathbb{P}_{p}^{(\alpha_{1},\ldots,\alpha_{\ell})}=\Big{\langle}{1,x,\ldots,x^{p-\ell},e^{\alpha_{1}x},\ldots,e^{\alpha_{\ell}x}}\Big{\rangle},\quad
p\geq\ell\geq 0,$

for $i=1,\ldots,m$. From Remark 10 we deduce that
$\mathbb{P}_{p}^{(\alpha_{1},\ldots,\alpha_{\ell})}$ is an ECT-space on the
whole real line, and in particular on the interval $[a,b]$. Moreover, the
derivative space of $\mathbb{P}_{p}^{(\alpha_{1},\ldots,\alpha_{\ell})}$ is
$\mathbb{P}_{p-1}^{(\alpha_{1},\ldots,\alpha_{\ell})}$ (a space without
constants when $p=\ell$), which is also an ECT-space on the whole real line.
We conclude that the corresponding spline space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ admits a TB-spline basis on any
partition of the interval $[a,b]$, without any restrictions.

######  Example 17.

Let $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ be a spline space as in (15) where

$\mathbb{T}_{p,i}=\Big{\langle}{1,x,\ldots,x^{p-\ell-2q},e^{\alpha_{1}x},\ldots,e^{\alpha_{\ell}x},\cos(\beta_{1}x),\sin(\beta_{1}x),\ldots,\cos(\beta_{q}x),\sin(\beta_{q}x)}\Big{\rangle},$

for $i=1,\ldots,m$, i.e., the space considered in (14). The case $q=0$ is
covered in Example 16. Hence, let us here focus on the case $q>0$ and without
loss of generality we assume $\beta_{1}>0,\ldots,\beta_{q}>0$. Since the
pieces of the spline functions belong to the null-space of a linear
differential operator and taking into account Remark 10, we conclude that the
space $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ admits a TB-spline basis on the
interval $[a,b]$ if the knot sequence $\bm{\xi}$ in (18) satisfies

$\xi_{k+p}-\xi_{k+1}<\frac{\pi}{\max\bigl{\\{}\beta_{1},\ldots,\beta_{q}\bigr{\\}}},\quad
k=1,\ldots,n.$

######  Example 18.

Let $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ be a spline space as in (15) where

$\mathbb{T}_{p,i}=\Big{\langle}{1,\cos(\beta x),\sin(\beta
x),\ldots,\cos(q\beta x),\sin(q\beta x)}\Big{\rangle},\quad\beta>0,\quad
p=2q\geq 0,$

for $i=1,\ldots,m$. As discussed in the previous example, such spline space
admits a TB-spline basis on the interval $[a,b]$ if the knot sequence
$\bm{\xi}$ in (18) satisfies $\xi_{k+p}-\xi_{k+1}<{\pi}/{(q\beta)}$,
$k=1,\ldots,n$. Actually, this bound can be improved: the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ admits a TB-spline basis if
$\xi_{k+p}-\xi_{k+1}<{\pi}/{\beta}$, $k=1,\ldots,n$; see [33, 39].

######  Example 19.

Let $\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ be a spline space as in (15) where

$\mathbb{T}_{p,i}=\Big{\langle}{1,x,\ldots,x^{p-2},\cos(\beta x),\sin(\beta
x)}\Big{\rangle},\quad\beta>0,\quad p\geq 2,$

for $i=1,\ldots,m$. As discussed in Example 17, such spline space admits a TB-
spline basis on the interval $[a,b]$ if the knot sequence $\bm{\xi}$ in (18)
satisfies $\xi_{k+p}-\xi_{k+1}<{\pi}/{\beta}$, $k=1,\ldots,n$; see also
Example 6. Actually, it is a generalized polynomial spline space and it can be
shown to have a GB-spline basis if $x_{i}-x_{i-1}<{\pi}/{\beta}$,
$i=1,\ldots,m$; see [26, Section 4].

In the rest of the paper we will only consider Tchebycheffian spline spaces as
in (15) obtained by considering the same ECT-space
$\mathbb{P}_{p}^{\scriptstyle{\mathcal{W}}}$, defined in (14), on each
interval of the partition $\Delta$. In view of the above discussion, this
setting ensures the existence of a TB-spline basis whenever the derivative
space of $\mathbb{P}_{p}^{\scriptstyle{\mathcal{W}}}$ is an ECT-space on
$[a,b]$. If this is not the case, a TB-spline basis can be obtained when a
sufficiently fine partition of the given interval is considered. For
notational convenience, we denote the corresponding TB-splines as

$N_{p,k}^{\scriptstyle{\mathcal{W}}},\quad k=1,\ldots,n,$

in order to stress the selected ECT-space. As numerically illustrated in
Section 4, this simplified Tchebycheffian setting allows already for such a
wide variety of ECT-spaces and related shape parameters to be extremely
profitable in the context of IgA.

######  Remark 20.

The recurrence relation in Definition 14 is the main theoretical tool to
define and analyze TB-splines. However, it is of limited practical interest
because it does not lead to an efficient evaluation strategy. Fortunately,
this deficiency has been recently overcome: efficient routines to construct
and manipulate TB-splines with pieces belonging to ECT-spaces as in Section
2.2 have been provided in [43] by implementing the approach described in [19].
This approach starts from the Tchebycheffian Bernstein functions on each
interval $[x_{i-1},x_{i}]$, $i=1,\ldots,m$, and, by using the smoothness
constraints at all breakpoints, generates an extraction matrix expressing the
TB-spline basis locally in terms of the Tchebycheffian Bernstein functions.
The latter can be computed directly through Hermite interpolation (see Example
22 for an illustration). No weights are required along the process. Although
the evaluation of TB-splines for extreme values of the shape parameters (i.e.,
the roots of the characteristic polynomial (8)), high degrees, high
smoothness, and highly non-uniform partitions still faces numerical challenges
[43], the above mentioned routines provide reliable tools for a large class of
TB-splines of interest in practical applications. A Matlab implementation is
publicly available through the ACM Collected Algorithms (CALGO) library [43].

######  Remark 21.

The TB-spline construction procedure developed in [19, 43] does not require
any knowledge of weights. In case the space
$\mathbb{S}_{p}^{\mathbf{r}}(\Delta)$ does not possess admissible weight
systems, the procedure simply returns a sequence of functions that might not
enjoy the properties listed in Proposition 15. In this perspective, the
computational process can also be used in practice to numerically test the
(non-)existence of a TB-spline basis for a given Tchebycheffian spline space.

######  Example 22.

The space

$\mathbb{P}_{2}^{(\mathrm{i},-\mathrm{i})}=\Big{\langle}{1,\cos(x),\sin(x)}\Big{\rangle}$

is an ECT-space of dimension $3$ on any interval of length less than $2\pi$;
see Example 6. However, the space only admits a Tchebycheffian Bernstein basis
on intervals of smaller length, namely less than $\pi$; see Examples 18 and
19. If such a basis exists, then the corresponding basis elements $B_{j,2}$,
$j=0,1,2$, can be computed by solving the following Hermite interpolation
problems:

$\displaystyle B_{0,2}(a)$ $\displaystyle=1,\quad B_{0,2}(b)=0,\quad
DB_{0,2}(b)=0;$ (19) $\displaystyle B_{1,2}(a)$ $\displaystyle=0,\quad
DB_{1,2}(a)=-DB_{0,2}(a),\quad B_{1,2}(b)=0;$ (20) $\displaystyle B_{2,2}(a)$
$\displaystyle=0,\quad DB_{2,2}(a)=0,\quad B_{2,2}(b)=1.$ (21)

These interpolation problems have a unique solution whenever $b-a<2\pi$, since
$\mathbb{P}_{2}^{(\mathrm{i},-\mathrm{i})}$ is an ECT-space then. However, the
unisolvency of the above interpolation problems does not ensure that the
corresponding functions $B_{j,2}$, $j=0,1,2$ enjoy all the properties listed
in Proposition 15, in particular their non-negativity. This is illustrated in
Figure 4 for the interval $[a,b]=[0,2]$ (so $b-a<\pi$) and the interval
$[a,b]=[0,4]$ (so $\pi<b-a<2\pi$).

Figure 4: The functions satisfying (19)–(21) in the space
$\mathbb{P}_{2}^{(\mathrm{i},-\mathrm{i})}=\Big{\langle}{1,\cos(x),\sin(x)}\Big{\rangle}$.
Left: interval $[0,2]$, Right: interval $[0,4]$.

##  3 Isogeometric analysis with Tchebycheffian B-splines

In this section we outline the isogeometric Galerkin method based on TB-
splines. For the sake of simplicity, we focus on second-order elliptic
differential problems with homogeneous Dirichlet boundary conditions.

###  3.1 Isogeometric Tchebycheffian Galerkin methods

Let ${\mathfrak{L}}$ be a linear second-order elliptic differential operator
on the domain $\Omega\subset\mathbb{R}^{d}$ with Lipschitz boundary
$\partial\Omega$. We consider the differential problem

$\begin{cases}{\mathfrak{L}}u={\rm f},&\text{in\ }\Omega,\\\ u=0,&\text{on\
}\partial\Omega,\end{cases}$ (22)

whose weak form reads as follow:

$\text{find }u\in\mathbb{V}\quad\text{ such that }\quad
a(u,v)=F(v),\quad\forall v\in\mathbb{V},$ (23)

where $\mathbb{V}$ is a suitable function space,
$a:\mathbb{V}\times\mathbb{V}\rightarrow\mathbb{R}$ is a bilinear form
depending on ${\mathfrak{L}}$, and $F:\mathbb{V}\rightarrow\mathbb{R}$ is a
linear form depending on ${\rm f}$.

The Galerkin approach to approximate the solution of (22) is based on the weak
form (23). We select a finite-dimensional approximation space on $\Omega$,

$\mathbb{W}\coloneqq\bigl{\langle}\varphi_{1},\varphi_{2},\ldots,\varphi_{n_{\mathbb{W}}}\bigr{\rangle}\subset\mathbb{V},\quad\dim(\mathbb{W})=n_{\mathbb{W}},$
(24)

and we look for

$u_{\mathbb{W}}\in\mathbb{W}\quad\text{ such that }\quad
a(u_{\mathbb{W}},w)=F(w),\quad\forall w\in\mathbb{W}.$

Taking

$u_{\mathbb{W}}=\sum_{i=1}^{n_{\mathbb{W}}}c_{i}\varphi_{i}$

gives rise to a linear system $A\mathbf{c}=\mathbf{F}$, where the matrix $A$
and the vector $\mathbf{F}$ are defined as

$A_{i,j}\coloneqq a(\varphi_{j},\varphi_{i}),\quad
i,j=1,\ldots,n_{\mathbb{W}},\qquad F_{i}\coloneqq F(\varphi_{i}),\quad
i=1,\ldots,n_{\mathbb{W}}.$ (25)

Different Galerkin methods correspond to different choices of the subspace
$\mathbb{W}$.

Figure 5: Isogeometric Galerkin method. Description of the physical domain
$\Omega$ by means of a global geometry map $\mathbf{G}$.

In the standard formulation of IgA, the physical domain $\Omega$ is
represented by means of an invertible global geometry map. We define the
geometry map $\mathbf{G}$ from the closure of the parametric domain
$\widehat{\Omega}\coloneqq(0,1)^{d}$ to the closure of the physical domain
$\Omega$ as (see Figure 5)

$\mathbf{G}(\hat{\mathbf{x}})\coloneqq\sum_{i=1}^{n_{\mathbf{G}}}{\mathbf{P}}_{i}\hat{\varphi}_{i}(\hat{\mathbf{x}}),\quad{\mathbf{P}}_{i}\in\mathbb{R}^{d},$
(26)

where the basis functions

$\bigl{\\{}\hat{\varphi}_{1},\ldots,\hat{\varphi}_{n_{\mathbf{G}}}\bigr{\\}}$
(27)

have to be selected so as to produce an exact (or at least very accurate)
representation of the geometry (and to allow imposition of homogeneous
boundary conditions). Following the isoparametric approach, the fields of
interest are described with the same basis functions as the geometry in IgA.
The space $\mathbb{W}$ in (24), which incorporates the homogeneous boundary
conditions, is then spanned by the isogeometric functions

$\varphi_{i}(\mathbf{x})\coloneqq\hat{\varphi}_{j_{i}}\circ\mathbf{G}^{-1}(\mathbf{x})=\hat{\varphi}_{j_{i}}(\hat{\mathbf{x}}),\quad
i=1,\ldots,n_{\mathbb{W}},\quad\ \mathbf{x}=\mathbf{G}(\hat{\mathbf{x}}),$
(28)

where $n_{\mathbb{W}}<n_{\mathbf{G}}$ and
$j_{i}:\\{1,\dots,n_{\mathbf{G}}\\}\rightarrow\\{1,\dots,n_{\mathbf{G}}\\}$
are a suitable rearrangement of the indices. In the most common isogeometric
formulation, the functions in (27) are chosen to be tensor-product B-splines
or NURBS. Nevertheless, B-splines/NURBS are not a requisite ingredient in the
isogeometric paradigm and several alternatives have been presented in the
literature.

The properties described in Proposition 15 show that TB-splines are plug-to-
plug compatible with classical (polynomial) B-splines, hence they can be used
as an alternative basis on the parametric domain for building isogeometric
discretization spaces according to (28). In this paper, we focus on TB-splines
from ECT-spaces of the form (14).

We now provide some notation that will be useful in the case studies reported
in Section 4. We deal with two-dimensional problems, where the bivariate TB-
splines on the parametric domain are constructed by taking the tensor product
of univariate TB-splines considering a single ECT-space on the interval
$[0,1]$. Let $\Delta_{1}\times\Delta_{2}$ be a rectangular grid in
$[0,1]^{2}$, let $\mathbf{r}_{\Delta_{1}},\mathbf{r}_{\Delta_{2}}$ be two
smoothness vectors at the corresponding breakpoints, let $p_{1},p_{2}\geq 1$
be two degrees, and let
${\scriptstyle{\mathcal{W}}}_{1},{\scriptstyle{\mathcal{W}}}_{2}$ be two root
vectors as in (10). Then, the tensor-product Tchebycheffian spline space is
given by

$\Big{\langle}N_{p_{1},k}^{{\scriptstyle{\mathcal{W}}}_{1}}\,N_{p_{2},l}^{{\scriptstyle{\mathcal{W}}}_{2}}:k=1,\dots,n_{1},\
l=1,\dots,n_{2}\Big{\rangle},$ (29)

where $n_{1}$ and $n_{2}$ are the dimensions of the two univariate
Tchebycheffian spline spaces. For notational convenience, the basis functions
are also shortly indexed as

$N_{i}\coloneqq
N_{p_{1},k}^{{\scriptstyle{\mathcal{W}}}_{1}}\,N_{p_{2},l}^{{\scriptstyle{\mathcal{W}}}_{2}},\quad
i\coloneqq(l-1)n_{1}+k,\quad k=1,\dots,n_{1},\quad l=1,\ldots,n_{2}.$

A point in (the closure of) the physical domain $\Omega$ is denoted as
$\mathbf{x}=(x,y)$ and a point in (the closure of) the parametric domain
$\widehat{\Omega}$ as $\hat{\mathbf{x}}=(s,t)$. The geometry map in (26) is
given by

$\mathbf{x}=\mathbf{G}(\hat{\mathbf{x}})=\sum_{i=1}^{n_{\mathbf{G}}}\mathbf{P}_{i}N_{i}(\hat{\mathbf{x}}),\quad\mathbf{P}_{i}\in\mathbb{R}^{2},\quad
n_{\mathbf{G}}=n_{1}n_{2},$ (30)

while the basis functions in (28) are

$\phi_{i}(\mathbf{x})\coloneqq
N_{j_{i}}\circ\mathbf{G}^{-1}(\mathbf{x})=N_{j_{i}}(\hat{\mathbf{x}}),$

where, in order to incorporate the homogeneous boundary conditions, the
tensor-product TB-splines

$N_{p_{1},k}^{{\scriptstyle{\mathcal{W}}}_{1}}\,N_{p_{2},l}^{{\scriptstyle{\mathcal{W}}}_{2}},\quad
k=1,n_{1},\quad l=1,\dots,n_{2},\qquad
N_{p_{1},k}^{{\scriptstyle{\mathcal{W}}}_{1}}\,N_{p_{2},l}^{{\scriptstyle{\mathcal{W}}}_{2}},\quad
k=1,\dots,n_{1},\quad l=1,n_{2},$

are excluded according to the last property in Proposition 15. Thus we have
$n_{\mathbb{W}}=(n_{1}-2)(n_{2}-2)$.

###  3.2 Selection of the ECT-spaces

As mentioned above, we are interested in testing the performance of
isogeometric Galerkin methods based on TB-splines from ECT-spaces of the form
(14). It is clear that such a class provides a large variety of combinations
of polynomial, exponential, and trigonometric functions equipped with a wide
spectrum of shape parameters. Of course, it is crucial to properly choose the
structure of the space (including different kinds of functions) and the
related shape parameters to fully exploit the functionality of ECT-spaces.
Particularly, there are two aspects that should be considered in such a
selection to profit from a non-polynomial structure in the reference ECT-
spaces.

  * 1.

The first is the geometric aspect of the problem based on the geometry mapping
at hand. Often, the boundary of the physical domain $\Omega$ is built from
(arcs of) conic sections. Therefore, it can be exactly described in terms of
polynomials, exponential, and trigonometric functions with suitable shape
parameters. In the pure isogeometric philosophy, the geometry of the physical
domain is exact from the beginning of the procedure and does not change along
the process when possible refinements are required. In this perspective, it is
then natural to select ECT-spaces (14) that allow for an exact representation
of the geometry. In its classical formulation, IgA is based on B-splines or
NURBS. Even though (arcs of) conic sections can be exactly represented in
terms of NURBS, it is not an arc-length parameterization and smooth
representations of (closed) conic sections require high degrees. As an
example, a $C^{r}$ NURBS parameterization of a circle requires at least degree
$2(r+1)$; see [2]. Moreover, NURBS behave poorly when dealing with
differentiation and integration. On the other hand, spaces of the form (14)
allow for natural smooth parameterizations of conic sections and their
“structure” does not change when dealing with differentiation and integration,
in a complete similarity with polynomial spaces; see also Remark 11.

  * 2.

The second, of interest in the numerical treatment of PDEs, is the fact that
functions in ECT-spaces of the form (14) are fundamental solutions of
important ordinary differential operators, meaning that the shape parameters
are immediately connected to the associated differential problems. Although
the extension to the multivariate case is not straightforward, the use of TB-
splines with pieces drawn from suitable ECT-spaces can still be beneficial. As
an example, they can lead to a significant reduction, or even a complete
removal, of extraneous oscillations in advection-dominated problems, without
the need of stabilization techniques. Once again, a proper selection of the
shape parameters is crucial but can be driven by the problem setting.

In the next section we present some case studies that show how TB-splines can
outperform classical polynomial B-splines in isogeometric Galerkin
discretizations of Poisson and advection-diffusion problems. As discussed
above, the selection of the underlying ECT-spaces has to be done according to
a problem-driven strategy; this will be detailed in each of the case studies.

##  4 Case studies

In this section we discuss the performance of the isogeometric Galerkin method
based on TB-splines through different case studies. In all the numerical
experiments we use the following setup.

Bivariate tensor-product Tchebycheffian spline spaces of the form (29) are
taken as discretization spaces in the isogeometric Galerkin method. The
corresponding ECT-spaces in the two directions will be specified in each
experiment by writing

$\mathbb{P}_{p_{1}}^{{\scriptstyle{\mathcal{W}}}_{1}}\otimes\mathbb{P}_{p_{2}}^{{\scriptstyle{\mathcal{W}}}_{2}},$

where $p_{1},p_{2}$ are the degrees and
${\scriptstyle{\mathcal{W}}}_{1},{\scriptstyle{\mathcal{W}}}_{2}$ are the root
vectors as in (10). We only consider uniform partitions
$\Delta_{1}=\Delta_{2}$ with breakpoints

$x_{i}=\frac{i}{m},\quad i=0,\dots,m,$

and denote by $h$ the (uniform) distance between two consecutive breakpoints,
i.e., $h=\frac{1}{m}$. We assume maximal smoothness at the breakpoints, i.e.,

$\mathbf{r}_{\Delta_{1}}=\\{p_{1}-1,\dots,p_{1}-1\\},\quad\mathbf{r}_{\Delta_{2}}=\\{p_{2}-1,\dots,p_{2}-1\\}.$

This amounts to univariate Tchebycheffian spline spaces of dimension
$n_{1}=m+p_{1}$ and $n_{2}=m+p_{2}$, respectively.

We focus on the general advection-diffusion problem of the form

$\begin{split}\begin{cases}-\nabla\cdot({\kappa}\nabla
u)+\textbf{a}\cdot\nabla u={\rm f},\quad&\text{in}\
\Omega=\mathbf{G}\big{(}(0,1)^{2}\big{)},\\\ u={\rm g},\quad&\text{in}\
\partial\Omega,\end{cases}\end{split}$ (31)

where $\mathbf{a}$ is the advection flow velocity, $\kappa$ the diffusivity,
${\rm f}$ the prescribed source function, and ${\rm g}$ the Dirichlet boundary
data. The physical domain $\Omega$ is described by the geometry map
$\mathbf{G}$. The dominance of advection over diffusion is defined by the
global Péclet number

$\bm{Pe}_{g}\coloneqq\frac{\|\bm{a}\|}{\kappa}.$ (32)

For simplicity, we fix the diffusion coefficient $\kappa=1$ in all the case
studies.

As usual, homogeneous boundary conditions are satisfied pointwise exactly. In
the non-homogeneous case, the boundary values are represented in the
underlying Tchebycheffian spline space by suitable approximation strategies
(e.g., least squares approximation) and subsequently the reduction to the
homogeneous case is considered, so dealing again with a special instance of
the problem (22).

For an accurate numerical approximation of the integrals needed in the
construction of the matrix and the vector in (25), we use elementwise Gaussian
quadrature rules. The selection of quadrature points for discretizations based
on TB-splines is not trivial and quadrature rules of higher order compared to
the algebraic polynomial splines are often required. In the presented case
studies we use $3p$ quadrature points in each element for TB-splines only
involving trigonometric functions (see Sections 4.1 and 4.2), and $5p$
quadrature points in each element when also exponential functions with large
shape parameters are considered (see Sections 4.3, 4.4 and 4.5).

###  4.1 Case study 1: Poisson problem on a domain bounded by circular arcs

Figure 6: Case study 1. Left: Geometry map presented in (33) as a combination
of curves with phases $\beta$ and $2\beta$. Right: Plot of the approximate
solution obtained by using tensor-product TB-splines identified by
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{4}$}}}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{4},-\mathrm{i}\frac{\pi}{4},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$
with $m=8$ and dof $=144$.

In this case study we address a Poisson problem by considering (31) with
$\kappa=1$, $\textbf{a}=\textbf{0}$, and ${\rm f}$ obtained from the exact
solution

$u(x,y)=5((x+4)^{2}+(y+3)^{2})+xy.$

The parametric domain $\widehat{\Omega}=(0,1)^{2}$ is mapped to the physical
domain $\Omega$ through the geometry map given by

$\begin{pmatrix}x\\\ y\end{pmatrix}=\textbf{G}\begin{pmatrix}s\\\
t\end{pmatrix}=(1-s)r\,\begin{pmatrix}\cos(2\beta t)\\\ \sin(2\beta
t)\end{pmatrix}+\ s\,\begin{pmatrix}-C+R\cos(\beta t)\\\ R\sin(\beta
t)\end{pmatrix},$ (33)

with

$r=1,\quad\beta=\frac{\pi}{4},\quad
C=\frac{2r\cos(\beta)}{1-\cos(\beta)},\quad R=C+2r,$

which leads to the domain depicted in Figure 6, left. It is evident that the
geometry can be exactly represented in the form (30) by considering in the
parametric $t$-direction Tchebycheffian splines identified by an ECT-space
that contains the trigonometric functions with phases $\beta$ and $2\beta$.
Hence we consider in the $t$-direction the space of degree $p_{2}=4$,

$\mathbb{P}_{4}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\mathrm{i}2\beta,-\mathrm{i}2\beta)}=\Biggl{\langle}\,1,\cos\biggl{(}\frac{\pi}{4}t\biggr{)},\sin\biggl{(}\frac{\pi}{4}t\biggr{)},\cos\biggl{(}\frac{\pi}{2}t\biggr{)},\sin\biggl{(}\frac{\pi}{2}t\biggr{)}\,\Biggr{\rangle}.$

The Tchebycheffian Bernstein basis of this ECT-space is illustrated in Figure
7. From Example 18 we know that the corresponding Tchebycheffian spline space
admits a TB-spline basis on any partition of the interval $[0,1]$ since
$\beta<\pi$. In the $s$-direction, we simply consider algebraic polynomials of
degree $p_{1}=4$, i.e.,

$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{4}$}}}=\Big{\langle}{1,s,s^{2},s^{3},s^{4}}\Big{\rangle}.$

This gives the tensor-product space

$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{4}$}}}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{4},-\mathrm{i}\frac{\pi}{4},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}.$

To obtain the control coefficients $\mathbf{P}_{i}$ of the geometry in the
form (30) we simply do interpolation on the grid obtained by the tensor
product of the Greville points of polynomial B-splines of the same degree $4$.
We remark that this geometry cannot be exactly represented by using
generalized polynomial B-splines (see [28]) with pieces in ECT-spaces of the
form (13).

Figure 7: Case study 1. Tchebycheffian Bernstein basis of
$\mathbb{P}_{4}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\mathrm{i}2\beta,-\mathrm{i}2\beta)}$
on the interval $[0,1]$.

We approximate the $L^{\infty}$ norm of the error by sampling the approximate
and exact solutions on a uniform grid in the parametric domain consisting of
$501$ points along each direction. Section 4.1 shows the $L^{\infty}$ error
obtained for different levels of refinements, starting with a very coarse mesh
of only one interval in each parametric direction, comparing the TB-splines
against the classical B-splines for spaces of the same dimension. We see an
improvement of about a factor $10$ that we achieve just by exactly
representing the geometry using TB-splines. The contour plot of the
approximate solution partitioning the domain in $8\times 8$ elements is
illustrated in Figure 6, right.

$m$ |  dof |  $\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{4}$}}}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{4},-\mathrm{i}\frac{\pi}{4},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$ |  $\mathbb{P}_{4}\otimes\mathbb{P}_{4}$
---|---|---|---
1 |  25 |  7.9294e-04 |  5.0629e-03
2 |  36 |  3.4854e-04 |  5.9308e-04
4 |  64 |  2.9655e-05 |  8.4760e-05
8 |  144 |  8.3556e-07 |  2.7812e-06
16 |  400 |  2.6336e-08 |  9.2305e-08
32 |  1296 |  8.4881e-10 |  2.9782e-09

Table 1: Case study 1. Comparison of the $L^{\infty}$ error of tensor-product
TB-splines and B-splines of degrees $p_{1}=p_{2}=4$ for different numbers of
intervals $m$ in each direction.

###  4.2 Case study 2: Poisson problem on a symmetric domain bounded by
circular arcs

Here we focus again on a Poisson problem by considering (31) with $\kappa=1$,
$\textbf{a}=\textbf{0}$, and ${\rm f}$ obtained from the exact solution u(x,y)
= 4((x+4)^2+y^2) + x(y-4), to exhibit the evolution of error under uniform
refinement and verify the consistency of the TB-splines for optimal
convergence. The geometry map from the parametric domain
$\widehat{\Omega}=(0,1)^{2}$ to the physical domain $\Omega$ is given by

$\begin{pmatrix}x\\\ y\end{pmatrix}=\textbf{G}\begin{pmatrix}s\\\
t\end{pmatrix}=(1-s)r\,\begin{pmatrix}\cos(\beta t)\\\ \sin(\beta
t)\end{pmatrix}+\ sR\,\begin{pmatrix}C+\cos(2\beta t-\gamma)\\\ C+\sin(2\beta
t-\gamma)\end{pmatrix},$ (34)

with r = 1, R = 2, C = 1/2, γ= π4, β= π2, which leads to the domain depicted
in Figure 8, left. Similar to the previous case study, the geometry consists
of two arcs with phases $\beta$ and $2\beta$, hence a suitable selection of
the ECT-space in the $t$-direction is

$\mathbb{P}_{p_{2}}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\dots,\mathrm{i}q\beta,-\mathrm{i}q\beta)}=\Big{\langle}{1,\cos(\beta
t),\sin(\beta t),\dots,\cos(q\,\beta t),\sin(q\,\beta t)}\Big{\rangle},\quad
p_{2}=2q\geq 4.$ (35)

The Tchebycheffian Bernstein basis of this ECT-space is illustrated in Figure
9 for $p_{2}=6$. From Example 18 we know that the corresponding Tchebycheffian
spline space admits a TB-spline basis on any partition of the interval $[0,1]$
since $\beta<\pi$. In the $s$-direction, we simply consider algebraic
polynomials of the same degree as in the $t$-direction. This gives the tensor-
product space P_p_1 ⊗P_p_2^(iπ2,-iπ2,…, iqπ2, -iqπ2), p_1=p_2=p=2q. Again,
this geometry cannot be exactly represented by using generalized polynomial
B-splines (see [28]) with pieces in ECT-spaces of the form (13).

Figure 8: Case study 2. Left: Geometry map presented in (34) as a combination
of curves with phases $\beta$ and $2\beta$. Right: Plot of the approximate
solution obtained by using tensor-product TB-splines identified by
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{6}$}}}\otimes\mathbb{P}_{6}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\mathrm{i}2\beta,-\mathrm{i}2\beta,\mathrm{i}{3\beta},-\mathrm{i}{3\beta})}$,
$\beta=\frac{\pi}{2}$ with $m=8$ and dof $=196$. Figure 9: Case study 2.
Tchebycheffian Bernstein basis of
$\mathbb{P}_{6}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\mathrm{i}2\beta,-\mathrm{i}2\beta,\mathrm{i}3\beta,-\mathrm{i}3\beta)}$,
$\beta=\frac{\pi}{2}$, on the interval $[0,1]$. Figure 10: Case study 2.
Convergence plot of the $L^{\infty}$ error of tensor-product TB-splines
identified by
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{p}$}}}\otimes\mathbb{P}_{p}^{(\mathrm{i}\beta,-\mathrm{i}\beta,\mathrm{i}2\beta,-\mathrm{i}2\beta,\dots,\mathrm{i}{q\beta},-\mathrm{i}{q\beta})}$,
$\beta=\frac{\pi}{2}$, for $p=2q$, $q=\\{1,2,3\\}$, and different resolution
of the mesh $h=1,\frac{1}{2},\frac{1}{2^{2}},\dots,\frac{1}{2^{5}}$.

Figure 10 illustrates the convergence of the error in $L^{\infty}$ norm, which
has been approximated by sampling the approximate and exact solutions on a
uniform grid in the parametric domain consisting of $501$ points along each
direction. The plot clearly shows that the approximate solution obtained by
TB-splines exhibits the optimal convergence of order $p+1$ as classical
polynomial B-splines in agreement with the theoretical expectations; see
Theorem 12. This establishes the reliability of the Tchebycheffian spline
spaces from the accuracy point of view. The contour plot of the approximate
solution with $p=6$ on a mesh consisting of $8\times 8$ elements is presented
in Figure 8, right.

###  4.3 Case study 3: Advection-diffusion problem with radial advection flow

In the experiments presented so far, the selection of the ECT-spaces has been
driven by their geometrical features. Now we will also take into account some
analytical features of the problem. Let us consider the advection-diffusion
problem (31) with

$\kappa=1,\quad\mathbf{a}\,(x,y)\coloneqq\mathfrak{a}\begin{pmatrix}\frac{x}{\sqrt{x^{2}+y^{2}}}\\\
\frac{y}{\sqrt{x^{2}+y^{2}}},\end{pmatrix},\quad\mathfrak{a}=100,\quad{\rm
f}(x,y)=1,$ (36)

and homogeneous boundary conditions, on the same geometry used in the previous
case study given by the geometry map (34).

It is well known that for problems with high Péclet number $(\bm{Pe}_{g}\gg
1)$ approximate solutions belonging to piecewise polynomial spaces tend to
exhibit spurious oscillations unless the discretization is fine enough. The
above mentioned issue can be efficiently addressed by turning to
Tchebycheffian splines built from exponential functions with suitable
parameters. Their use for solving advection-dominated advection-diffusion
problems reduces, and sometimes even eliminates, the spurious oscillations
while still providing an accurate localization of the sharp layer(s) without
demanding extremely fine meshes and so keeping the number of degrees of
freedom (dof) reasonably low.

Figure 11: Case study 3. Illustration of the selection strategy of the shape
parameters $\\{\alpha_{i}:i=1,\ldots,\ell\\}$ of the exponential functions to
accommodate for different distances traveled by the advection flow for
$\ell=3$ in (37). The extremes of the chosen shape parameters
$\alpha_{1}=\alpha_{min}=\mathfrak{a}\,(\sqrt{2}R-r)$ and
$\alpha_{3}=\alpha_{max}=\mathfrak{a}\,(2R-r)$ are evidently related to the
maximum and minimum radial distances of the domain, where $\mathfrak{a}$ is
the modulus of the advection coefficient.

For an advection-dominated system with advection flow in the radial direction
as in (36), we can exploit exponential functions in the parametric
$s$-direction with suitable shape parameters to represent the sharp layer
generated near the homogeneous boundary. To construct effective Tchebycheffian
splines with exponential functions in the $s$-direction, we consider the ECT-
space

$\mathbb{P}_{p_{1}}^{(\alpha_{1},\alpha_{2},\ldots,\alpha_{\ell})}=\Big{\langle}{1,s,\ldots,s^{p_{1}-\ell},e^{\alpha_{1}s},\dots
e^{\alpha_{\ell}s}}\Big{\rangle},\quad\quad p_{1}\geq\ell\geq 1,$

where for $i=1,\ldots,\ell,$

$\alpha_{i}\coloneqq\begin{cases}\frac{1}{2}(\alpha_{max}+\alpha_{min}),&\ell=1,\\\\[2.84544pt]
\alpha_{min}+\frac{i-1}{\ell-1}(\alpha_{max}-\alpha_{min}),&\ell\geq
2,\end{cases}$ (37)

and α_min≔a (2R - r) , α_max≔a (2R - r). The governing factors in the
selection of the shape parameters here are the global Péclet number (32),
which is $\bm{Pe}_{g}=\mathfrak{a}$, and the range of distances traveled by
the advection flow in the radial direction. Figure 11 illustrates this
approach for the selection of the shape parameters by considering $\ell=3$
different radial distances traveled by the advection flow. The Tchebycheffian
Bernstein basis of this ECT-space is illustrated in Figure 12 for $p_{1}=4$.
From Example 16 we know that the corresponding Tchebycheffian spline space
admits a TB-spline basis on any partition of the interval $[0,1]$.

Figure 12: Case study 3. Tchebycheffian Bernstein basis of
$\mathbb{P}_{4}^{(\alpha_{1},\alpha_{2},\alpha_{3})}$ on the interval $[0,1]$,
with a magnified view of the functions on the interval $[0.95,1]$.

As discussed in Section 4.2, to exactly reproduce the geometry, we consider in
the $t$-direction the ECT-space (35) of degree $p_{2}=2q\geq 4$, which allows
for a Tchebycheffian spline space equipped with a TB-spline basis on the
interval $[0,1]$. The resulting tensor-product space is P_p_1^(α_1,α_2,…,α_ℓ)
⊗P_p_2^(iπ2,-iπ2,…, iqπ2, -iqπ2), p_2=2q, with both spaces of the form (14).

(a) TB-spline space:
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{4}$}}}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2},\mathrm{i}\pi,-\mathrm{i}\pi)}$

(b) TB-spline space:
$\mathbb{P}_{4}^{(\alpha_{1})}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2},\mathrm{i}\pi,-\mathrm{i}\pi)}$

(c) TB-spline space:
$\mathbb{P}_{4}^{(\alpha_{1},\alpha_{2})}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2},\mathrm{i}\pi,-\mathrm{i}\pi)}$

(d) TB-spline space:
$\mathbb{P}_{4}^{(\alpha_{1},\alpha_{2},\alpha_{3})}\otimes\mathbb{P}_{4}^{(\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2},\mathrm{i}\pi,-\mathrm{i}\pi)}$

Figure 13: Case study 3. Plots of the approximate solution to the advection-
diffusion problem with global Péclet number $\bm{Pe}_{g}=100$, obtained by
using tensor-product TB-splines with $p_{1}=p_{2}=4$, $m=8$, dof = $144$.
Increasing the number of exponential functions with shape parameters
$\\{\alpha_{i}:i=1,\ldots,\ell\\}$ according to (37) in the TB-spline space
results into a reduction of oscillations in the solution.

In Figure 13 we present some results obtained by using tensor-product TB-
splines of degrees $p_{1}=p_{2}=4$ on a mesh consisting of $8\times 8$
elements. To exhibit the effect of the presence of exponential functions in
the TB-spline space for addressing the advection-dominant problem (31) with
(36), we vary the number of shape parameters in the selection from (37). The
outcome of using only algebraic polynomial splines (so without any exponential
functions) in the $s$-direction is shown in Figure 13(a). Even though the
geometry is described exactly with the trigonometric functions present in the
$t$-direction, the absence of the exponential functions results in spurious
oscillations in the solution. However, taking one exponential function in the
TB-spline space somewhat reduces the oscillations as shown in Figure 13(b),
and taking two exponential functions in the space almost removes all the over-
and undershoots as shown in Figure 13(c). Following the trend, with three
exponential functions in the TB-spline space, the oscillations are completely
eliminated from the solution (at least visually) as illustrated in Figure
13(d).

###  4.4 Case study 4: Advection-diffusion problem with tangential advection
flow

The ECT-spaces used in the previous case studies do not exploit the full
functionality of mixing exponential and trigonometric functions with
polynomials, as in (14). Hence, now we consider a problem with analytical and
geometrical features in the same parametric direction. Let us consider the
advection-diffusion problem (31) with

$\kappa=1,\quad\mathbf{a}\,(x,y)\coloneqq\mathfrak{a}\begin{pmatrix}\frac{-y}{\sqrt{x^{2}+y^{2}}}\\\
\frac{x}{\sqrt{x^{2}+y^{2}}},\end{pmatrix},\quad\mathfrak{a}=100,\quad{\rm
f}(x,y)=1,$ (38)

and homogeneous boundary conditions, on a classical quarter of an annulus
given by the geometry map

$\begin{pmatrix}x\\\ y\end{pmatrix}=\textbf{G}\begin{pmatrix}s\\\
t\end{pmatrix}=((1-s)r+sR)\begin{pmatrix}\cos(\frac{\pi}{2}t)\\\
\sin(\frac{\pi}{2}t)\end{pmatrix},\quad r=1,\quad R=2.$ (39)

which leads to the domain depicted in Figure 5, right.

Figure 14: Case study 4. Illustration of the selection strategy of the shape
parameters $\\{\alpha_{i}:i=1,\ldots,\ell\\}$ of the exponential functions to
accommodate for different distances traveled by the advection flow for
$\ell=3$ in (41). The extremes of the chosen shape parameters
$\alpha_{1}=\alpha_{min}=\mathfrak{a}\,\frac{\pi}{2}r$ and
$\alpha_{3}=\alpha_{max}=\mathfrak{a}\,\frac{\pi}{2}R$ are evidently related
to the two extreme arc lengths of the domain, where $\mathfrak{a}$ is the
modulus of the advection coefficient.

To build suitable TB-spline spaces for the above problem, we need to consider
both the geometry map in (39) and the advection coefficient $\bm{a}$ in (38).
The geometry map requires trigonometric functions in the $t$-direction.
Moreover, since we have a tangential flow along the $t$-direction, see Figure
14, it is natural to add exponential functions in that direction as well to
accurately represent the sharp layer. The ECT-space in the $t$-direction can
be chosen as

$\mathbb{P}_{p_{2}}^{(\alpha_{1},\alpha_{2},\ldots,\alpha_{\ell},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}=\Biggl{\langle}\,1,t,\ldots,t^{p_{2}-2-\ell},e^{\alpha_{1}t},\ldots
e^{\alpha_{\ell}t},\cos\biggl{(}\frac{\pi}{2}t\biggr{)},\sin\biggl{(}\frac{\pi}{2}t\biggr{)}\,\Biggr{\rangle},\quad
p_{2}\geq\ell+2\geq 3,$ (40)

where, similar to (37), for $i=1,\dots,\ell$,

$\alpha_{i}\coloneqq\begin{cases}\frac{1}{2}(\alpha_{max}+\alpha_{min}),&\ell=1,\\\\[2.84544pt]
\alpha_{min}+\frac{i-1}{\ell-1}(\alpha_{max}-\alpha_{min}),&\ell\geq
2,\end{cases}$ (41)

and α_min≔a π2r, α_max≔a π2R. The governing factors in the selection of the
shape parameters here are again the global Péclet number (32), which is
$\bm{Pe}_{g}=\mathfrak{a}$, and the range of distances traveled by the
advection flow along circular arcs. Figure 14 illustrates this approach for
the selection of the shape parameters by considering $\ell=3$ different
distances along arcs traveled by the advection flow. The Tchebycheffian
Bernstein basis of this ECT-space is illustrated in Figure 15 for $p_{2}=6$.
In the $s$-direction, we simply consider algebraic polynomials of degree
$p_{1}=2$. The resulting tensor-product space is P_p_1
⊗P_p_2^(α_1,α_2,…,α_ℓ,iπ2,-iπ2).

Figure 15: Case study 4. Tchebycheffian Bernstein basis of
$\mathbb{P}_{6}^{(\alpha_{1},\alpha_{2},\alpha_{3},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$
on the interval $[0,1]$, with a magnified view of the functions on the
interval $[0.95,1]$.

(a) GB-spline space:
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{2}$}}}\otimes\mathbb{P}_{6}^{(\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$

(b) TB-spline space:
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{2}$}}}\otimes\mathbb{P}_{6}^{(\alpha_{1},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$

(c) TB-spline space:
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{2}$}}}\otimes\mathbb{P}_{6}^{(\alpha_{1},\alpha_{2},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$

(d) TB-spline space:
$\settoboxheight{}{}{\mathbb{P}}{\raisebox{0.0pt}{\leavevmode\hbox{\set@color${}_{2}$}}}\otimes\mathbb{P}_{6}^{(\alpha_{1},\alpha_{2},\alpha_{3},\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}$

Figure 16: Case study 4. Plots of the approximate solution to the advection-
diffusion problem with global Péclet number $\bm{Pe}_{g}=100$, obtained by
using tensor-product TB-splines with $p_{1}=2$, $p_{2}=6$, $m=6$, dof = $96$.
Increasing the number of exponential functions with different shape parameters
$\\{\alpha_{i}:i=1,\ldots,\ell\\}$ according to (41) in the TB-spline space
results into a reduction of oscillations in the solution.

For the sake of comparison, we also consider in the $t$-direction the
generalized polynomial space obtained by adding to the polynomial space the
pair of trigonometric functions necessary for an exact representation of the
geometry, see (13), i.e.,

$\mathbb{P}_{p_{2}}^{(\mathrm{i}\frac{\pi}{2},-\mathrm{i}\frac{\pi}{2})}=\Biggl{\langle}\,1,t,\dots,t^{{p_{2}}-2},\cos\biggl{(}\frac{\pi}{2}t\biggr{)},\sin\biggl{(}\frac{\pi}{2}t\biggr{)}\,\Biggr{\rangle},$
(42)

so that the resulting tensor-product space used for the comparison is P_p_1
⊗P_p_2^(iπ2,-iπ2).

Note that both the ECT-spaces in (40) and (42) are of the form (14). From
Example 17 we know that the corresponding Tchebycheffian spline spaces admit a
TB-spline basis on any partition of the interval $[0,1]$.

In Figure 16 we present some results obtained by using tensor-product TB-
splines of degrees $p_{1}=2$ and $p_{2}=6$ on a mesh consisting of $6\times 6$
elements. To exhibit the effect of the presence of exponential functions in
the TB-spline space for addressing the advection-dominant problem (31) with
(38), we vary the number of shape parameters in the selection from (41). We
can see in Figure 16(a) that the GB-spline space with only trigonometric
functions leads to spurious oscillations in the solution. However, taking just
one exponential function in the TB-spline space already substantially reduces
the oscillations as shown in Figure 16(b), and with two exponential functions
all the over- and undershoots almost disappear as shown in Figure 16(c). The
trend continues to follow and increasing the number of exponential functions
to $3$ ultimately gives us a solution with no oscillations (at least visually)
as illustrated in Figure 16(d).

###  4.5 Case study 5: Advection-diffusion problem with internal sharp layer
on a square

In the previous two case studies the advection flow travels exactly along one
of the parametric directions. Here we consider a classical benchmark problem
with advection flow not parallel to any parametric direction [20, 29, 44].

Figure 17: Case study 5. The domain with the Dirichlet boundary conditions,
the sharp internal layer (red) generated along the advection flow at angle
$\theta$ and the boundary layers (red).

We want to solve the advection-diffusion problem (31) on a square, hence the
geometry map is just the identity, with κ= 1, a = a (cos(θ), sin(θ))^T, θ=π4,
a= 10^4, f(x,y) = 0, and the discontinuous Dirichlet boundary conditions as
shown in Figure 17. The solution exhibits sharp boundary and internal layers
depicted in red in the same figure. In particular, the sudden jump in the
boundary conditions at the point $(0,0.2)$ generates an inner sharp layer
aligned with the advection flow direction identified by
$(\cos(\theta),\sin(\theta))$.

As already mentioned in Section 4.3, since the problem has a large global
Péclet number, namely $\bm{Pe}_{g}=\mathfrak{a}=10^{4}$, approximate solutions
belonging to piecewise polynomial spaces exhibit spurious oscillations until
the discretization is fine enough to resolve the sharp layers featured by the
exact solution. To overcome this issue, it is common to rely on stabilization
methods such as the SUPG and GLS method; see [10, 21]. Typically,
stabilization methods remove spurious oscillations but at the same time
“smooth out” the layers featured by the exact solution. Their efficiency
deeply depends on a careful choice of some parameters appearing in the various
stabilization methods.

Tchebycheffian splines can offer a flexible alternative for the treatment of
the above advection-dominated problem without the need for stabilization. To
catch the sharp layers, it is natural to consider Tchebycheffian spline spaces
that contain exponential functions with suitable parameters. Since the
advection flow velocity is constant here, we can simply take its two
components as parameters in the two (parameter) directions:

$\mathbb{P}_{p_{1}}^{(\mathfrak{a}\cos(\theta))}=\Big{\langle}{1,x,\dots,x^{p_{1}-1},e^{\mathfrak{a}\cos(\theta)x}}\Big{\rangle},\quad\mathbb{P}_{p_{2}}^{(\mathfrak{a}\sin(\theta))}=\Big{\langle}{1,y,\dots,y^{p_{2}-1},e^{\mathfrak{a}\sin(\theta)y}}\Big{\rangle}.$
(43)

Both these spaces are of the form (14) and, since they just contain polynomial
and exponential functions, from Example 16 we know that the corresponding
Tchebycheffian spline spaces admit a TB-spline basis on any partition of the
interval $[0,1]$. Note that, since the physical domain is simply a square, the
parameter selection in (43) is not governed by the geometry.

It is important to remark that a proper treatment of the boundary conditions
is imperative due to its effect on the approximate solution. In order to avoid
oscillations along the boundary, instead of using a least squares
approximation as in the previous case studies, we have considered quasi-
interpolation. More precisely, in order to exploit the shape-preserving
property of the (T)B-spline representation, following the approach in [20], we
approximate the boundary function through the Schoenberg operator (see [25,
40]) along each edge. In other words, we consider the linear combination of
the boundary TB-splines whose coefficients are obtained by evaluating the
boundary function at the corresponding Greville abscissae (under the
assumptions that linear polynomials belong to the ECT-space of interest).

(a) TB-spline space:
$\mathbb{P}_{6}^{(\mathfrak{a}\cos(\theta))}\otimes\mathbb{P}_{6}^{(\mathfrak{a}\sin(\theta))}$
with $m=50$, dof $=3136$.

(b) B-spline space: $\mathbb{P}_{6}\otimes\mathbb{P}_{6}$ with SUPG
stabilization and $m=50$, dof $=3136$.

(c) B-spline space: $\mathbb{P}_{6}\otimes\mathbb{P}_{6}$ with $m=1500$, dof
$=2268036$.

Figure 18: Case study 5. Plots of the approximate solution to the advection-
diffusion problem with advection skew at $\theta=45^{\circ}$ and global Péclet
number $\bm{Pe}_{g}=10^{4}$, obtained by using tensor-product spline spaces
with $p_{1}=p_{2}=6$; the number of elements $m$ in each direction and the dof
are mentioned in the subcaptions.

In Figure 18(a) we illustrate the approximate solution obtained by using
tensor-product TB-splines, identified by (43), of degrees $p_{1}=p_{2}=6$ on a
mesh consisting of $50\times 50$ elements. If we use classical tensor-product
B-splines of the same degrees on a mesh of only $50\times 50$ elements, then
the result will exhibit huge spurious oscillations near the layers. Hence, it
makes sense to compare the stabilized solution of the B-splines against the
TB-splines (without stabilization). Figure 18(b) shows the result obtained by
using tensor-product B-splines with SUPG stabilization, where the
stabilization constant is given by

$\tau=\frac{h_{\mathbf{a}}}{2\|\mathbf{a}\|},\quad\text{with }\quad
h_{\mathbf{a}}=\frac{h}{\max\bigl{\\{}\cos(\theta),\sin(\theta)\bigr{\\}}}.$

Note that $h_{\mathbf{a}}$ is the element length in the direction of the
advection flow; see [20] and [15, Section 9.3.2]. As a reference solution we
are considering the one obtained by classical tensor-product B-splines of the
same degrees on a very fine mesh of $1500\times 1500$ elements, as shown in
Figure 18(c).

TB-spline space |  SUPG stabilization |  $m$ |  dof |  max |  min
---|---|---|---|---|---
$\mathbb{P}_{6}^{(\mathfrak{a}\cos(\theta))}\otimes\mathbb{P}_{6}^{(\mathfrak{a}\sin(\theta))}$ |  No |  50 |  3136 |  1.0016 |  $-1.6385\times 10^{-3}$
$\mathbb{P}_{6}\otimes\mathbb{P}_{6}$ |  Yes |  50 |  3136 |  1.0009 |  $-9.4671\times 10^{-4}$
$\mathbb{P}_{6}\otimes\mathbb{P}_{6}$ |  No |  1500 |  2268036 |  1.0022 |  $-1.8761\times 10^{-6}$
Table 2: Case study 5. Comparison of maximum and minimum of the TB-spline
solution with exponential functions in the space against the polynomial
B-spline solution with SUPG stabilization on the same mesh and the polynomial
B-spline solution on a really high number of intervals with $p_{1}=p_{2}=6$
for $\mathfrak{a}=10^{4}$.

Table 2 collects the values of maximal over- and undershoot in the
neighbourhood of the layers evaluated on a uniform grid of $3001$ points along
each direction for the three considered spline spaces, together with the used
number of elements and the total number of dof.

When comparing the TB-spline solution with the stabilized B-spline solution,
we observe that both approaches undoubtedly eliminate the spurious
oscillations. However, stabilization produces a too smooth solution and the
localization of the internal and boundary layers is less accurate than in the
case of the TB-spline solution. On the other hand, the unstabilized B-spline
solution on a very fine mesh is able to localize accurately the layers but at
a very high cost (i.e., a large number of dof). In conclusion, TB-splines on
coarse meshes are able to predict the layers precisely, resulting in very few
required dof while still giving stable results and eliminating the need of
stabilization.

##  5 Conclusion

Tchebycheffian splines are smooth piecewise functions whose pieces are drawn
from ECT-spaces, a natural generalization of algebraic polynomial spaces.
Under suitable assumptions, they admit a basis of functions equipped with all
the nice properties of classical polynomial B-splines, the so-called TB-
splines. Theoretically, TB-splines are completely plug-to-plug compatible with
polynomial B-splines. Unfortunately, TB-splines in their most generality still
lack practical algorithms for their evaluation and manipulation. For this
reason, despite their great potential, they did not gain much attention in
applications so far.

On the contrary, for the restricted class of TB-splines with pieces belonging
to ECT-spaces that are null-spaces of constant-coefficient linear differential
operators, efficient manipulation routines have been developed recently and
made publicly available in a Matlab toolbox [43], so that they can be easily
incorporated in any software library supporting polynomial B-splines to enrich
its capability. This subclass of TB-splines already provides a large variety
of combinations of polynomial, exponential, and trigonometric functions
equipped with a wide spectrum of shape parameters. They allow for an exact
representation of profiles of interest in applications, they behave nicely
with respect to differentiation and integration and, by construction, they
include fundamental solutions of interesting differential operators.
Therefore, they offer a valid alternative to classical polynomial B-splines
and NURBS in Galerkin isogeometric methods.

It turns out that TB-splines can outperform polynomial B-splines whenever
appropriate problem-driven selection strategies for the underlying ECT-spaces
are applied. Typically, there are two aspects that can be considered to
properly select and profit from a non-polynomial structure in the reference
ECT-spaces: the geometric aspect of the problem based on an exact geometry
mapping by using suitable ECT-spaces, and the analytical behavior of the
solution that can be better captured, avoiding spurious oscillations, by
including in the discretization spaces functions related to fundamental
solutions of the given differential operator.

The proposed setting is promising but also presents some theoretical,
computational, and numerical issues that require more attention and deserve
further investigation in future research.

  * 1.

The existence of a TB-spline basis is ensured when the different pieces are
taken from a $(p+1)$-dimensional ECT-space on the entire interval $[a,b]$
which contains constants and its derivative space is a $p$-dimensional ECT-
space on the same interval. Since our interest is confined to null-spaces of
constant-coefficient linear differential operators, we can weaken this
restriction and we only need to check that the critical length of the
derivative space we are dealing with exceeds $\xi_{k+p}-\xi_{k+1}$ for
$k=1,\ldots,n$. Besides the results concerning the critical length recalled in
Sections 2.1 and 2.2, we remark that a numerical procedure to estimate the
critical length for a given ECT-space is provided in [4]. In the more general
setting where the various pieces of the Tchebycheffian splines are drawn from
different ECT-spaces, the existence of a TB-spline basis requires some
theoretical assumptions. A numerical characterization for the existence can be
found in [3]. Alternatively, for practical purposes one can simply use the
Matlab toolbox in [43] and directly check whether the produced set of
functions possess the properties of interest (non-negativity and partition of
unity).

  * 2.

The Matlab toolbox in [43] provides efficient evaluation routines for a wide
range of TB-splines, with pieces belonging to ECT-spaces that are null-spaces
of constant-coefficient linear differential operators, of interest in
practical applications. However, as discussed and illustrated in [43, Section
6.2], possible numerical instabilities can arise when extreme values for shape
parameters (i.e., of the roots of the characteristic polynomial (8)), high
degrees, high smoothness, and highly non-uniform partitions are considered. A
“fidelity” check for the obtained basis is recommended when dealing with
extreme configurations. Such choices are anyway less relevant in most
practical applications. Nevertheless, providing more robust implementations is
an interesting but challenging line of further research.

  * 3.

As for any Galerkin method, accurate quadrature rules are imperative for
isogeometric Galerkin methods based on TB-splines. When using classical
elementwise Gaussian quadrature, the same order of quadrature used for
polynomial B-splines produces satisfactory results only for mild values of the
shape parameters. For larger values of the shape parameters, however, the
amount of quadrature points needs to be significantly increased (up to five
times). The construction of tailored quadrature rules for TB-splines, in the
spirit of [12] and references therein, is definitely worth to be investigated.

Finally, we observe that, due to the complete structural similarity with
polynomial B-splines, the most popular local tensor-product structures
supporting local refinement (such as T-splines [42], (truncated) hierarchical
B-splines [18], and LR-splines [16]) can be easily extended to the
Tchebycheffian setting; see [7] for an example.

Although the large variety of ECT-spaces that are null-spaces of constant-
coefficient linear differential operators already offers an extreme flexible
environment for applications, an additional extension is provided by the so-
called multi-degree Tchebycheffian B-splines (MDTB-splines); see [19, 38, 43].
MDTB-splines allow for ECT-spaces of different dimensions on different
intervals of the given partitions; they are also supported by the Matlab
toolbox in [43]. However, finding (automatic) strategies to select both the
dimension and the structure for different ECT-spaces seems to be a very
challenging task.

## Acknowledgements

This work was supported by the MSCA-ITN-2019 program through the project
GRAPES (contract n. 860843 – CUP E54I19002020006) and by the MIUR Excellence
Department Project awarded to the Department of Mathematics, University of
Rome Tor Vergata (CUP E83C18000100006). The last two authors are members of
Gruppo Nazionale per il Calcolo Scientifico, Istituto Nazionale di Alta
Matematica.

## References

  * [1] A. Aimi, M. Diligenti, M. L. Sampoli, and A. Sestini, _Non-polynomial spline alternatives in isogeometric symmetric Galerkin BEM_ , Applied Numerical Mathematics 116 (2017), 10–23.
  * [2] C. Bangert and H. Prautzsch, _Circle and sphere as rational splines_ , Neural, Parallel & Scientific Computations 5 (1997), 153–162.
  * [3] C. V. Beccari, G. Casciola, and M. L. Mazure, _Design or not design? A numerical characterisation for piecewise Chebyshevian splines_ , Numerical Algorithms 81 (2019), 1–31.
  * [4] , _Critical length: An alternative approach_ , Journal of Computational and Applied Mathematics 370 (2020), 112603.
  * [5] C. V. Beccari, G. Casciola, and L. Romani, _A practical method for computing with piecewise Chebyshevian splines_ , Journal of Computational and Applied Mathematics 406 (2022), 114051\.
  * [6] D. Bister and H. Prautzsch, _A new approach to Tchebycheffian B-splines_ , Curves and Surfaces with Applications in CAGD (A. Le Méhauté, C. Rabut, and L. L. Schumaker, eds.), Vanderbilt University Press, Nashville, 1997, pp. 387–394.
  * [7] C. Bracco, T. Lyche, C. Manni, F. Roman, and H. Speleers, _Generalized spline spaces over T-meshes: Dimension formula and locally refined generalized B-splines_ , Applied Mathematics and Computation 272 (2016), 187–198.
  * [8] , _On the dimension of Tchebycheffian spline spaces over planar T-meshes_ , Computer Aided Geometric Design 45 (2016), 151–173.
  * [9] M. Brilleaud and M. L. Mazure, _Design with L-splines_ , Numerical Algorithms 65 (2014), 91–124.
  * [10] A. N. Brooks and T. J. R. Hughes, _Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations_ , Computer Methods in Applied Mechanics and Engineering 32 (1982), 199–259.
  * [11] B. Buchwald and G. Mühlbach, _Construction of B-splines for generalized spline spaces generated from local ECT-systems_ , Journal of Computational and Applied Mathematics 159 (2003), 249–267.
  * [12] F. Calabrò, G. Loli, G. Sangalli, and M. Tani, _Quadrature rules in the isogeometric Galerkin method: State of the art and an introduction to weighted quadrature_ , Advanced Methods for Geometric Modeling and Numerical Simulation (C. Giannelli and H. Speleers, eds.), Springer INdAM Series, vol. 35, Springer International Publishing, 2019, pp. 43–55.
  * [13] J. M. Carnicer, E. Mainar, and J. M. Peña, _Critical length for design purposes and extended Chebyshev spaces_ , Constructive Approximation 20 (2003), no. 1, 55–71.
  * [14] , _Critical lengths of cycloidal spaces are zeros of Bessel functions_ , Calcolo 54 (2017), 1521–1531.
  * [15] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs, _Isogeometric analysis: Toward integration of CAD and FEA_ , John Wiley & Sons, 2009.
  * [16] T. Dokken, T. Lyche, and K. F. Pettersen, _Polynomial splines over locally refined box-partitions_ , Computer Aided Geometric Design 30 (2013), 331–356.
  * [17] N. Dyn and A. Ron, _Recurrence relation for Tchebycheffian B-splines_ , Journal d’Analyse Mathématique 51 (1988), 118–138.
  * [18] C. Giannelli, B. Jüttler, and H. Speleers, _THB-splines: The truncated basis for hierarchical splines_ , Computer Aided Geometric Design 29 (2012), 485–498.
  * [19] R. R. Hiemstra, T. J. R. Hughes, C. Manni, H. Speleers, and D. Toshniwal, _A Tchebycheffian extension of multidegree B-splines: Algorithmic computation and properties_ , SIAM Journal on Numerical Analysis 58 (2020), 1138–1163.
  * [20] T. J. R. Hughes, J. A. Cottrel, and Y. Basilev, _Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement_ , Computer Methods in Applied Mechanics and Engineering 194 (2005), 4135–4195.
  * [21] T. J. R. Hughes, L. P. Franca, and G. M. Hulbert, _A new finite element formulation for computational fluid dynamics: VIII. the Galerkin/least-squares method for advective-diffusive equations_ , Computer Methods in Applied Mechanics and Engineering 73 (1989), 173–189.
  * [22] S. Karlin, _Total positivity_ , Stanford University Press, Stanford, 1968.
  * [23] S. Karlin and W. J. Studden, _Tchebycheff systems: With applications in analysis and statistics_ , Interscience Publishers, New York, 1966.
  * [24] T. Lyche, _A recurrence relation for Chebyshevian B-splines_ , Constructive Approximation 1 (1985), 155–173.
  * [25] T. Lyche, C. Manni, and H. Speleers, _Foundations of spline theory: B-splines, spline approximation, and hierarchical refinement_ , Splines and PDEs: From Approximation Theory to Numerical Linear Algebra (T. Lyche, C. Manni, and H. Speleers, eds.), Lecture Notes in Mathematics, vol. 2219, Springer International Publishing, 2018, pp. 1–76.
  * [26] , _Tchebycheffian B-splines revisited: An introductory exposition_ , Advanced Methods for Geometric Modeling and Numerical Simulation (C. Giannelli and H. Speleers, eds.), Springer INdAM Series, vol. 35, Springer International Publishing, 2019, pp. 179–216.
  * [27] T. Lyche and L. L. Schumaker, _Total positivity properties of LB-splines_ , Total Positivity and Its Applications (M. Gasca and C. A. Micchelli, eds.), Mathematics and Its Applications, vol. 359, Klüwer, 1996, pp. 35–46.
  * [28] C. Manni, F. Pelosi, and M. L. Sampoli, _Generalized B-splines as a tool in isogeometric analysis_ , Computer Methods in Applied Mechanics and Engineering 200 (2011), 867–881.
  * [29] , _Isogeometric analysis in advection–diffusion problems: Tension splines approximation_ , Journal of Computational and Applied Mathematics 236 (2011), 511–528.
  * [30] C. Manni, F. Pelosi, and H. Speleers, _Local hierarchical $h$-refinements in IgA based on generalized B-splines_, Mathematical Methods for Curves and Surfaces 2012 (M. Floater, T. Lyche, M. L. Mazure, K. Mørken, and L. L. Schumaker, eds.), Lecture Notes in Computer Science, vol. 8177, Springer–Verlag, Heidelberg, 2014, pp. 341–363.
  * [31] C. Manni, A. Reali, and H. Speleers, _Isogeometric collocation methods with generalized B-splines_ , Computers & Mathematics with Applications 70 (2015), 1659–1675.
  * [32] C. Manni, F. Roman, and H. Speleers, _Generalized B-splines in isogeometric analysis_ , Approximation Theory XV: San Antonio 2016 (G. E. Fasshauer and L. L. Schumaker, eds.), Springer Proceedings in Mathematics & Statistics, vol. 201, Springer International Publishing, 2017, pp. 239–267.
  * [33] M. L. Mazure, _Chebyshev spaces and Bernstein bases_ , Constructive Approximation 22 (2005), 347–363.
  * [34] , _Extended Chebyshev piecewise spaces characterised via weight functions_ , Journal of Approximation Theory 145 (2007), 33–54.
  * [35] , _Finding all systems of weight functions associated with a given extended Chebyshev space_ , Journal of Approximation Theory 163 (2011), 363–376.
  * [36] , _How to build all Chebyshevian spline spaces good for geometric design?_ , Numerische Mathematik 119 (2011), 517–556.
  * [37] , _Constructing totally positive piecewise Chebyshevian B-spline bases_ , Journal of Computational and Applied Mathematics 342 (2018), 550–586.
  * [38] G. Nürnberger, L. L. Schumaker, M. Sommer, and H. Strauss, _Generalized Chebyshevian splines_ , SIAM Journal on Mathematical Analysis 15 (1984), 790–804.
  * [39] J. Sánchez-Reyes, _Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials_ , Computer Aided Geometric Design 15 (1998), 909–923.
  * [40] I. J. Schoenberg, _On spline functions_ , Inequalities (O. Shisha, ed.), Academic Press, New York, 1967, pp. 255–286.
  * [41] L. L. Schumaker, _Spline functions: Basic theory_ , third ed., Cambridge University Press, Cambridge, 2007.
  * [42] T. W. Sederberg, J. Zheng, A. Bakenov, and A. Nasri, _T-splines and T-NURCCs_ , ACM Transactions on Graphics 22 (2003), 477–484.
  * [43] H. Speleers, _Algorithm 1020: Computation of multi-degree Tchebycheffian B-splines_ , ACM Transactions on Mathematical Software 48 (2022), 1–31.
  * [44] H. Speleers, C. Manni, F. Pelosi, and M. L. Sampoli, _Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems_ , Computer Methods in Applied Mechanics and Engineering 221–222 (2012), 132–148.