\section{Introduction}
Electromagnetic field computations and, more generally, the numerical discretization of equations enjoying a relevant geometric structure, is one of the most interesting challenge of numerical analysis for PDEs  and several results have been obtained in the last decade. Indeed, only for Galerkin methods, three Acta Numerica overview papers have been  published:  by Hipmair \cite{HIP02a}, by Arnold, Falk and Winther \cite{AFW06}, and by Boffi \cite{Boff10}, addressing different aspects of the problem.

 On the one hand, discrete schemes have to preserve the geometric structure of the underlying PDEs in order to avoid spurious behaviors, instability or non-physical solutions (see e.g., the pioneering paper \cite{BFGP}). \color{black} For electromagnetics, as it is clear from the references above, numerical schemes have to be related with a discrete De Rham complex. \color{black} On  the other hand, especially in view of high frequency computations, numerical schemes need to be efficient and accurate. This requires many features, and among others it requires adaptivity, or at least local mesh refinement capability, in order to capture the strong singularities of the electromagnetic field,   possibly driven by \emph{a-posteriori} error estimator as, e.g., in \cite{BrSc08}.

In this paper we present and analyse discretization techniques for electromagnetic fields based on splines and  generalizations of splines, as NURBS (\cite{Piegl}) or T-splines (\cite{Sederberg_Zheng} or below). Our work originates from IsoGeometric Analysis (IGA),  \cite{IGA-book}. Isogeometric analysis has been introduced in 2005 by Hughes and co-authors in the seminal paper \cite{HCB05} to solve structural mechanic problems directly on the geometry output by a CAD system, and has set the paradigm to use splines, NURBS or their generalization as generating functions for the construction of Galerkin spaces. This idea has been proved to be extremely effective and IGA is spreading very fast across different scientific communities: structural mechanics (see e.g., \cite{IGA-book}, \cite{MR2443159}, \cite{Cottrell_Reali_Bazilevs_Hughes}, \cite{Echter2010374}, \cite{Lipton_Evans_Bazilevs}, \cite{rank12}, \cite{MR2835758}), geometric modeling (see e.g., \cite{DQS10}, \cite{Martin_Cohen_Kirby}, \cite{CoMa10} and also \cite{Zhang2012}) and numerical analysis (see e.g., \cite{BBCHS06}, \cite{Beirao_Cho_Sangalli} , \cite{BRSV11}, \cite{nwidth-iga}, \cite{Vuong_giannelli_juttler_simeon}, \cite{BePa12}, \cite{kleissieti}).

 In this paper we present the recent advances in the use of the isogeometric paradigm and  spline-based  methods for  electromagnetic wave computations. This research has started with the two papers \cite{Buffa_Sangalli_Vazquez} and \cite{BRSV11} and can likely be considered as still in infancy (see also \cite{ratnani} for the applications of this results). This paper aims at showing the potential impact  of these techniques in the electromagnetic community by addressing several aspects: from the geometric structure of the proposed methods, to local refinement strategies.

  We introduce the \textbf{spline complex} studied in \cite{BRSV11} (see \eqref{eq:3D-parametric-diagram-bc} and \eqref{eq:diag-discrete}) and we present its properties: we construct canonical bases so that the matrices representing differential operators are the incidence matrices of the underlying meshes, and this enlightens the relation between the spline complex and the geometry of the underlying meshes. We show that for different choices of the  degree of splines, the spline complex is isomorphic to the co-chain complex or to the chain complex of the underlying mesh.
Besides this interesting fact, we also introduce the concept of \emph{control fields} in analogy to control points which are ubiquitous in spline theory (see e.g., \cite{cohen2001geometric} or \cite{DeBoor}) which provide the correct physical interpretation of degrees of freedom.
Finally, we extend the results of \cite{BRSV11} to \textbf{multi-patch} geometries, i.e., geometries which are piecewise spline or NURBS mapping of the unit cube. We refer the reader to \cite{kleissieti} for a detailed description of this class of geometries.

The second major contribution in this paper is a step towards
\textbf{adaptivity} for spline-based methods. Leveraging on the recent
work on T-splines, we design a two dimensional \textbf{T-spline
  complex} where meshes with T-junctions can be used to allow for
adaptivity.  T-splines are  the most attractive way to break the
tensor product structure of splines while keeping their structure and
their accuracy. T-splines have been introduced in
\cite{Sederberg_Zheng} and \cite{Sederberg_Cardon} and their use as a fundamental tool to enhance isogeometric analysis with adaptivity has been proposed in \cite{Bazilervs_Calo_Cottrell_Evans}.  A series of papers has followed \cite{Buffa_Cho_Sangalli}, \cite{Scott_Li_Sederberg_Hughes}, \cite{Wang2011477}, \cite{Beirao_Buffa_Cho_Sangalli},  together with the relevant class of Analysis Suitable T-splines \cite{LZSHS12}, \cite{BBCS12} which we use in our construction. The two dimensional T-spline complex is used to treat three dimensional problems with symmetry. We should also mention that  the definition and use of T-splines in three dimension are not yet well understood, but object of an intensive study. Their use will allow, on a longer time perspective, to design full adaptive algorithms, on very general geometries parametrized on  totally unstructured meshes. We refer the reader to \cite{Scott} for a monograph on the modern use of T-splines in geometric modeling.

\color{black} Finally, we should remark that the spline spaces we study in this
paper have a wide domain of applications and can be applied
successfully to the discretization of  other problems than
electromagnetics. In fact, they can be used to solve the Darcy flows
equations or more generally the Hodge laplacian operator as detailed
in \cite{AFW06} and  \cite{AFW-2}. Moreover, thanks to the regularity
of spline spaces, their use in fluids is very promising. In the paper
\cite{Buffa_deFalco_Sangalli} they are used for the first time to
solve the Stokes equations, \color{black} in  \cite{EvHu12-4} the
Stokes eigenvalue problem  is addressed, and \color{black} in the sequence of three papers \cite{EvHu12}, \cite{EvHu12-2}, \cite{EvHu12-3} they are applied to solve Stokes and Brinkman equations, steady and unsteady Navier-Stokes equations, providing impressive results.
\color{black}

The outline of the paper is the following. In Section \ref{sec:notation} we set up the notation for the problems we address, in Section \ref{sec:prel-splin-nurbs} we present known results about splines and NURBS in a self-contained way; in Section \ref{sec:spline-complex} we present the spline complex and all the related results while in Section \ref{sec:TT} we introduce the T-spline complex and analyse its properties. Finally, in Section \ref{sec:num} we present numerical results: the first ones are two and three dimensional, academic tests aiming at demonstrating the validity of the proposed approach. As a last example, we compute the propagation in  a waveguide with geometric inhomogeneity, on a three dimensional locally refined mesh.

\section{Notation}
\label{sec:notation}

In this section we present the notation that we need to describe
the time-harmonic Maxwell problem.
Let $\Omega \subset \mathbb{R}^3$ be a bounded Lipschitz domain. We denote by
$L^2(\Omega)$ the space of complex square integrable functions on
$\Omega$, endowed with standard $L^2$ norm $\| \cdot
\|_{L^2(\Omega)}$, and by ${\bf L}^2(\Omega)$ their vectorial
counterparts. The Hilbert space $H^1(\Omega)$ contains functions
of $L^2(\Omega)$ such that their first order derivatives also belong
to $L^2(\Omega)$. We denote by $H^1_0(\Omega) \subset H^1(\Omega)$ the subspace of functions with homogeneous boundary condition.
 We will also make use of the space $\,{\bf H}({\mathbf{curl}};\Omega)\,$, constituted by all
 functions in ${\bf L}^2(\Omega)$ such that their curl also belongs to
 ${\bf L}^2(\Omega)$, and $\,{\bf H}({\rm div};\Omega)\,$, the space of functions in
 ${\bf L}^2(\Omega)$ such that their divergence belongs to $L^2(\Omega)$.
 Moreover, we denote by $\,{\bf H}_0(\mathbf{curl};\Omega)\,$
 (resp. $\,{\bf H}_0({\rm div};\Omega)\,$) the subspace of $\,{\bf H}({\mathbf{curl}};\Omega)\,$
 (resp. $\,{\bf H}({\rm div};\Omega)\,$) of fields with vanishing tangential
 (resp. normal) component.

For the sake of simplicity, we assume that the domain $\Omega$, referred to as physical domain in
the following,  is bounded Lipschitz and simply connected, and that
its boundary $\partial \Omega$ is connected. We also assume that it is
defined through a \color{black} continuously differentiable \color{black}  parametrization with \color{black} continuously differentiable \color{black}  inverse which we
denote as  ${\bf F}: \widehat{\Omega} \longrightarrow \Omega$, where $\widehat{\Omega}$
will be referred to as the \emph{parametric domain}. Further assumptions on the
geometrical mapping ${\bf F}$ will be given later.

Some notation will be borrowed from the context of differential forms:
first of all, we define the spaces
 \begin{eqnarray*}
\hat{X}^0 := \,{H}^1(\Omegaref)\,, \; \hat{X}^1 := \,{\bf H}({\mathbf{curl}};{\Omegaref})\,, \; \hat{X}^2 := \,{\bf H}({\rm div};{\Omegaref})\,, \; \hat{X}^3 := \,{L}^2(\Omegaref)\, ,\\
X^0 := \,{H}^1(\Omega)\,, \; X^1 := \,{\bf H}({\mathbf{curl}};\Omega)\,, \; X^2 := \,{\bf H}({\rm div};\Omega)\,, \; X^3 := \,{L}^2(\Omega)\,;
\end{eqnarray*}
Since the parametrization ${\bf F}$ and its inverse are smooth, we can
define the pullbacks that relate these spaces as (see
\cite[Sect.~2.2]{HIP02a}):
\begin{equation}
\begin{array}{ll}
 \iota^0(\phi) := \phi \circ {\bf F} , &\quad \phi \in X^0 , \\
 \iota^1(\bu) := (D{\bf F})^T(\bu \circ {\bf F}) , &\quad \bu \in X^1 , \\
 \iota^2(\bv) := \det(D{\bf F}) (D{\bf F})^{-1} (\bv \circ {\bf F}) , & \quad \bv \in X^2 , \\
 \iota^3(\varphi) := \det(D{\bf F}) (\varphi \circ {\bf F}) , & \quad \varphi \in X^3 , \label{def:iota3}
\end{array}
\end{equation}
where $D{\bf F}$ is the Jacobian matrix of the mapping ${\bf F}$. Then, due to the curl and divergence conserving properties of $\iota^1$ and $\iota^2$, respectively (see \cite[Sect.~3.9]{Monk}, for instance), the following commuting De~Rham diagram is satisfied (see \cite[Sect.~2.2]{HIP02a}):
\begin{equation} \label{eq:diag1}
\begin{CD}
\mathbb{R} @>>> \hat{X}^0 @>\hat \,{\bf grad}\,>> \hat{X}^1 @>\hat \,{\bf curl}\,>> \hat{X}^2 @>\hat \,{\rm div}\,>> \hat{X}^3 @>>> 0 \\
@. @A\iota^0AA @A\iota^1AA @A\iota^2AA @A\iota^3AA \\
\mathbb{R} @>>> X^0 @>\,{\bf grad}\,>> X^1 @>\,{\bf curl}\,>> X^2 @>\,{\rm div}\,>> X^3 @>>> 0 . \\
\end{CD}
\end{equation}

We are also interested in  spaces with boundary
conditions,  denoted with the subindex $0$,
\begin{eqnarray*}
\hat{X}^0_0 := \,{H}^1_0(\Omegaref)\,, \; \hat{X}^1_0 := \,{\bf H}_0(\mathbf{curl};{\Omegaref})\,, \; \hat{X}^2_0 := \,{\bf H}_0({\rm div};{\Omegaref})\,, \;
\hat{X}^3_0 :=   \,{L}^2(\Omegaref)\,  , \\
X^0_0 := \,{H}^1_0(\Omega)\,, \; X^1_0 := \,{\bf H}_0(\mathbf{curl};\Omega)\,, \; X^2_0 := \,{\bf H}_0({\rm div};\Omega)\,, \;
X^3_0 := \,{L}^2(\Omega)\, ,
\end{eqnarray*}
for which  the De~Rham diagram reads
\begin{equation} \label{eq:diag2}
\begin{CD}
0 @>>> \hat{X}^0_0 @>\hat \,{\bf grad}\,>> \hat{X}^1_0 @>\hat \,{\bf curl}\,>> \hat{X}^2_0 @>\hat
\,{\rm div}\,>> \hat{X}^3_0 @>\int >> \mathbb{R}  \\
@. @A\iota^0AA @A\iota^1AA @A\iota^2AA @A\iota^3AA \\
0 @>>> X^0_0 @>\,{\bf grad}\,>> X^1_0 @>\,{\bf curl}\,>> X^2_0 @>\,{\rm div}\,>> X^3_0 @> \int>> \mathbb{R} , \\
\end{CD}
\end{equation}
which also expresses the integral preserving property of
$\iota^3$.

\color{black}
\begin{remark}
As it is well known, the exactness of the sequences \eqref{eq:diag1} and \eqref{eq:diag2} relies on the assumption that $\Omega$ (and $\hat\Omega$) has a  trivial topology. All what we develop in this paper applies in principle also to the case of arbitrary topology but we do not present all the details here.
\end{remark}
\color{black}

\section{Preliminaries on splines and NURBS}
\label{sec:prel-splin-nurbs}

We give here a brief overview on B-splines and, in the spirit of \cite{Bazilervs_Calo_Cottrell_Evans}, we also introduce some concepts that will be needed in the definition of T-splines. For more details on B-splines we  refer \color{black}  the reader to \cite{HCB05,Piegl}.

\subsection{Univariate B-splines}
\subsubsection{Knot vector and B-spline functions, refinement, spline derivatives} \label{sect:knot_vector}
Given two positive integers $p$ and $n$, we say that $\Xi := \{\xi_1,\dots, \xi_{n+p+1}\}$ is a $p$-open knot vector if
$$
\xi_1 =\ldots=\xi_{p+1} < \xi_{p+2} \le \ldots \le
\xi_{n} < \xi_{n+1}=\ldots=\xi_{n+p+1},
$$
where repeated knots are allowed and denote by  $m_j$  the multiplicity of the knot
 $\xi_j$. We assume  $ m_j \leq p+1$ for all internal knots.\color{black}

From the knot vector $\Xi$, B-spline functions of degree $p$ are defined following the well-known Cox-DeBoor recursive formula: we start with piecewise constants ($p=0$):
\begin{equation}\label{eq:cox-deboor-1}
N_{i,0}(\zeta) = \left \{
\begin{array}{ll}
1 & \text{if } \xi_i \leq \zeta < \xi_{i+1}, \\
0 & \text{otherwise},
\end{array}
\right.
\end{equation}
and for $p \ge 1$ the \emph{B-spline}  functions are defined by the recursion
\begin{equation}\label{eq:cox-deboor-2}
N_{i,p}(\zeta) = \frac{\zeta - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\zeta) + \frac{\xi_{i+p+1} - \zeta}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\zeta).
\end{equation}
 This gives a set of $n$ B-splines that form a basis of the space of
 \emph{splines}, that is, piecewise polynomials of degree $p$ with
 $p-m_j$ continuous derivatives at the internal knots $\xi_j$, for $j
 = p+2, \ldots, n$.  We denote this univariate spline space by
 \begin{equation}
   \label{eq:spline-univariate-space}
   \sp_p(\Xi) = \text{span} \{ N_{i,p}, \, i=1,\ldots, n \}
 \end{equation}
  An example of some B-splines is given in Figure~\ref{fig:1}.
\begin{figure}[pt!]
\centering
\includegraphics[width=.5\textwidth]{univariate_deg2.pdf} \hspace{-1cm}\includegraphics[width=.5\textwidth]{univariate_deg3.pdf}
\caption{Example of B-splines of degree $2$ (left) and $3$ (right).} \label{fig:1}
\end{figure}
 Notice that the B-spline function $N_{i,p}$ is supported in the
 interval $[\xi_i, \xi_{i+p+1}]$, and in fact its definition only
 depends on the knots within that interval. For this reason, we define
 the \emph{local knot vector} $\Xi_{i,p} = \{\xi_i, \xi_{i+1}, \ldots,
 \xi_{i+p+1}\}$, and we will sometimes denote $N_{i,p}(\zeta) \equiv
 N[\Xi_{i,p}](\zeta)$.

In the context of splines, three kinds of refinement are possible, as explained in \cite{HCB05}:
\begin{enumerate}
\item  \emph{$k$-refinement} which corresponds to successive application of the Cox-DeBoor formula
  \eqref{eq:cox-deboor-1}--\eqref{eq:cox-deboor-2}. Regularity is
  raised together with the degree:  therefore, the spaces
  \eqref{eq:spline-univariate-space} are not nested under
  {$k$-refinement} but, at each step (degree and regularity elevation),
  the dimension of the space increases by $1$.  The name {$k$-refinement} has been introduced in \cite{HCB05};
\item \emph{$h$-refinement} which corresponds to mesh refinement and is obtained by knot insertion. Let $\bar\Xi := \{\xi_1,\dots,\xi_{k}, \bar \xi, \xi_{k+1}, \xi_{n+p+1}\}$ be the knot vector after inserting the knot $\bar \xi$ in $\Xi$. Then, the new B-spline functions $\{ \bar N_{1,p} (\zeta),\ldots, \bar N_{n+p+2,p} (\zeta)\} $  are constructed as follows:
  \begin{equation}
    \label{eq:href}
  \bar N_{i,p} (\zeta)  = \alpha_i N_{i,p} (\zeta) + (1-\alpha_i) N_{i-1,p} (\zeta)
  \end{equation}
  where  $\alpha_i =1  $, if  $1\leq i\leq k-p$, $\alpha_i =
  \frac{\displaystyle \strut  \bar \xi - \xi_i }{\displaystyle \strut
    \xi_{i+p} - \xi_i}$, if $k-p+1\leq i\leq k$ and, $\alpha_i =0$ for
  $k+1\leq i \leq n+p+2$.  When $\bar \xi$ is  equal to $\xi_k$ or
  $\xi_{k+1}$ or to both, the  knot insertion  corresponds to reduction of the inter-element regularity at $\bar \xi $. \color{black}
\item \emph{$p$-refinement} which corresponds to the degree raising with fixed interelement regularity, and generates a sequence of nested spaces.
\end{enumerate}

Assuming the maximum multiplicity of the internal knots is less than
or equal to $p$,  i.e., the B-spline functions are at least continuous, the derivative of the B-spline $N_{i,p}$ is a spline
as well. Indeed, it belongs to the spline  space  $S_{p-1}(\Xi')$,
where $\Xi' = \{ \xi_2, \ldots, \xi_{n+p} \}$ is a ($p-1$)-open knot
vector. Obviously, the regularity of splines in $\sp_{p-1}(\Xi')$ is one less than the
regularity in $\sp_p(\Xi)$.

 In the following we assume that $\xi_1 = 0$ and $\xi_{n+p+1} = 1$. The domain  $(0,1)$ of definition of the spline functions is the one-dimensional  \emph{parametric} domain. On it, the knot vector
$\Xi$ induces a partition  of the interval $(0,1)$ that we denote by
$\mathcal{M}$.  Precisely, we define  $\mathcal{M}$ as the set of the knot spans
$(\xi_i,\xi_{i+1})$, $i=\lceil p/2\rceil +1 ,...,n + \lfloor
p/2\rfloor$, that can also be empty due to knot multiplicity greater
than $1$. Empty intervals still play a role in the definition of
B-splines and are graphically represented as points close one to the
other, as proposed in \cite{Buffa_Cho_Sangalli}.
Note that in this representation of $\mathcal{M}$, the number of
lines is the knot multiplicity with one exception: for each
boundary knot (at 0 or 1) of an open knot vector in $\mathcal{M}$ we represent
only a  multiplicity  of  $\lfloor p/2\rfloor + 1$, which is  $(p+1)/2$
lines if $p$ is odd, and $p/2+1$ lines if $p$ is even (see Figure~\ref{fig:1}). \color{black}  The reason for
this construction of $\mathcal{M}$ will be motivated in the next
section.

 Finally, it is worth  noting the relationship between the space
$\sp_p(\Xi)$ and the space of derivatives $\sp_{p-1}(\Xi')$, and their
respective meshes $\mathcal{M}$ and $\mathcal{M}'$. The meshes $\mathcal{M}$ and $\mathcal{M}'$ may differ only as regards the number of points at the boundary.
Indeed, according to the definition above,  if $p$ is odd both meshes
coincide, and if $p$ is even the number of elements of $\mathcal{M}'$ with
respect to $\mathcal{M}$ is reduced by two, one on each side.
\color{black}

\subsubsection{Anchors and Greville sites}\label{sect:anchors}

In this section we present the concept of anchors and of Greville sites as points
in the parametric space $(0,1)$  which may be associated to each B-spline function.
Greville sites,  which are also known as knot averages, \color{black} are classical and can be found for instance in \cite{DeBoor},
while the concept of anchors has been introduced recently in
\cite{Bazilervs_Calo_Cottrell_Evans}.

Since splines are not interpolatory, the association of functions to
points (or, as we will see, other geometric entities)  is somehow more
arbitrary than with Lagrangian  finite elements. Anchors and Greville
sites are  two different choices, and we present here both.  Anchors
are very useful
when dealing with non-tensor product extensions of splines as
T-meshes, while Greville sites (and related geometric entities) carry

degrees of freedom in a more natural way.

  Given a B-spline function  $N_{i,p}(\zeta)$,  and its local knot vector  $\Xi_{i,p} =\{  \xi_i, \xi_{i+1}, \ldots, \xi_{i+p+1}\}$, we set: if $p$ is odd, the anchor $A$   associated with $N_{i,p}(\zeta)$ is  the central knot of $\Xi_{i,p} $. If $p$ is even, the anchor $A$   associated with $N_{i,p}(\zeta)$  is chosen to be the midpoint of the central knot span of $\Xi_{i,p} $, namely: $\zeta^A := \frac{\displaystyle \strut \xi_{i+p/2} + \xi_{i + p/2+1}}{\displaystyle \strut 2}$. The position of the anchors for degrees $p=2$ and $p=3$ are represented in Figure~\ref{fig:1}.

Note that obviously the correspondence between anchors and B-splines functions is one to one, but different anchors $A\neq \bar{A}$  may lie at the same position $\zeta^A=\zeta^{\bar{A}}$. A remedy to this abuse of notation,
at the cost of more complex definition,  is proposed  in \cite{BBSV12}
where  the use of both an \emph{index} and a \emph{parametric} domain is proposed.

The set of anchors is denoted as $\mathcal{A}_p=\mathcal{A}_p(\Xi)$. When $p$ is odd
anchors are located at all knots \color{black}   of the partition $\mathcal{M}$ (which
may be repeated), while when $p$ is even  anchors are located at midpoints
 of all elements in $\mathcal{M}$ (including the ones of zero
area). Indeed, this fact is the reason for the definition of $\mathcal{M}$ in
particular as regards to the multiplicity of boundary knots.\color{black}

Most often, we will use anchors to index functions and local knot vectors. Namely, for an anchor $A\in \mathcal{A}_p$,  $\Xi^A_p$ and $B^A_p(\zeta)$ will denote the corresponding local knot vector and B-spline function, respectively.  When no confusion occurs, the subindex may be removed.

\begin{remark}\label{rem:interpolatory-B-splines}
The B-splines are, in general,  not interpolatory at the anchor $A \in \mathcal{A}_p(\Xi)$, while they are interpolatory at knots having multiplicity $p$.
This always happens at $\zeta=0$ and $\zeta=1$, and happens in the interior of
the parametric domain where the basis is $C^0$ continuous, i.e., at knots with multiplicity $p$.  See e.g., Figure~\ref{fig:1}(left).
\end{remark}

Given $A\in \mathcal{A}_p$,  and $\Xi^A_p =\{  \xi_i,.., \xi_{i+p+1} \}$ for some $i$, then the \emph{Greville site} is defined as:
\begin{equation}
  \label{eq:greville1D}
  \gamma^A = \frac{\xi_{i+1} + \ldots + \xi_{i+p}}{p}.
\end{equation}
Unlike anchor positions, Greville sites are all different one from the
other, when the multiplicities $m_j$ verify $m_j\leq p$ and thus
B-splines are all continuous. The Greville sites induce a partition
of the unit interval, referred as Greville mesh and denoted $\mathcal{M}_G$.
These concepts are ubiquitous in spline theory and geometry representation. Greville sites are intimately related to
 \emph{control points} and \emph{control polygon} whose
properties we briefly recall in the next section.

\subsubsection{B-spline curves}
A B-spline curve $\Gamma$  in $\mathbb{R}^3$ is defined by a parametrization in the interval $(0,1)$, in the form
\begin{equation}\label{eq:Bspline-curve}
{\bf F}(\zeta) = \sum_{A \in \mathcal{A}_p} {\bf C}^A B^A_p(\zeta), \qquad 0 < \xi < 1,
\end{equation}
where ${\bf C}^A \in \mathbb{R}^3$ are called the \emph{control points}. Control
points are in a one-to-one correspondence with B-spline basis
functions. The piecewise linear interpolation of the control points
gives the control polygon $\Gamma_C$. See Figure~\ref{fig:ctrl_polygon} for an example.

The control points ${\bf C}^A $ have an important role not only in the
definition of the spline parametrization \eqref{eq:F-Bspline}, but also in
the visualization and interaction with spline geometries within CAD
softwares.
Indeed, it is common in CAD softwares to represent, together with the
parametrized curve  $\Gamma$, the control points ${\bf C}^A$ and the
associated control polygon $\Gamma_C$. Typically,
the CAD user defines or interacts with  the control points in order to
input and modify the geometry.
Since the B-splines are not in general interpolatory (recall Remark
\ref{rem:interpolatory-B-splines}), then the control polygon
$\Gamma_C$  differs from
$\Gamma$, but it is ``close'' to it.   Precisely,  $\Gamma_C$
converges to $\Gamma$  under $h$-refinement.
This convergence is proved, e.g., in \cite{DeBoor}  and
discussed here below.

We introduce the usual Lagrangian basis for piecewise linear polynomials on the Greville mesh $\mathcal{M}_G$, denoted by $\lambda^A(\cdot)$:
\begin{equation}\label{eq:lambda-1D}
\lambda^A(\gamma^{A'})=\left \{
  \begin{aligned}
    1&\text{ if } A=A',\\
    0&\text{ if } A\neq A'.\\
  \end{aligned}
\right .
\end{equation}
The control polygon $\Gamma_C$  is then parametrized by  the mapping ${\bf F}_C:[0,1] \rightarrow
\Gamma_C $ defined by
\begin{equation}\label{eq:Fh-Bspline-1D}
 {\bf F}_C(\zeta) = \sum_{A \in \mathcal{A}_{p}} {\bf C}^A
\lambda^A(\zeta), \qquad 0 < \xi < 1,
\end{equation}
\color{black}
that is, ${\bf F}_C$ and  ${\bf F}$ share the same
control points.  When ${\bf F}$ is smooth enough, the following approximation estimate holds (see, e.g.,
\cite[Ch. XI]{DeBoor}):
\begin{equation}
  \label{eq:-control-polygon-convergence}
  \underset{\zeta \in [0,1] }{\sup} \| {\bf F}(\zeta)  -   {\bf F}_C(\zeta) \| \simeq h^2,
\end{equation}
$h$ denoting the mesh-size. In other words, $\Gamma_C$ approximates
$\Gamma$ up to an error $O(h^2)$ under $h$-refinement. A graphical representation of this convergence can be seen in Figure~\ref{fig:ctrl_polygon}.

\begin{figure}[pt!]
\centering
\includegraphics[width=.48\textwidth]{ctrl_polygon1.pdf} \includegraphics[width=.48\textwidth]{ctrl_polygon2.pdf} \vspace{-1cm}
\caption{B-spline curve and its control polygon (left), and the same curve after one step of $h$-refinement (right).} \label{fig:ctrl_polygon}
\end{figure}

\color{black}

\subsection{Multivariate B-splines}\label{sect:bivariate}
Multivariate B-splines are defined from univariate B-splines by tensor
product, see for instance \cite{Piegl,DeBoor}. Anchors are defined
in a similar way. We give here a quick overview.

\subsubsection{Knot vectors, B-spline functions, anchors, Greville sites}\label{sect:b-anchors}
Let $d$ be the space dimensions (in practical cases,
$d=2,3$). Assume  $n_\ell \in \mathbb N$, the degree $p_\ell \in \mathbb N$ and  the
\mbox{$p_\ell$-open} knot vector  $\Xi_\ell = \{\xi_{\ell,1}, \ldots,
\xi_{\ell,n_\ell + p_\ell + 1} \}$ are given, for $\ell =
1,\ldots,d$. These  knot vectors define a tensor product mesh $\mathcal{M}$ in
the parametric domain $\hat\Omega = (0,1)^d$
where, as in Section
\ref{sect:anchors}, we have to  take into account the knot
multiplicity.  The multiplicity of a knot vector in $\Xi_\ell$ is
represented graphically by lines ($d=2$) or planes ($d=3$) one close to the other, while the boundary is treated exactly as in one dimension.

The set of anchors is defined on $\mathcal{M}$ as the Cartesian product
$\mathcal{A}_{p_1,\ldots,p_d}(\Xi_1,\ldots,\Xi_d) = \mathcal{A}_{p_1}(\Xi_1) \times
\ldots \times  \mathcal{A}_{p_d}(\Xi_d)$. Considering, for example, the
trivariate case ($d=3$) and recalling the definitions from Section
\ref{sect:anchors} for the univariate case, we have that: if all
$p_\ell$  are odd the anchors lie at the vertices of the mesh; if both
$p_1$ and $p_2$ are odd and $p_3$ is even, then the anchors are
middle-points of the vertical edges  of  $\mathcal{M}$;  if both $p_1$ and $p_2$ are even and $p_3$ is odd, then the anchors are
centers of the horizontal faces of  $\mathcal{M}$;   if  all $p_\ell$  are
even the anchors lie at the center of the elements  of  $\mathcal{M}$, and so
on.  As in the univariate case, the anchors may be located at the
center of zero length edges or zero area faces or empty elements,
according to  knot repetition.
Also,  the computation of the local knot vectors for each anchor
follows from the univariate case.  Given an anchor $A =(A_1, A_2, A_3)\in
\mathcal{A}_{p_1,p_2,p_3} \equiv \mathcal{A}_{p_1,p_2,p_3}(\Xi_1,\Xi_2,\Xi_3)$, we have that its coordinates are
$(\zeta_1^A, \zeta_2^A, \zeta_3^A) = (\zeta^{A_1}, \zeta^{A_2},
\zeta^{A_3})$.
The three  local knot vectors (one in each coordinate direction)  corresponding to $A$ are defined as $\Xi_i^A := [\Xi_i]_{p_i}^{A_i}$ for $i=1,2,3$ and
 the B-spline  associated to $A$ is constructed by tensor product as:
\begin{equation}\label{eq:Bsplines-trivariate}
B^A_{p_1,p_2,p_3}({\boldsymbol{\zeta}}) = B^{A_1}_{p_1}(\zeta_1)  B^{A_2}_{p_2}(\zeta_2)  B^{A_3}_{p_3}(\zeta_3).
\end{equation}
with ${\boldsymbol{\zeta}} = (\zeta_1, \zeta_2, \zeta_3) \in \hat\Omega = (0,1)^3$.

The B-spline functions \eqref{eq:Bsplines-trivariate} span the space $\sp_{p_1,p_2,p_3}(\Xi_1,\Xi_2,\Xi_3)$ (or simply
$\sp_{p_1,p_2,p_3}$), which is the {\color{black} space of piecewise
polynomials of degree $p_\ell$ in the $x_\ell$ direction on $\mathcal{M}$}, whose continuity at the internal mesh plane $\zeta_\ell =
\xi_{\ell,k}$  is $C^{p_\ell-m_{\ell,k}}$,  $m_{\ell,k} $ being the
multiplicity of $\xi_{\ell,k} $ in the knot vector $ \Xi_{\ell}$.

To each anchor $A \in\mathcal{A}_{p_1,p_2,p_3}(\Xi_1,\Xi_2,\Xi_3)$ (or, equivalently, to each B-spline
function \eqref{eq:Bsplines-trivariate}) we also associate a Greville
site in the natural way, that is
\begin{equation}
  \label{eq:greville3D}
 \boldsymbol  \gamma^A = (\gamma_1^A,\gamma_2^A,\gamma_3^A)
\end{equation}
where each $ \gamma_i^A$ is defined as in  \eqref{eq:greville1D}, from
the local knot vector  $\Xi^A_i$.
Connecting adjacent  Greville sites,  we obtain the \emph{Greville
  mesh} $ \mathcal{M}_G$, which is a regular tensor product  mesh with
all elements of positive volume

\subsubsection{Spline and NURBS geometries, multi-patch domains}

\label{sec:spline-geometries}
Analogously to spline curves, a trivariate single-patch spline parametrization of the domain $\Omega \subset\mathbb{R}^3$ is  ${\bf F}:\hat  \Omega \rightarrow
\Omega $ defined as a linear combination of B-splines,
\begin{equation}\label{eq:F-Bspline}
{\bf F}({\boldsymbol{\zeta}}) = \sum_{A \in \mathcal{A}_{p_1,p_2 ,p_3}} {\bf C}^A
B^A_{p_1,p_2,p_3}({\boldsymbol{\zeta}}), \qquad \text{with } {\boldsymbol{\zeta}} \in \hat \Omega ,
\end{equation}
where ${\bf C}^A \in \mathbb{R}^3$ are called \emph{control points}. In a similar
way, it is also possible to define bivariate spline domains in $\mathbb{R}^2$
 or surfaces in $\mathbb{R}^3$, which are commonly used  in CAD (see, e.g.,  \cite{Piegl, cohen2001geometric}).

The control points ${\bf C}^A $  have again the same  important role in
the visualization and interaction with geometries within CAD
softwares. Now, the concept of control polygon generalizes to the control mesh $\mathsf{M}_C$, that is, the mesh connecting the
control points. Figure~\ref{fig:pipe-ctrl-net} shows an example geometry, with its control points and
control mesh. The control mesh defines a
polyhedral domain, denoted $ \Omega_C$, which is an approximation of
$\Omega$: again,  the control domain
$ \Omega_C$ converges to $\Omega$  under $h$-refinement.
\begin{figure}[h!]
\centering
\includegraphics[width=.4\textwidth]{pipe_ctrl_net_deg3.pdf}
\vspace*{-1cm}
\caption{Representation of a geometry (green), with its control points (blue)  and control mesh (black lines) for splines of degree $3$.} \label{fig:pipe-ctrl-net}
\end{figure}
  This is
stated as in the univariate case: we introduce the usual Lagrangian
basis for piecewise trilinear polynomials on the tridimensional
Greville mesh $\mathcal{M}_G$, still denoted  by $\lambda^A(\cdot)$, for the
sake of brevity,
\begin{displaymath}
\lambda^A(\boldsymbol \gamma^{A'})=\left \{
  \begin{aligned}
    1&\text{ if } A=A',\\
    0&\text{ if } A\neq A'.\\
  \end{aligned}
\right .
\end{displaymath}
The control mesh $\mathsf{M}_C$ is the image of the Greville mesh $\mathcal{M}_G$ through the piecewise trilinear mapping
 ${\bf F}_C:\hat \Omega \rightarrow
 \Omega_C$,
\begin{equation}\label{eq:Fh-Bspline}
 {\bf F}_C({\boldsymbol{\zeta}}) = \sum_{A \in \mathcal{A}_{p_1,p_2 ,p_3}} {\bf C}^A
\lambda^A({\boldsymbol{\zeta}}), \qquad \text{with } {\boldsymbol{\zeta}} \in \hat \Omega,
\end{equation}
 which is a parametrization of  $\Omega_C$, since
\begin{displaymath}
  {\bf F}_C ( \boldsymbol  \gamma^A) =  {\bf C}^A.
\end{displaymath}
When ${\bf F}$ is
 smooth enough, as for \eqref{eq:-control-polygon-convergence}, we have
\begin{equation}
  \label{eq:-control-mesh-convergence}
  \underset{{\boldsymbol{\zeta}} \in \hat \Omega }{\sup} \| {\bf F}({\boldsymbol{\zeta}})  -   {\bf F}_C({\boldsymbol{\zeta}}) \| \simeq h^2.
\end{equation}

The control mesh plays a fundamental role in structural mechanics applications where the unknowns are sought as displacements of control points. In our work, we will show how this interpretation can be used also in other contexts.
\begin{remark}\label{rem:control-points-linear-case}
  When $p_1=p_2=p_3=1$ (and all anchors have multiplicity one) the Greville sites coincide with the anchor representations, i.e.,
  $ \boldsymbol \gamma^A = {\boldsymbol{\zeta}}^A$, and $   {\bf F}({\boldsymbol{\zeta}}) = {\bf F}_C ({\boldsymbol{\zeta}})$, $\forall {\boldsymbol{\zeta}} \in   \hat \Omega $, that is,  $ \Omega_C$ and  $\Omega$
  coincide.
\end{remark}

 In CAD and isogeometric analysis the geometry is often
parametrized by Non Uniform Rational B-splines (NURBS).
NURBS are generated  from projective
transformations of splines (see \cite{Piegl}).   A trivariate single-patch NURBS
parametrization of the domain $\Omega \subset\mathbb{R}^3$ is a function  ${\bf F}:\hat  \Omega \rightarrow
\Omega $ defined   as quotient of  linear combination of  B-splines,
\begin{equation}\label{eq:F}
{\bf F}({\boldsymbol{\zeta}}) = \frac{\sum_{A \in \mathcal{A}_{p_1,p_2 ,p_3}}   {\bf C}^A w^A B^A_{p_1,p_2 ,p_3}({\boldsymbol{\zeta}})}{\sum_{A'\in \mathcal{A}_{p_1,p_2 ,p_3}} w^{A'} B^{A'}_{p_1,p_2 ,p_3}({\boldsymbol{\zeta}})}, \qquad {\boldsymbol{\zeta}} \in \hat\Omega,
\end{equation}
where $ {\bf C}^A $ are the NURBS control points  and  $w^A$ the positive
NURBS  weights.

In order to enhance flexibility and allow for more complex  geometries, the  definition of tensor-product spline and NURBS
parametrized domain can be generalized to domains that are union of
$N$ images of  cubes; precisely
\begin{equation}
  \label{eq:F-multi-patch}
  \text{closure } (\Omega )= \bigcup_{k=1,\ldots,N} \text{closure } (\Omega_k)
\end{equation}
where the $ \Omega_k = {\bf F}_k (\hat \Omega)$ are referred to as \emph{patches}, and are assumed to be disjoint. Each patch has its own
parametrization $ {\bf F}_k $, defined on its own spline or NURBS space. The
whole $\Omega$ is then referred  to as a \emph{multi-patch} domain.
For the construction of discrete fields on  a \emph{multi-patch}
domain $\Omega$ we will introduce in Section
\ref{sec:multi-patch-complex-continuity} suitable  \emph{conformity} assumptions. These will
restrict the framework to configurations where it is easy to implement
the proper continuity of the fields at the patches interface.

In this paper,  $\Omega$ is assumed to be  parametrized  either  by
spline or NURBS functions but the unknown fields are always
constructed by splines. This means that, in case of NURBS
geometries,  we  leave the isoparametric concept which is a
fundamental assumption for isogeometric methods in the context of continuum mechanics (see \cite{IGA-book}).

\section{The spline complex}
\label{sec:spline-complex}
This section is devoted to present the \mbox{spline} spaces that are
compatible with the  De~Rham complex. The definition of the spaces is
taken from \cite{BRSV11}, and is given in three dimensions (though the
same construction is generalizable  to arbitrary dimension). We first
recall the construction on the parametric domain $\widehat{\Omega}$, and then
the discrete spaces  on the physical domain $\Omega$ are obtained
 by the {\it push-forward} mapping associated to \eqref{def:iota3}. As shown
in \cite{BRSV11}, it is
also possible to complement these spaces with commuting and continuous
projectors, in the setting of the so called Finite Element Exterior
Calculus (see \cite{AFW06}), however  this issue is not discussed
here. Instead, we discuss  the selection of a suitable
basis for the implementation of the proposed spaces, and  the meaning
of the associated degrees-of-freedom. We will see that the proposed spline spaces
are a natural high-order extension of classical low-order N\'ed\'elec hexahedral
finite elements of the first family (see \cite{NED80}),  obtained in this setting for degree $p_1=p_2=p_3=1$, and that in a natural way they  are related to
cochain or chain complexes of the mesh where they are defined.

\subsection{Complex on the parametric domain $\hat \Omega$}
\label{sec:compatible -splines-parametric-domain}

We recall the following property of univariate splines, from Section
\ref{sect:knot_vector}: the derivative  of a (continuous) spline is a
spline, and in particular
\begin{equation}
  \label{eq:1D-diagram}
        \begin{CD}
{S_{p}(\Xi)} @ > \frac{d}{d\zeta}>>  {S_{p-1}(\Xi')},
\end{CD}
\end{equation}
where $\Xi'$ is the $(p-1)$-open knot vector that coincides with
 the $p$-open knot vector  $\Xi$ except for the boundary knot
repetitions. Moreover, we have that the
derivative of the B-spline associated to an anchor $A$ in
$\mathcal{A}_p(\Xi)$ is a linear combination of the B-splines
associated to the previous and next  adjacent anchors $A^-$ and $A^+$ in $\mathcal{A}_{p-1}(\Xi')$ (only one
adjacent anchor for the first and last $A \in \mathcal{A}_p(\Xi)$); precisely,
denoting by $\Xi_p^A $ the local knot vector (formed by  $p+2$ knots) of $A$
and  by  $ \Xi_{p-1}^{A^\pm} $ the local knot vectors  (formed by  $p+1$
knots) of $A^\pm$, and  with the general  notation of  Section
\ref{sect:knot_vector}, we have
\begin{equation}  \label{eq:derivative-of-splines}
\frac{d}{d\zeta}  N[\Xi_p^A](\cdot) = \frac{p}{|\Xi_{p-1}^{A^-} |}
N[\Xi_{p-1}^{A^-}](\cdot)   - \frac{p}{|\Xi_{p-1}^{A^+} |} N[\Xi_{p-1}^{A^+}](\cdot) ,
\end{equation}
where $ |\Xi_{p-1}^{A^\pm} |$ are the length of the support of the
$(p-1)$-degree B-splines $
N[\Xi_{p-1}^{A^\pm}]$, that is, the difference of the last and first
knots in the local knot vectors $ \Xi_{p-1}^{A^\pm} $. An example is given in Figure~\ref{fig:derivative}. When $A$ is the first
(resp., last) anchor, \eqref{eq:derivative-of-splines} holds with
$N[\Xi_{p-1}^{A^-}]=0 $ (resp., $N[\Xi_{p-1}^{A^+}]=0$). This is a well known
property of B-splines (see \cite{Piegl,DeBoor}) that also suggests the
following scaling of the basis functions of $S_{p}(\Xi)$ and  $S_{p-1}(\Xi')$

\begin{equation}
  \label{eq:spline-basis}
  S_{p}(\Xi) = \text{span} \left \{ B^A_p(\cdot) \equiv
    N[\Xi_p^A](\cdot):\, A \in \mathcal{A}_p(\Xi) \right \} ,
\end{equation}
\begin{equation}
  \label{eq:spline-derivative-basis}
  S_{p-1}(\Xi') = \text{span}\left \{ D^A_{p-1}(\cdot)
    =\frac{p}{|\Xi_{p-1}^{A}|} N[\Xi_{p-1}^A](\cdot):\, A \in \mathcal{A}_{p-1}(\Xi')
  \right \} .
\end{equation}
The  scaling in \eqref{eq:spline-derivative-basis} gives the Curry-Schoenberg B-splines (see
\cite[Ch.~IX]{DeBoor}), that have been already used in isogeometric
analysis in \cite{ratnani}.
Indeed, with the bases \eqref{eq:spline-basis} and
\eqref{eq:spline-derivative-basis}, the matrix associated to the
operator $ \frac{d}{d\zeta}$ is the \emph{edge-vertex incidence} matrix
related to the mesh $\mathcal{M}$,
when $p$ is odd, or the \emph{vertex-edge incidence} matrix related to $\mathcal{M}$,
when $p\geq 2$ is even. We recall that $\mathcal{M}$ also contains zero length
edges and repeated vertices.

\begin{figure}[h!]
\centering
\includegraphics[width=.5\textwidth]{univariate_deriv.pdf}
g\caption{Derivative of the spline associated to the anchor $A$ as a linear combination of the splines associated to $A^-$ and $A^+$.} \label{fig:derivative}
\end{figure}

The observations above are the key ingredients of the trivariate
construction. Following \cite{BRSV11}, and using the notation of
Section \ref{sect:bivariate},  we introduce the following discrete
\color{black}  spaces \color{black} on the parametric domain $\hat \Omega$

\begin{equation}
    \label{eq:discrete-forms-on-omegahat}
  \begin{aligned}
\hat{X}^0_h &:= S_{p_1,p_2,p_3} (\Xi_1, \Xi_2, \Xi_3) , \\
\hat{X}^1_h &:= S_{p_1-1,p_2,p_3} (\Xi'_1, \Xi_2, \Xi_3) \times S_{p_1,p_2-1,p_3} (\Xi_1, \Xi'_2, \Xi_3)\times S_{p_1,p_2,p_3-1} (\Xi_1, \Xi_2, \Xi'_3) ,  \\
\hat{X}^2_h &:= S_{p_1,p_2-1,p_3-1} (\Xi_1, \Xi'_2, \Xi'_3) \times S_{p_1-1,p_2,p_3-1} (\Xi'_1, \Xi_2, \Xi'_3)\times S_{p_1-1,p_2-1,p_3} (\Xi'_1, \Xi'_2, \Xi_3) , \\
\hat{X}^3_h &:= S_{p_1-1,p_2-1,p_3-1} (\Xi'_1, \Xi'_2, \Xi'_3)  .
\end{aligned}
\end{equation}
From \eqref{eq:1D-diagram},  they form a De Rham complex:
\begin{equation}\label{eq:3D-parametric-diagram}
\begin{CD}
 \mathbb{R} @>>>\hat{X}^0_h @>\hat \,{\bf grad}\,>> \hat{X}^1_h @>\hat \,{\bf curl}\,>> \hat{X}^2_h @>\hat \,{\rm div}\,>> \hat{X}^3_h@>>> 0.
\end{CD}
\end{equation}

Moreover, we have the following result.
\begin{theorem}\label{th:exact}
  The sequence \eqref{eq:3D-parametric-diagram} is exact.
\end{theorem}
\begin{proof}
This result has been already presented in \cite{BRSV11}. We present an alternative proof, that will be useful Section~\ref{sect:T-compatible-parametric}.

We have to show that in  \eqref{eq:3D-parametric-diagram} it holds
\begin{equation}
  \label{eq:exactenss-splines-1}
\mathbb{R} =  \mathrm{ker}( \hat \,{\bf grad}\,),
\end{equation}
\begin{equation}
  \label{eq:exactenss-splines-2}
  \mathrm{im}(\hat \,{\bf grad}\,) =  \mathrm{ker}( \hat\,{\bf curl}\,),
\end{equation}
\begin{equation}
  \label{eq:exactenss-splines-3}
  \mathrm{im}(\hat\,{\bf curl}\,)=   \mathrm{ker}( \hat\,{\rm div}\,),
\end{equation}
\begin{equation}
  \label{eq:exactenss-splines-4}
  \mathrm{im}(\hat\,{\rm div}\,)=   \hat{X}^3_h.
\end{equation}
In particular, we have to prove the inclusion $\supseteq$  in
\eqref{eq:exactenss-splines-1}--\eqref{eq:exactenss-splines-4}, since
the other inclusion $\subseteq$ is trivial in all cases. It is also trivial that
\begin{displaymath}
  \mathbb{R} \supseteq \mathrm{ker}( \hat \,{\bf grad}\,).
\end{displaymath}

Let  $\hat {\bf u} = (\hat u_1,\hat u_2,\hat u_3)\in\hat{X}^1_h $, then define
\begin{equation}
  \label{eq:phi-h-poincare}
  \hat \phi(\zeta_1,\zeta_2,\zeta_3) = \int _0^{\zeta_1}  \hat
  u_1(\eta,0,0)\, d \eta + \int _0^{\zeta_2}  \hat
  u_2(\zeta_1,\eta,0)\, d \eta  + \int _0^{\zeta_3}  \hat
  u_3(\zeta_1,\zeta_2,\eta)\, d \eta;
\end{equation}
it is easy to check that  $ \hat {\bf u} = \hat \,{\bf grad}\, \hat
\phi $ when  $\hat \,{\bf curl}\,  {\bf u} =
\boldsymbol 0  $, and that   $\hat\phi \in \hat{X}^0_h$; then
\begin{displaymath}
   \mathrm{im}(\hat \,{\bf grad}\,)\supseteq  \mathrm{ker}( \hat\,{\bf curl}\,).
\end{displaymath}

Consider $\hat \varphi \in\hat{X}^3_h $, and define  $\hat {\bf v} = (\hat
v_1,0,0)\in\hat{X}^1_h$  such that
\begin{equation}
  \label{eq:phi-h-poincare2}
  \hat v_1(\zeta_1,\zeta_2,\zeta_3) = \int _0^{\zeta_1}  \hat
  \varphi(\eta,\zeta_2,\zeta_3)\, d \eta,
\end{equation}
as before, it is easy to check that  $ \hat \varphi = \hat \,{\rm div}\, \hat
{\bf v} $ and that   $\hat{\bf v} \in \hat{X}^2_h$; then
\begin{displaymath}
  \mathrm{im}(\hat\,{\rm div}\,)\supseteq   \hat{X}^3_h.
\end{displaymath}

In order to complete the proof we need to show that
\begin{displaymath}
    \mathrm{im}(\hat\,{\bf curl}\,)\supseteq   \mathrm{ker}( \hat\,{\rm div}\,),
\end{displaymath}
which is implied by
\begin{equation}
  \label{eq:temp2}
    \mathrm{dim}  (\mathrm{im} (\hat\,{\bf curl}\,)) =  \mathrm{dim}(\mathrm{ker}( \hat\,{\rm div}\,)).
\end{equation}
To count dimensions recall from Section \ref{sect:knot_vector} that  $
\mathrm{dim}(S_{p_\ell}(\Xi_\ell))= n_\ell$, $
\mathrm{dim}(S_{p_\ell-1}(\Xi'_\ell))= n_\ell-1$; then from Section
\ref{sect:bivariate} and from  \eqref{eq:discrete-forms-on-omegahat}
we get
\begin{equation}
    \label{eq:dim-of-discrete-forms-on-omegahat}
  \begin{aligned}
\mathrm{dim}(\hat{X}^0_h)& = n_1 n_2 n_3 , \\
\mathrm{dim}(\hat{X}^1_h)&= (n_1 -1 )  n_2 n_3 +  n_1 ( n_2-1) n_3 + n_1 n_2 (n_3-1),  \\
\mathrm{dim}(\hat{X}^2_h) &= n_1 (n_2-1)  (n_3-1) + (n_1 -1 ) n_2 (n_3-1)+(n_1 -1 )( n_2-1)  n_3, \\
\mathrm{dim}(\hat{X}^3_h) &= (n_1 -1 )( n_2-1)(n_3-1) .
\end{aligned}
\end{equation}
Then, by \eqref{eq:exactenss-splines-1}--\eqref{eq:exactenss-splines-2},
\begin{equation*}
  \label{eq:dimension-count-1}
  \begin{aligned}
    \mathrm{dim}  (\mathrm{im} (\hat\,{\bf curl}\,))  &  =  \mathrm{dim}
    (\hat{X}^1_h) -   \mathrm{dim}  (\mathrm{ker} (\hat\,{\bf curl}\,)) \\
  &  =  \mathrm{dim}  (\hat{X}^1_h) -   \mathrm{dim}  (\mathrm{im} (\hat\,{\bf grad}\,))\\
&  =  \mathrm{dim}  (\hat{X}^1_h) -   \mathrm{dim}  (\hat{X}^0_h) +
\mathrm{dim}  (\mathbb{R})\\
& = 2 n_1 n_2 n_3 - n_2 n_3 - n_1 n_3  -  n_1 n_2 + 1
  \end{aligned}
\end{equation*}
and by \eqref{eq:exactenss-splines-4}
\begin{equation*}
  \label{eq:dimension-count-2}
  \begin{aligned}
     \mathrm{dim}(\mathrm{ker}( \hat\,{\rm div}\,))  &  =  \mathrm{dim}
    (\hat{X}^2_h) -   \mathrm{dim}  (\mathrm{im} (\hat\,{\rm div}\,)) \\  &  =  \mathrm{dim}
    (\hat{X}^2_h) -   \mathrm{dim}  (\hat{X}^3_h)  \\
& = 2 (n_1-1) (n_2 -1)(n_3 -1) +   (n_2 -1)(n_3 -1) + (n_1-1) (n_3 -1)  +  (n_1-1)  (n_2 -1)\\& = 2 n_1 n_2 n_3 - n_2 n_3 - n_1 n_3  -  n_1 n_2 + 1,
  \end{aligned}
\end{equation*}
which gives \eqref{eq:temp2}, and as a consequence \eqref{eq:exactenss-splines-3}.

\end{proof}

We now show how suitable basis functions for the spaces can be constructed and as well associated
to geometric entities of the mesh $\mathcal{M}$ by using the concept of anchors.  We focus on basis functions first, and inspired by
\eqref{eq:spline-basis}--\eqref{eq:spline-derivative-basis}
we define them as follows: \color{black}

\begin{equation}
  \label{eq:0-forms-basis}
  \hat{X}^0_h = \text{span} \left \{ \boldsymbol \zeta \mapsto
    B^{A_1}_{p_1}(\zeta_1)  B^{A_2}_{p_2}(\zeta_2)
    B^{A_3}_{p_3}(\zeta_3):\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1,p_2,p_3}(\Xi_1,\Xi_2,\Xi_3) \right \} ,
\end{equation}
\begin{equation}
  \label{eq:1-forms-basis}
  \begin{aligned}
     \hat{X}^1_h & = \text{span }  I\cup II \cup III, \text{ with}\\
 I & = \left \{ \boldsymbol \zeta \mapsto  D^{A_1}_{p_1-1}(\zeta_1)  B^{A_2}_{p_2}(\zeta_2)
    B^{A_3}_{p_3}(\zeta_3)  \hat {\bf e}_1:\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1-1,p_2,p_3}(\Xi'_1,\Xi_2,\Xi_3) \right \} ,\\
 II & = \left \{ \boldsymbol \zeta \mapsto  B^{A_1}_{p_1}(\zeta_1)  D^{A_2}_{p_2-1}(\zeta_2)
    B^{A_3}_{p_3}(\zeta_3) \hat{\bf e}_2:\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1,p_2-1,p_3}(\Xi_1,\Xi'_2,\Xi_3) \right \} ,\\
 III& = \left \{ \boldsymbol \zeta \mapsto
B^{A_1}_{p_1}(\zeta_1)  B^{A_2}_{p_2}(\zeta_2)
    D^{A_3}_{p_3-1}(\zeta_3)  \hat{\bf e}_3:\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1,p_2,p_3-1}(\Xi_1,\Xi_2,\Xi'_3) \right \} ,
  \end{aligned}
\end{equation}

 \begin{equation}
  \label{eq:2-forms-basis}
  \begin{aligned}
     \hat{X}^2_h & = \text{span }  I \cup II \cup III , \text{ with}\\
 I  & = \left \{ \boldsymbol \zeta \mapsto  B^{A_1}_{p_1}(\zeta_1)  D^{A_2}_{p_2-1}(\zeta_2)
   D^{A_3}_{p_3-1}(\zeta_3) \hat {\bf e}_1 :\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1,p_2-1,p_3-1}(\Xi_1,\Xi'_2,\Xi'_3) \right \} ,\\
 II & = \left \{ \boldsymbol \zeta \mapsto D^{A_1}_{p_1-1}(\zeta_1)  B^{A_2}_{p_2}(\zeta_2)
    D^{A_3}_{p_3-1}(\zeta_3)  \hat {\bf e}_2   :\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1-1,p_2,p_3-1}(\Xi'_1,\Xi_2,\Xi'_3) \right \} ,\\
 III & = \left \{ \boldsymbol \zeta \mapsto D^{A_1}_{p_1-1}(\zeta_1)  D^{A_2}_{p_2-1}(\zeta_2)
    B^{A_3}_{p_3}(\zeta_3) \hat {\bf e}_3  :\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1-1,p_2-1,p_3}(\Xi'_1,\Xi'_2,\Xi_3) \right \} ,
  \end{aligned}
\end{equation}
\begin{equation}
  \label{eq:3-forms-basis}
  \hat{X}^3_h = \text{span} \left \{ \boldsymbol \zeta \mapsto
    D^{A_1}_{p_1-1}(\zeta_1)  D^{A_2}_{p_2-1}(\zeta_2)
    D^{A_3}_{p_3-1}(\zeta_3):\, A= (A_1, A_2, A_3) \in
\mathcal{A}_{p_1-1,p_2-1,p_3-1}(\Xi'_1,\Xi'_2,\Xi'_3) \right \} ,
\end{equation}
where $\{\hat {\bf e}_\ell\}_{\ell=1,2,3}$ denote the canonical basis of
$\mathbb{R}^3$. We remark that all basis functions defined in
\eqref{eq:0-forms-basis}-\eqref{eq:3-forms-basis} are non-negative.

 We discuss now the association of the anchors of the bases
\eqref{eq:0-forms-basis}--\eqref{eq:3-forms-basis} to the mesh $\mathcal{M}$
that is associated to $\hat{X}^0_h$, that is, obtained from the knot
vectors $ \Xi_1,\Xi_2,\Xi_3$.
We focus on the relevant case $p=p_1=p_2=p_3$ and consider two possible choices: $p$ is odd, or $p$ is even. \color{black}

  When $p$ is odd, as an immediate consequence of the definition of anchors in one space dimension, we have that:
\begin{itemize}
\item anchors associated with $\hat{X}^0_h$ are $\mathcal{A}_{p,p,p}(\Xi_1,\Xi_2,\Xi_3)$, which are located at the vertices of $\mathcal{M}$, i.e., there is one degree of freedom per vertex;
\item anchors associated with $\hat{X}^1_h$ are located at edges of $\mathcal{M}$
  and there is one degree of freedom per edge. Indeed, e.g., anchors
  associated with the first component of the space $\hat{X}^1_h$, which is
  $S_{p-1,p,p} (\Xi'_1, \Xi_2, \Xi_3)$, are
  $\mathcal{A}_{p-1,p,p}(\Xi'_1,\Xi_2,\Xi_3)$ and are located at the edges
  oriented as $\hat {\bf e}_1$. This means that to each edge of the mesh a is associated a basis function tangential to the edge.
\item anchors associated with $\hat{X}^2_h$ are located at faces and there
  is one anchor per face. More in detail, if we consider the first
  component of $\hat{X}^2_h$, which is  $S_{p,p-1,p-1} (\Xi_1, \Xi'_2,
  \Xi'_3)$, the corresponding anchors  are
  $\mathcal{A}_{p,p-1,p-1}(\Xi_1,\Xi'_2,\Xi'_3)$  and are located at the
  barycenter of the faces $f$ such that $f$ is orthogonal to
  $\hat {\bf e}_1$.  This means that a basis functions normal to the
  face is associated to the face.
\item anchors associated with $\hat{X}^3_h$ are $\mathcal{A}_{p-1,p-1,p-1}(\Xi'_1,\Xi'_2,\Xi'_3)$  and located at barycentres of all elements of $\mathcal{M}$.
\end{itemize}

We now turn to the case when of even degree $p\geq 2$,
$p_1=p_2=p_3=p$.  Note that, according to our definition, and as explained in Section~\ref{sect:knot_vector}, the
meshes corresponding to  the spaces $\hat{X}^1_h$, $\hat{X}^2_h$ and $\hat{X}^3_h$
differ from $\mathcal{M}$ due to the different number of repeated lines at the
boundary.  Instead of working with different meshes for different
spaces, equivalently, we represent in this case too   the anchors of all
spaces on the mesh $\mathcal{M}$ of $\hat{X}^0_h$, keeping into account only the
interior geometrical entities for the representations of anchors of
$\hat{X}^1_h$, $\hat{X}^2_h$ and $\hat{X}^3_h$. \color{black}

 Using the definition of anchors we immediately deduce the following:
\begin{itemize}
\item  anchors associated with $\hat{X}^0_h$ are at the barycentres of all elements in $\mathcal{M}$;
\item anchors associated with $\hat{X}^1_h$ are attached to the barycentres
  of \emph{internal} faces of $\mathcal{M}$ and the corresponding vector basis function is normal to the face;
\item anchors associated with $\hat{X}^2_h$ are attached to
  \emph{internal} edges of $\mathcal{M}$ and the corresponding vector basis function is tangent to its corresponding edge;
\item anchors associated with $\hat{X}^3_h$ are attached to
  \emph{internal}  vertices of $\mathcal{M}$.
\end{itemize}

Clearly,  the positivity of the bases  induces an orientation of edges
and faces of the mesh $\mathcal{M}$. \color{black}

 With the bases \eqref{eq:0-forms-basis}--\eqref{eq:3-forms-basis}, the discrete differential operators in
\eqref{eq:3D-parametric-diagram} are represented by incidence matrices
for the corresponding geometrical entities. If  $p \equiv
p_1 = p_2 = p_3$ is odd, then the operator $\hat \,{\bf grad}\,$ is represented by the
\emph{edge-vertex} incidence matrix of $\mathcal{M}$ and when $p\geq 2$ is
even, by  the \emph{face-element} incidence  matrix of $\mathcal{M}$.
 We observe that, unlike in compatible finite elements, the matrices
representing the differential operators in the selected bases
\eqref{eq:0-forms-basis}--\eqref{eq:3-forms-basis} are
essentially independent of the degree.

The fundamental consequence of the  observations above is stated in
the following proposition.

\begin{proposition}
  \label{p:cochain-and-chain}
The following holds:
  \begin{itemize}
\item The spline complex \eqref{eq:3D-parametric-diagram} for odd degree $p$ is isomorphic to the cochain complex associated with the partition $\mathcal{M}$.
 \item The spline complex \eqref{eq:3D-parametric-diagram} for even
   degree $p$ is isomorphic to the chain complex associated with the
   partition $\mathcal{M}$  \color{black}  without  its boundary, that is, when only the
   interior geometrical entities (faces,
   edges and vertices) are taken    into account, as seen above. \color{black}
\end{itemize}
\end{proposition}
\color{black}

As a matter of fact, this observation, together with the structure of
the matrix representation of differential operators, makes the geometry
of the spline complex for odd degree $p$ very similar, if not equal,
to the one of the finite element complex of lowest order. However  the
spline complex for $p\geq 1$ delivers an approximation which is far
superior than the one of low order finite element.

For $p$ even  we have instead  a chain complex without explicitly
constructing the dual mesh, which has no analogue in the finite
element framework.

Moreover, the use of
anchors and the structure of the mesh at the boundary  guarantee  that in both the chain and cochain complex the
boundary is treated in a simple and canonical way. In the case of
finite elements this is not case (see e.g. \cite{BuCh07}, \cite{GeWiCl04,ClThWe99}, \cite{MR1872728}
or \cite{MR2143847} ) and, moreover,  these features can hardly be obtained in
conjunction with high-order finite element techniques.
Discretization methods based on the use of both chain and cochain complexes in the framework of isogeometric methods are very promising and object of on-going research.\color{black}

We conclude this section with a remark on boundary conditions.
 Consider the case when homogeneous boundary conditions are
 imposed on the whole boundary $\partial \hat {\Omega}$,  leading  to the definition of the discrete spaces  $\hat{X}^0_{0,h} := \hat{X}^0_h \cap
  \,{H}^1_0(\Omegaref)\,$, $\hat{X}^1_{0,h} := \hat{X}^1_h \cap \,{\bf H}_0(\mathbf{curl};{\Omegaref})\,$, $\hat{X}^2_{0,h} := \hat{X}^2_h \cap
  \,{\bf H}_0({\rm div};{\Omegaref})\,$ and $\hat{X}^3_{0,h} :=  \hat{X}^3_h $. These spaces are constructed
  as usual, by removing the functions with non-null trace at the boundary, because univariate  B-spline functions are interpolatory at
  the boundary, as we have discussed in Remark
\ref{rem:interpolatory-B-splines}. The associated De Rham
  complex
\begin{equation}\label{eq:3D-parametric-diagram-bc}
\begin{CD}
0 @>>>\hat{X}^0_{0,h} @>\hat \,{\bf grad}\,>> \hat{X}^1_{0,h} @>\hat \,{\bf curl}\,>> \hat{X}^2_{0,h} @>\hat \,{\rm div}\,>> \hat{X}^3_{0,h}@> \int>> \mathbb{R}
\end{CD}
\end{equation}
is exact, as easily follows by a variation of the argument of
Theorem~\ref{th:exact}.
 The same holds in more general cases, for example  when
the boundary conditions are imposed on a part of  $\partial \hat
{\Omega}$ formed by the union of some faces of the cube $\widehat{\Omega}$.
Since boundary conditions do not represent a conceptual difficulty, in order to keep the presentation as clear as
possible often in our presentation we will not take them into the
framework.

\subsection{Push-forward to the single-patch physical  domain $\Omega$}
\label{sec:compatible -splines-phys-domain}

Following Section \ref{sec:spline-geometries}, we suppose
that the domain $\Omega$ is obtained from $\hat \Omega$ through a
spline or NURBS single-patch mapping ${\bf F}$. Clearly, we need to choose the space for ${\bf F}$.
\begin{assumption}[Isogeometric mapping]\label{ass:F}
  We assume that ${\bf F}$ is either a spline function in
  $[\hat{X}^0_h]^3$,
  or ${\bf F}$ is  a NURBS function as in \eqref{eq:F}, with numerator in
$[\hat{X}^0_h]^3$ and weight denominator in   $\hat{X}^0_h$.
\end{assumption}
Assumption \ref{ass:F} is indeed very natural  in the context of
isogeometric methods: it means that the discrete fields are
constructed from the geometry  knot vectors and bases, possibly after
refinement.

We denote by $\mathsf{M}$ the image of $\mathcal{M}$ through the mapping ${\bf F}$.
$\mathsf{M}$ is then a partition of the physical domain $\Omega$, similar to
the finite element mesh, even though  it  contains elements of
zero area due to knot  multiplicity.  \color{black}
\color{black}

The discrete spaces $X^0_h, \ldots ,
X^3_h $ on the physical domain $\Omega$ can be  defined from the
spaces \eqref{eq:discrete-forms-on-omegahat}  on the parametric
domain $\widehat{\Omega}$ by {\it push-forward}, that is, the inverse of the
transformations  defined in \eqref{def:iota3}, that commute with the
differential operators (as given by the diagrams \eqref{eq:diag1} and
\eqref{eq:diag2}):
\begin{equation} \label{eq:diag-discrete}
\begin{CD}
\mathbb{R} @>>> \hat{X}^0_h @>\hat \,{\bf grad}\,>> \hat{X}^1_h @>\hat \,{\bf curl}\,>> \hat{X}^2_h @>\hat \,{\rm div}\,>> \hat{X}^3_h @>>> 0 \\
@. @A\iota^0AA @A\iota^1AA @A\iota^2AA @A\iota^3AA \\
\mathbb{R} @>>> X^0_h @>\,{\bf grad}\,>> X^1_h @>\,{\bf curl}\,>> X^2_h @>\,{\rm div}\,>> X^3_h @>>> 0, \\
\end{CD}
\end{equation}
 that is, the discrete spaces in the physical domain are defined as
\begin{equation}\label{eq:push-forwarded-discrete-diff-forms}
\begin{aligned}
X^0_h := \{ \phi : \iota^0(\phi) \in \hat{X}^0_h \} , \\
X^1_h := \{ \bu : \iota^1(\bu) \in \hat{X}^1_h \}, \\
X^2_h := \{ \bv : \iota^2(\bv) \in \hat{X}^2_h \}, \\
X^3_h := \{ \varphi : \iota^3(\varphi) \in \hat{X}^3_h \} .
\end{aligned}
\end{equation}
We remark that the space $X^1_h$, which is a discretization of $\,{\bf H}({\mathbf{curl}};\Omega)\,$, is defined through the curl conserving transformation $\iota^1$, and that the space $X^2_h$, which is a discretization of $\,{\bf H}({\rm div};\Omega)\,$, is defined through the divergence conforming transformation $\iota^2$. These are equivalent to the curl and divergence preserving
transformations that are used to define edge and face elements,
respectively (see \cite[Sect. 3.9]{Monk}).

Thanks to the properties of the operators  \eqref{def:iota3} the
push-forwarded spaces $X^0_h, \ldots ,
X^3_h$  inherit the same fundamental properties of $\hat{X}^0_h, \ldots ,
\hat{X}^3_h$, that we have discussed in the previous section:
\begin{itemize}
\item they form an exact
De Rham complex without boundary conditions
\begin{equation}\label{eq:3D-physical-diagram}
\begin{CD}
\mathbb{R}  @>>>X^0_h @>  \,{\bf grad}\,>> X^1_h @>  \,{\bf curl}\,>> X^2_h @>  \,{\rm div}\,>> X^3_h@>>> 0,
\end{CD}
\end{equation}
 or with boundary conditions
\begin{equation}\label{eq:3D-physical-diagram-bc}
\begin{CD}
0 @>>>X^0_{0,h} @>  \,{\bf grad}\,>> X^1_{0,h} @>  \,{\bf curl}\,>> X^2_{0,h} @>  \,{\rm div}\,>> X^3_{0,h}@> \int>> \mathbb{R}.
\end{CD}
\end{equation}
\item the basis functions for $X^0_h, \ldots ,
X^3_h$ are defined by push-forward of the basis functions of $\hat{X}^0_h, \ldots, \hat{X}^3_h$, similarly to
\eqref{eq:push-forwarded-discrete-diff-forms}, and are in  one-to one relation
with  the images of the anchors through ${\bf F}$. See Figure \ref{fig:pipe_ancore}.
\begin{figure}[h!]
\centering
\includegraphics[width=.4\textwidth]{pipe_ancore.pdf}
\vspace*{-1cm}
\caption{We show the mesh $\mathsf{M}$ and the image of the anchors related to the space $X^0_h$ on an example geometry.} \label{fig:pipe_ancore}
\end{figure}

\item since \eqref{eq:diag-discrete}, the
matrices associated with the differential operators $\,{\bf grad}\,$, $\,{\bf curl}\,$
and $\,{\rm div}\,$ on $\Omega$ are the same as the matrices of  $\hat \,{\bf grad}\,$, $\hat\,{\bf curl}\,$
and $\hat\,{\rm div}\,$ on $\widehat{\Omega}$, that is, incidence matrices of
the mesh $\mathsf{M}$.
\item when $p$ is odd (even, respectively), the complex ($X^0_h,\ldots,X^3_h$) is isomorphic to the cochain (chain, respectively) complex associated to the partition $\mathsf{M}$.
\end{itemize}

 Finally, the discrete spaces $X^0_h, \ldots ,
X^3_h$ inherit from their pull-back $\hat{X}^0_h, \ldots ,
\hat{X}^3_h $ optimal approximation properties, if the geometrical
mapping ${\bf F}$ satisfies Assumption~\ref{ass:F} and its inverse is smooth enough (see \cite{BRSV11} for details).

\subsection{Control fields and  degrees-of-freedom interpretation}

\label{sec:dofs}
\color{black}   In this section, we introduce the concept of \emph{control fields},
thanks to which we give an interpretation of the degrees-of-freedom of
the isogeometric fields defined in Section \ref{sec:compatible
  -splines-parametric-domain}--\ref{sec:compatible
  -splines-phys-domain}. The control fields are for the B-spline isogeometric
fields what  the control mesh is for the the B-spline geometry. \color{black}
 We recall that from the geometry control points  we define ${\bf F}_C$ (see
\eqref{eq:Fh-Bspline}),  the piecewise trilinear
function on the  Greville
mesh $\mathcal{M}_G$. The  image of ${\bf F}_C$  is the so-called
control domain $\Omega_C$. The
so-called control mesh  $\mathsf{M}_C$ (which is a partition of
$\Omega_C$) is the image through ${\bf F}_C$ of the Greville
mesh $\mathcal{M}_G$.

As described in Section \ref{sec:spline-geometries},  the standard
way to manipulate a spline parametrization  ${\bf F}$ is by moving its control
points, that is,  the vertices of the control mesh $\mathsf{M}_C$. The
parametrization ${\bf F}_C$ or, equivalently, the control mesh  $\mathsf{M}_C$, carries the {degrees-of-freedom} for the
geometry.  The distance between the two parametrizations ${\bf F}$ and
${\bf F}_C$  is at most $O(h^2)$, as in \eqref{eq:-control-mesh-convergence}.
\color{black}  We now apply the same rationale  for the complex \color{black}  $(X^0_h,..., X^3_h \color{black})$.
Let us first focus on scalar functions on the parametric domain $\hat
\Omega$,
i.e., on the space $\hat{X}^0_h$. Given a spline
\begin{equation}\label{eq:zero-form-field}
  \hat \phi({\boldsymbol{\zeta}}) = \sum_{A \in \mathcal{A}_{p,p ,p}} c^A
 B^A_{p,p,p}({\boldsymbol{\zeta}}), \qquad \text{with } {\boldsymbol{\zeta}} \in \hat \Omega,
\end{equation}
where $c^A$ are its control variables, we associate the
piecewise trilinear function defined on the mesh $\mathcal{M}_G$:
\begin{equation}\label{eq:trilinear-zero-form-field}
  \hat \phi_C({\boldsymbol{\zeta}}) = \sum_{A \in \mathcal{A}_{p,p ,p}} c^A
  \lambda^A({\boldsymbol{\zeta}}), \qquad \text{with } {\boldsymbol{\zeta}} \in \hat \Omega,
\end{equation}
which  carries the same degrees-of-freedom for $\hat \phi$ and indeed is close to $\hat \phi$ (the distance between the two functions
is at most $O(h^2)$,  analogously to  \eqref{eq:-control-polygon-convergence}).  By this relation,
we can interpret the degrees-of-freedom \color{black}  $c^A$ of  $\hat \phi$ as \color{black}  the values
of $ \hat \phi_C $ at each Greville site in $\mathcal{M}_G$.

It should be noted that, if the values of  these degrees-of-freedom are chosen wisely,
splines deliver approximation error of order  $O(h^{p})$ in
the norm of $H^1(\hat \Omega)$, where $p$ is the degree of the spline,
while the corresponding trilinear function can only provide
approximation errors of order $O(h)$. \color{black}

\smallskip
Let  now the geometry come into play. Using \eqref{eq:push-forwarded-discrete-diff-forms}, we set:
   \begin{equation}
\label{eq:anna1}
\phi \circ {\bf F} = \hat \phi  \text{ and } \phi_C \circ {\bf F}_C = \hat \phi_C.
\end{equation}
 The degrees-of-freedom for $\phi$ are \color{black}   the values
of $ \phi_C $ at the vertices of $\mathsf{M}_C$, that is, at the control
points. Or, we can say that  the field  $\phi_C$ determines, or
\emph{controls}, $\phi$,  and its degrees-of-freedom are the values of
$\phi_C$ at control points. \color{black} In
 Figure \ref{fig:pipe-ctrl-net}, the location of control points (blue
 bullet) is shown on an example geometry.  The
 field $\phi_C$  plays the role of \emph{control field} for $\phi$.
As for the parametric space, there are wise choices of the degrees-of-freedom
which ensure an approximation error of order  $O(h^{p})$ in
the norm of $H^1(\Omega)$, while the corresponding trilinear function can only provide
approximation errors of order $O(h)$.

The same reasoning can be applied to the whole complex
$(X^0_h,..., X^3_h \color{black})$,  defined in  Section
\ref{sec:compatible -splines-parametric-domain} and
\ref{sec:compatible -splines-phys-domain},  from degrees  $ p_\ell$ and
knot vectors $ \Xi_\ell$.  \color{black}  Indeed, we  introduce the
\emph{control complex} $(Z^0_h,..., Z^3_h \color{black})$ which is
obtained, still following Section \ref{sec:compatible -splines-parametric-domain} and
\ref{sec:compatible -splines-phys-domain},  with  the following
choices for \color{black}  $\hat Z^0_h $: \color{black}

\begin{itemize}
\item degrees in all directions equal to $1$;
\item the knot vector in the $\ell$ direction is the ordered
  collection of points $\{ \gamma_\ell^A  \ :\ A\in \mathcal{A}_{p_\ell}
  (\Xi_\ell) \}$, $\ell=1,2,3$,  along with repeated 0 and 1 to make the knot vectors open,
\end{itemize}
and {\color{black} replacing the geometric mapping ${\bf F}$ with ${\bf F}_C$ in the
pullbacks} \eqref{def:iota3}.
The complex $(Z^0_h,...,  Z^3_h \color{black})$ corresponds to the low
order finite element complex defined on the control mesh $\mathsf{M}_C$ and
it is immediate to see that if $\phi$ in \eqref{eq:anna1} belongs to
$X^0_h$, then the corresponding $\phi_C$ belongs to $Z^0_h$.
We denote by $I^0_h: X^0_h \to Z^0_h$ the operator which
associates  $\phi$ to $\phi_C$, \color{black} and in an analogous  way,  we
define the operators  $I^j_h : X^j_h \to Z^j_h$,
$j=0,\ldots,3$. These operators \color{black} are represented by  identity matrices when
the spaces are endowed with the bases described in Section
\ref{sec:compatible -splines-parametric-domain}.
It is not difficult to see that, in view of the structure of the matrices associated to differential operators,  the following diagram commutes:
\begin{equation}
\label{eq:cd-XZ}
\begin{CD}
\mathbb{R} @>>> X^0_h @> \,{\bf grad}\,>> X^1_h @>\,{\bf curl}\,>> X^2_h  @>\,{\rm div}\,>> X^3_h @>>> 0 \\
@. @V{I^0_h}VV @V{I^1_h}VV @V{I^2_h}VV @V{I^3_h}VV \\
\mathbb{R} @>>> Z^0_h @>\,{\bf grad}\,>> Z^1_h @>\,{\bf curl}\,>> Z^2_h @>\,{\rm div}\,>> Z^3_h @>>> 0. \\
\end{CD}
\end{equation}

 Let us comment about the meaning of the diagram
 \eqref{eq:cd-XZ}. First of all,  it says that the geometric structure
 of the spline complex $(X^0_h,..., X^3_h \color{black})$  is the one of
 the  low order finite element complex $(Z^0_h,..., Z^3_h \color{black})$
 on the control mesh.  The discrete fields in  $X^j_h$ can be
 associated to \emph{control fields} in  $Z^j_h$, through the operator
 $I^j_h$, as we have discussed for $j=0$ above. \color{black}  For example, we can
 say that there is a N\'ed\'elec field ${\bf u}_C$ which controls ${\bf u}$
 and the degrees-of-freedom are, in this case,  its circulation on the
 edges of the control mesh $\mathsf{M}_C$.   Moreover,
 following a reasoning similar to the one in Section
 \ref{sec:spline-geometries},  the operators $I^j_h$ are point-wise
 converging to the identity when $h$ goes to zero.
 We stress again that the order of approximation of the complex
$(X^0_h,..., X^3_h \color{black})$ is $O(h^p)$ while the control complex
$(Z^0_h,..., Z^3_h \color{black})$ only exhibits  first order convergence {\color{black} in the norm of $X^i$}.

\begin{figure}[h!]
\centering
\includegraphics[width=.4\textwidth]{pipe_ctrl_net_frecce_deg3.pdf}
\vspace*{-1cm}
\caption{Representation of the degrees-of-freedom location for the space $X^1_h$, on the green geometry for degrees $p_1=p_2=p_3=3$.} \label{fig:pipe-ned}
\end{figure}

 Finally, it should be noted that, as it is well known,  the complex $(Z^0_h,..,Z^3_h)$ is always isomorphic to the cochain complex of  the partition $\mathsf{M}_C$, while for the complex  $(X^0_h,..,X^3_h)$ Proposition \ref{p:cochain-and-chain} holds. This is in accordance with the fact that, when $p_\ell$ are all even, the control mesh $\mathsf{M}_C$ can be interpreted as a partition dual to $\mathsf{M}$, in the sense that the chain of $\mathsf{M}$ is isomorphic to the cochain complex of  $\mathsf{M}_C$.

\subsection{Push-forward to the multi-patch  physical  domain
  $\Omega$}
\label{sec:multi-patch-complex-continuity}
\color{black} In this section we construct the spline complex on a multi-patch geometry by addressing the questions of
how conformity is imposed at the interfaces between patches.\color{black}

We consider now a multi-patch domain $\Omega$ which is parametrized
from a reference patch $\hat \Omega$ through the
spline or NURBS mappings ${\bf F}_k$, $k=1,\ldots,N$, as in
\eqref{eq:F-multi-patch}.  Each patch is endowed with a (possibly
different) spline space and therefore for each  $k=1,\ldots,N$ we
can define discrete spaces $[\hat{X}^0_h]_k, \ldots, [\hat{X}^3_h]_k$ such
that a De Rham complex \eqref{eq:3D-parametric-diagram} holds.
Assuming each  ${\bf F}_k$ verifies Assumption \ref{ass:F}, then, as shown in
Section \ref{sec:compatible -splines-phys-domain} we  push-forward
patch-by-patch the discrete spaces $[\hat{X}^0_h]_k, \ldots, [\hat{X}^3_h]_k$
and obtain, on each $ \Omega_k = {\bf F}_k (\hat \Omega) $, the discrete
spaces $[X^0_h]_k, \ldots, [X^3_h]_k$ that fulfill the De Rham
complex \eqref{sec:compatible -splines-phys-domain} on each patch.

Then, the last and main step is
to assemble the spaces $X_h^j  \subset \bigoplus_{k=1,\ldots,N}
[X_h^j]_k$, and add the relevant continuity conditions  at the
inter-patches boundaries: trace continuity for $X^0_h$, tangential  trace
continuity for $X^1_h$, normal trace continuity for $X^2_h$, no continuity for $X^3_h$.
For this purpose, we introduce a \emph{conformity} condition as,
e.g., in \cite{kleissieti}. This condition guarantees that the
geometry parametrizations of the patches are equivalent at the
patch interfaces and, since Assumption \ref{ass:F}, it can be stated on
the spaces $[\hat{X}^0_h]_k$.

\begin{assumption}[Geometrical conformity]
  \label{ass:IETI}
On each non-empty patch interface $\Gamma= \partial \Omega_{k}
\cap \partial \Omega_{k'} $, with $k \neq k'$,  the spaces
$ [X^0_h]_{k} | _{\Gamma}$ and $ [X^0_h]_{k'} | _{\Gamma}$
coincide, as the corresponding bases do.
\end{assumption}

This means that the meshes $\mathsf{M}_k$ and
$\mathsf{M}_{k'}$ match on $\Gamma$, and therefore \color{black}
\begin{displaymath}
  \mathsf{M} = \bigcup_{k=1,\ldots,N} \mathsf{M}_k
\end{displaymath} \color{black}
is  a locally structured but globally unstructured mesh $\mathsf{M}$  on
$\Omega$. In a similar way, the patch control meshes $[\mathsf{M}_C]_k$ match
conformally and \color{black}
\begin{displaymath}
  \mathsf{M}_C = \bigcup_{k=1,\ldots,N} [\mathsf{M}_C]_k
\end{displaymath} \color{black}
is a locally structured but globally unstructured mesh $\mathsf{M}_C$ of hexahedra
on $\Omega_C$, the union of the patch control domains
$[\Omega_C]_k$.

Assumption \ref{ass:IETI}  corresponds to the \emph{full-matching}
conditions of  \cite{kleissieti}, to which we refer for further
details.

Having conformity we can implement the continuity conditions
easily. Indeed, due to the definitions in Sections~\ref{sec:compatible
  -splines-parametric-domain}--\ref{sec:dofs}, the needed continuity  holds
for $(X^0_h,..., X^3_h)$ if and only if it holds for $(Z^0_h,..., Z^3_h
)$ on the global mesh $\mathsf{M}_C$. Continuity for $Z^0_h, Z^1_h$ and
$Z^2_h$ is imposed by merging the coincident degrees-of-freedom at
the interfaces, which in the case of $Z^1_h$ and $Z^2_h$ also
requires to take into account the orientation; see Figure~\ref{multi-patch}. This is however the
same as in finite elements (indeed, the control fields are classical
finite elements). This merging automatically gives the degrees-of-freedom for fields in
$(X^0_h, ..., X^3_h)$.

\begin{figure}[h]
\begin{subfigure}[Control variables of the two patches before merging.]
{\includegraphics[width=0.48\textwidth]{multipatch1.pdf}} \label{fig:multipatch1}
\end{subfigure}
\begin{subfigure}[Control variables of the merged patch.]
{\includegraphics[width=0.41\textwidth]{multipatch2.pdf}} \label{fig:multipatch2}
\end{subfigure}
\caption{Implementing continuity for $ X^1_h$ on a two-patch domain. The orientation of the edges at the interface is chosen as that of the {\color{black} lower} patch.}
\label{multi-patch}
\end{figure}

\section{Beyond the tensor product structure: T-splines}
\label{sec:TT}

In this section, we generalize the definition of tensor-product
B-splines to  T-splines
\cite{Sederberg_Zheng,Sederberg_Cardon,Bazilervs_Calo_Cottrell_Evans}. The
theory of T-splines is well developed in two dimensions (see the very
recent papers \cite{LZSHS12, BBCS12, BBSV12, Li_Scott})
while it is still incomplete in three dimensions (though some recent
important advances have been recently proposed  in \cite{Zhang2012}). For this reason,  we only present in Sections
\ref{sect:T-mesh} and~\ref{sect:AS}, T-splines in two
dimensions and construct, in Section \ref{sect:T-compatible-parametric}, a
discrete T-spline based complex. The extension to three
dimensions is given in Section~\ref{sect:3D-T-splines} by tensor-product of the
two-dimensional T-spline spaces with B-spline one-dimensional spaces.
\color{black}
\subsection{T-mesh}\label{sect:T-mesh}
Let $n_\ell \in \mathbb N$ and the degree $p_\ell \in \mathbb N$, and let $\Xi_\ell
= \{\xi_{\ell,1}, \ldots, \xi_{\ell,n_\ell+p_\ell+1}\}$ be a
$p_\ell$-open knot vector for $\ell = 1,2$. A  T-mesh is  a
rectangular tiling of the unit square $[0,1]^2$, such that all
\emph{vertices}  belong to $\Xi_1 \times \Xi_2$. A T-mesh may contain
 interior
vertices that connect only three edges,  called \emph{T-junctions},
that break the tensor product structure of the mesh (see
Figure~\ref{fig:Tmesh-p3}).  We will say that a T-junction is horizontal (respectively, vertical) if
the missing edge is horizontal (resp. vertical). By an abuse of notation,
we still denote a T-mesh by $\mathcal{M}$.

As before,  we represent the  knot multiplicities by repeated lines
close to each other, with now  the line multiplicity possibly
varying along lines  (see
\cite[Section~4.3]{Bazilervs_Calo_Cottrell_Evans}). The only exception are the boundary lines, that maintain the same multiplicity all along the line. As in B-spline
meshes, the vertical (resp. horizontal) lines at the boundaries have
multiplicity $\lfloor p_1/2 \rfloor + 1$ (resp. $\lfloor p_2/2 \rfloor + 1$).

\subsection{Analysis suitable T-meshes and  T-splines.}
\label{sect:AS}
We define, for a horizontal (resp. vertical) T-junction $T$, the
$k$-bay face-extension as the horizontal (resp. vertical) closed
segment that  extends from $T$  in the direction of the missing
edge until it intersects $k$ lines of the mesh $\mathcal{M}$. The $k$-bay
edge-extension is defined analogously extending the segment in the
opposite direction.

Following \cite{BBSV12}, we define the extension of a horizontal
(resp. vertical) T-junction $T$  the union of the
$\lfloor(p_1+1)/2\rfloor$-bay face-extension, and the
$\lfloor(p_1-1)/2\rfloor$-bay edge-extension (resp.  the union of the $\lfloor(p_2+1)/2\rfloor$-bay face-extension, and the $\lfloor(p_2-1)/2\rfloor$-bay edge-extension). More precisely, if $p_\ell$ is odd we extend $(p_\ell+1)/2$ bays in the direction of the missing edge, and $(p_\ell-1)/2$ bays in the opposite direction; if $p_\ell$ is even we extend $p_\ell/2$ bays in both directions. An example is given in Figure~\ref{fig:extensions}.

\begin{figure}[h!]
\centering
\includegraphics[width=0.4\linewidth]{extensions.pdf}
\caption{Extensions for degree $p_1 = 2$ (horizontal) and $p_2 = 3$ (vertical). Dashed lines represent the face extensions.} \label{fig:extensions}
 \end{figure}

\begin{definition}\label{def:analysis-suitable}
A T-mesh $\mathcal{M}$ is analysis suitable for degrees $p_1$ and $p_2$,
denoted  $\mathcal{M} \in \mathsf{AS}_{p_1,p_2}$, if vertical extensions and horizontal extensions do not intersect.
\end{definition}

Analysis suitable T-meshes were first identified in \cite{LZSHS12} in
the bi-cubic case, and generalized to arbitrary degree in
\cite{BBSV12}.  Despite their very geometric definition, analysis suitable T-meshes and T-splines enjoy fundamental properties  which make their use in isogeometric analysis really promising. Some of these properties will be discussed in what follows.

As in the case of B-splines,
anchors are inferred from the T-mesh $\mathcal{M}$ and their position depends upon the parity of $p_1$ and $p_2$.
For example, if both $p_1$ and $p_2$ are odd, anchors are at the vertices of the $\mathcal{M}$, if they are even, anchors are at the barycenters of elements and so on  (see \cite{BBSV12}). We will denote the set of
anchors by $\mathcal{A}_{p_1,p_2}(\mathcal{M})$, or simply $\mathcal{A}_{p_1,p_2}$.

T-spline basis functions are constructed as B-splines associated to the anchors
$\mathcal{A}_{p_1,p_2}(\mathcal{M})$, and defined from two \emph{local knot vectors}, $\Xi_1^A = \{\xi_{1,i_1},
\ldots, \xi_{1,i_{p_1+2}} \} \subset \Xi_1$ and $\Xi_2^A = \{\xi_{2,i_1},
\ldots, \xi_{2,i_{p_2+2}} \} \subset \Xi_2$. To construct the horizontal knot vector $\Xi_1^A$ we trace the horizontal line through $A$ and select its intersections with vertical lines of $\mathcal{M}$, depending on the degree $p_1$: if $p_1$ is even we choose the first $(p_1+2)/2$ intersections to the left of $A$, and the first $(p_1+2)/2$ to the right; if $p_1$ is odd we first select the coordinate of the anchor $A$, and then the first $(p_1+1)/2$ intersections to the left and to the right of $A$. In the case that we arrive at the boundary, we add the value 0 or 1 as many times as needed to complete the $p_1+2$ entries of $\Xi_1^A$. The construction of $\Xi_2^A$ is analogous, and depends on $p_2$. Examples are shown in Figure~\ref{fig:local_knot_vector} for $p_1 = 2$ and $p_2 = 3$. For more details we refer to \cite{Bazilervs_Calo_Cottrell_Evans}.

\begin{figure}
\centering
\includegraphics[scale=0.4]{knot_vector23.pdf}
\caption{\color{black} Representation \color{black} of the local knot vectors for degrees $p_1 = 2$ and $p_2 = 3$. The local knot vectors are $\Xi_1^{A^1} = \{0, 0, 1/6, 2/6\}$, $\Xi_2^{A^1} = \{0, 0, 0, 1/6, 2/6\}$, and $\Xi_1^{A^2} = \{3/6, 4/6, 5/6, 1\}$, $\Xi_2^{A^2} = \{0, 2/6, 3/6, 4/6, 5/6\}$.} \label{fig:local_knot_vector}
\end{figure}

The T-spline function associated to the anchor $A$ is denoted as:
\begin{equation}\label{eq:T-spline-function}
B^A_{p_1,p_2}({\boldsymbol{\zeta}}) = N[\Xi_1^A] (\zeta_1) N[\Xi_2^A] (\zeta_2),
\quad {\boldsymbol{\zeta}} = (\zeta_1, \zeta_2) \in (0,1)^2,
\end{equation}
they are linearly independent (see \cite{BBSV12}) and by definition span the
T-spline space $T_{p_1,p_2}= T_{p_1,p_2}(\mathcal{M})$:
\begin{equation}\label{Tsplinespace}
 T_{p_1,p_2}(\mathcal{M}) := \mbox{span} \{B_{p_1,p_2}^{A} \ : \ A \in \mathcal{A}_{p_1,p_2}(\mathcal{M}) \} .
\end{equation}

Definition \ref{def:analysis-suitable} guarantees fundamental
properties of the T-spline space \eqref{Tsplinespace}.  In
\cite{BBCS12,BBSV12} it is defined a dual basis for
the T-spline functions constructed from an analysis suitable T-mesh,
thus proving the linear independence of \eqref{eq:T-spline-function} (see also \cite{LZSHS12}) and good approximation properties
for the  space \eqref{Tsplinespace}.  We also remark that the construction of the
local knot vectors described above is analogous to the one in
  \cite{BBSV12}  for analysis suitable T-meshes.

Finally, we define the \emph{extended T-mesh} of $\mathcal{M}$, and denote it
by $\mathcal{M}_{\rm ext}$, as the T-mesh obtained from $\mathcal{M}$ by adding all the
T-junction extensions.
 The extended T-mesh, sometimes also called B\'ezier mesh, is the
 minimal mesh such that the functions \eqref{eq:T-spline-function} restricted to the
 non-empty elements are bivariate polynomials of degree
 $(p_1,p_2)$. The importance of the extended mesh for implementation
 \cite{Scott2011126}, local refinement \cite{Li_Scott} and approximation
 properties of T-splines \cite{BBSV12} is already known.  In particular,  for the implementation,
and in order to ensure accuracy, integration has to be performed on the elements of $\mathcal{M}_{\rm ext}$ and
this means that the data structure is constructed based on $\mathcal{M}_{\rm ext}$.

Finally, a key result which is useful in the construction of compatible T-spline discretizations,  is the characterization stated in the following
proposition.
\begin{proposition}\label{prop:characterization}
  Given an analysis suitable T-mesh $\mathcal{M}$, if furthermore no T-junction
  extensions of any kind intersect each other or intersect mesh
  lines with multiplicity greater than one,  then the T-spline  space \eqref{Tsplinespace}
  is the space of all piece-wise bivariate polynomials of degree
 $(p_1,p_2)$ on   $\mathcal{M}_{\rm ext}$ with the same continuity of the T-spline
 functions \eqref{eq:T-spline-function} at the mesh lines.
\end{proposition}
\begin{proof}
  The case $p_1=p_2=3$ has been covered
  in \cite{Li_Scott}, while the
  general case is a work in progress by A. Bressan in
  \cite{bressan-in-preparation}.
Related works are also \cite{Mourr10}, and \cite{LiCh11} which show the
mathematical complexity of the problem.
 The condition that the extensions do not intersect lines with
multiplicity greater than one can be removed at the price of a more complex
statement, which we do not consider here for the sake of simplicity. \color{black}
\end{proof}

\subsection{Two-dimensional De Rham complex with T-splines  on the parametric domain $(0,1)^2$} \label{sect:T-compatible-parametric}
The aim of this section is to introduce a two-dimensional  T-spline
based De Rham complex, thus generalizing the tensor-product
construction  of  \color{black} section~\ref{sec:spline-complex}. \color{black} Throughout this
Section we will assume, for the sake of simplicity, $p_1 = p_2 =
p$. The results are also valid in the general case, but the proofs
become more intricate.

As for B-splines, T-splines spaces are
constructed by a suitable selection  of the polinomial  degree in the two
directions and by a suitable design of the mesh, that is, the knot vectors. The main
difference now is that we need to modify the mesh $\mathcal{M}$, depending of the
form degree, not only at the boundary but also around T-junctions.

Let $p \in \mathbb N$, let $\Xi_1,\Xi_2$ be two $p$-open knot
vectors, and let $\mathcal{M} \in \mathsf{AS}_{p,p}$ be a T-mesh with knot repetitions,
as defined in Section~\ref{sect:T-mesh}.
The starting mesh is ${\mathcal M}^{0} \equiv \mathcal{M}$, on which we define the space of scalar fields:
\begin{equation}
    \label{eq:T-splines-0-forms-on-parametric}
 \hat Y^0_h := T_{p,p}({\mathcal M}^{0}).
\end{equation}

 The \color{black}  T-splines  vector fields \color{black}  are defined on the two T-meshes ${\mathcal M}^{1}_1$ and ${\mathcal M}^{1}_2$.
 If $p$ is odd, ${\mathcal M}^{1}_1$ is obtained from $\mathcal{M}$ by adding the first-bay
  face-extension of all  horizontal T-junctions. If $p$ is even,
 ${\mathcal M}^{1}_1$  is  equal to $\mathcal{M}$ everywhere but on the boundary where, due to the definition of $\mathcal{M}$, and recalling Section~\ref{sect:knot_vector},  the first and the last
 vertical columns of elements of $\mathcal{M}$ are removed. We define analogously ${\mathcal M}^{1}_2$,
 reasoning in the vertical direction: if $p$ is odd  ${\mathcal M}^{1}_2$  is
 defined by adding the first-bay face-extension of all the vertical
 T-junctions, and if $p$ is even, it is defined by removing the first
 and last horizontal rows of elements of $\mathcal{M}$. Then the \color{black} vector
 fields \color{black}   and the {\color{black} rotated} \color{black} vector
 fields \color{black}  are defined as

 \begin{equation}
    \label{eq:T-splines-1-forms-on-parametric}
 \hat Y^1_h := T_{p-1,p}({\mathcal M}^{1}_1) \times  T_{p,p-1}({\mathcal M}^{1}_2).
\end{equation}

\begin{equation}
    \label{eq:T-splines-1*-forms-on-parametric}
 \hat Y^{1*}_h := T_{p,p-1}({\mathcal M}^{1}_2) \times  T_{p-1,p}({\mathcal M}^{1}_1).
\end{equation}

Finally,  the \color{black}  last space $\hat Y^{2}_h$ \color{black} is  defined on the T-meshes ${\mathcal M}^{2}$: if $p$ is
odd, ${\mathcal M}^{2}$ is obtained from $\mathcal{M}$ by adding all the first-bay
face-extensions (horizontal and vertical), and if $p$ is even it is
defined by removing the first and last rows and columns of elements in
$\mathcal{M}$.

Then
\begin{equation}
    \label{eq:T-splines-2-forms-on-parametric}
 \hat Y^{2}_h := T_{p-1,p-1}({\mathcal M}^{2}).
\end{equation}

An example of the sequence of meshes is shown in
Figure~\ref{fig:Tmesh-p3} for $p=3$, and in Figure~\ref{fig:Tmesh-p2}
for $p=2$. We notice that whenever ${\mathcal M}^{0} = \mathcal{M}$ is a tensor product
mesh, the construction is equivalent to the one presented in
Section~\ref{sec:compatible -splines-parametric-domain} for
B-splines. Indeed, for odd $p$ the four meshes are equal to $\mathcal{M}$,
because there are no T-junctions, and for even $p$ they only differ in
the number of line repetitions on the boundary.

\begin{figure}[h!]
\begin{center}
\begin{subfigure}[Mesh ${\mathcal M}^{0}$]{\includegraphics[width=.47\textwidth]{tsp_cochain_33.pdf}}
\end{subfigure}
\end{center}
\begin{subfigure}[Mesh ${\mathcal M}^{1}_1$]{\includegraphics[width=.47\textwidth]{tsp_cochain_23.pdf}}
\end{subfigure}
\begin{subfigure}[Mesh ${\mathcal M}^{1}_2$]{\includegraphics[width=.47\textwidth]{tsp_cochain_32.pdf}}
\end{subfigure}
\begin{center}
\begin{subfigure}[Mesh ${\mathcal M}^{2}$]{\includegraphics[width=.47\textwidth]{tsp_cochain_22.pdf}}
\end{subfigure}
\end{center}
\caption{Sequence of meshes for the spline complex, with their respective anchors, for $p=3$.} \label{fig:Tmesh-p3}
\end{figure}

\begin{figure}[h!]
\begin{center}
\begin{subfigure}[Mesh ${\mathcal M}^{0}$]{\includegraphics[width=.47\textwidth]{tsp_chain_22.pdf}}
\end{subfigure}
\end{center}
\begin{subfigure}[Mesh ${\mathcal M}^{1}_1$]{\includegraphics[width=.47\textwidth]{tsp_chain_12.pdf}}
\end{subfigure}
\begin{subfigure}[Mesh ${\mathcal M}^{1}_2$]{\includegraphics[width=.47\textwidth]{tsp_chain_21.pdf}}
\end{subfigure}
\begin{center}
\begin{subfigure}[Mesh ${\mathcal M}^{2}$]{\includegraphics[width=.47\textwidth]{tsp_chain_11.pdf}}
\end{subfigure}
\end{center}
\caption{Sequence of meshes for the spline complex, with their respective anchors, for $p=2$.} \label{fig:Tmesh-p2}
\end{figure}

The choice of these meshes becomes clear when computing the derivatives. For instance, let $p=3$ and consider the simple example of a mesh with only one horizontal T-junction, as in Figure~\ref{fig:oneT0}. Choosing the anchor $A \in \mathcal{A}_{p,p}({\mathcal M}^{0})$ located at the T-junction, it is clear from \eqref{eq:derivative-of-splines} that $\frac{\partial B^A_{p,p}}{\partial x}$ is a linear combination of $B^{A^l}_{p-1,p}$ and $B^{A^r}_{p-1,p}$, with $A^l, A^r \in \mathcal{A}_{p-1,p}({\mathcal M}^{1}_1)$ as in Figure~\ref{fig:oneT11}. Hence, $\frac{\partial B^A_{p,p}}{\partial x} \in T_{p-1,p}({\mathcal M}^{1}_1)$ (see \eqref{eq:T-splines-1-forms-on-parametric}), but since $A^l \not \in \mathcal{A}_{p-1,p}({\mathcal M}^{0})$, we have $\frac{\partial B^A_{p,p}}{\partial x} \not \in T_{p-1,p}({\mathcal M}^{0})$. The argument is analogous for the  partial derivative with respect to the $y$ direction, with vertical T-junctions.

\begin{figure}
\begin{subfigure}[Anchor $A \in \mathcal{A}_{p,p}({\mathcal M}^{0})$]{\includegraphics[width=.47\textwidth]{mesh_oneT.pdf} \label{fig:oneT0}}
\end{subfigure}
\begin{subfigure}[Anchors $A^l, A^r \in \mathcal{A}_{p-1,p}({\mathcal M}^{1}_1)$]{\includegraphics[width=.47\textwidth]{mesh_oneT2.pdf} \label{fig:oneT11}}
\end{subfigure}
\caption{The first partial derivative $\frac{\partial B^A_{p,p}}{\partial x}$ is a linear combination of $B^{A^l}_{p-1,p}$ and $B^{A^r}_{p-1,p}$. \label{fig:oneT}}
\end{figure}

We have the following result.
\begin{proposition}\label{prop:extension}
Assuming    $\mathcal{M}^0 \in \mathsf{AS}_{p,p}$, \color{black}  it holds that ${\mathcal M}^{1}_1 \in \mathsf{AS}_{p-1,p}$, ${\mathcal M}^{1}_2 \in \mathsf{AS}_{p,p-1}$, and ${\mathcal M}^{2} \in \mathsf{AS}_{p-1,p-1}$.
 Moreover,  $(\mathcal{M}^0)_{{\rm ext}} = ({\mathcal M}^{1}_1)_{\rm ext} = ({\mathcal M}^{1}_2)_{{\rm ext}} =
 ({\mathcal M}^{2})_{{\rm ext}}$.

\end{proposition}
\begin{proof}
The result is an immediate consequence of $\mathcal{M} \in \mathsf{AS}_{p,p}$, and the length of the extensions specified in Section~\ref{sect:AS}.
\end{proof}
\begin{remark} Note that, although the four meshes are different, all integral
  computations are carried out in the extended T-mesh, which is the
  same for all the spaces. As a consequence the four spaces can be
  implemented within the same data structure, which is based on one
  single mesh, but with different basis functions for each space. This
  is also what occurs  with  standard finite elements.

\end{remark}
\color{black}

The bases for $ \hat Y^{0}_h, \ldots   \hat Y^{2}_h $, are formed by
T-spline functions \eqref{eq:T-spline-function}  with a scaling as in
Section~\ref{sec:compatible -splines-parametric-domain}. Precisely, introducing the notation  $D[\Xi^A_i](\zeta_i) := \frac{p}{|\Xi^A_i|}N[\Xi^A_i](\zeta_i)$, we have:

\begin{equation}
  \label{eq:0-forms-basis-Tsplines}
  \hat Y^{0}_h = \text{span} \left \{ (\zeta_1, \zeta_2)  \mapsto
   N[\Xi_1^A] (\zeta_1) N[\Xi_2^A] (\zeta_2) :\, A \in \mathcal{A}_{p,p}(\mathcal{M}^0)  \right \} ,
\end{equation}

\begin{equation}
  \label{eq:1-forms-basis-Tspline}
  \begin{aligned}
    \hat Y^{1}_h & = \text{span }  I \cup II, \text{ with}\\
 I & =\left \{ (\zeta_1, \zeta_2)  \mapsto
    D[\Xi_1^A] (\zeta_1) N[\Xi_2^A] (\zeta_2) \hat {\bf e}_1 :\, A \in \mathcal{A}_{p-1,p}(\mathcal{M}^1_1)  \right \} ,\\
 II & =\left \{ (\zeta_1, \zeta_2)  \mapsto
    N[\Xi_1^A] (\zeta_1) D[\Xi_2^A] (\zeta_2) \hat {\bf e}_2 :\, A \in \mathcal{A}_{p,p-1}(\mathcal{M}^1_2)  \right \}  .
  \end{aligned}
\end{equation}

\begin{equation}
  \label{eq:1*-forms-basis-Tspline}
  \begin{aligned}
    \hat Y^{1*}_h & = \text{span }  I \cup II, \text{ with}\\
 I & =\left \{ (\zeta_1, \zeta_2)  \mapsto
    N[\Xi_1^A] (\zeta_1) D[\Xi_2^A] (\zeta_2) \hat {\bf e}_1 :\, A \in \mathcal{A}_{p,p-1}(\mathcal{M}^1_2)  \right \} ,\\
 II & =\left \{ (\zeta_1, \zeta_2)  \mapsto
    D[\Xi_1^A] (\zeta_1) N[\Xi_2^A]
   (\zeta_2) \hat {\bf e}_2 :\, A \in \mathcal{A}_{p-1,p}(\mathcal{M}^1_1)  \right \} .
  \end{aligned}
\end{equation}
\begin{equation}
  \label{eq:2-forms-basis-Tsplines}
  \hat Y^{2}_h = \text{span} \left \{ (\zeta_1, \zeta_2)  \mapsto
    D[\Xi_1^A] (\zeta_1) D[\Xi_2^A] (\zeta_2) :\, A \in \mathcal{A}_{p-1,p-1}(\mathcal{M}^2)  \right \} ,
\end{equation}

The main result of this section is the following.
\begin{theorem}\label{thm:T-spline-complexes}
Under the assumptions of Proposition~\ref{prop:characterization}, the
following two-dimensional complexes
\begin{equation}\label{eq:Tsplines-complex}
\begin{CD}
\mathbb{R} @>>>\hat Y^{0}_h  @>\hat \,{\bf grad}\,>>  \hat Y^{1}_h @>\hat {\,{\rm rot}\,}>>
\hat Y^{2}_h @> >> 0,
\end{CD}
\end{equation}
\begin{equation}\label{eq:Tsplines-complex-*}
\begin{CD}
\mathbb{R} @>>>\hat Y^{0}_h  @>\hat \,{\bf rot}\,>>  \hat Y^{1*}_h @>\hat {\,{\rm div}\,}>>
\hat Y^{2}_h @> >> 0,
\end{CD}
\end{equation}
where  $\hat {\,{\rm rot}\,} {\bf u} = (\partial_1 u_2 - \partial_2 u_1)$ is the scalar rotor and $\,{\bf rot}\, u =  (\partial_2 u,  - \partial_1 u)^T$ is the vector rotor, are well defined and exact.
\end{theorem}
\begin{proof}
  In the proof we only consider \eqref{eq:Tsplines-complex}, since
  \eqref{eq:Tsplines-complex-*} is equivalent. The well posedness of
  the complex follows from
  \begin{equation}
    \label{eq:proof-step-1}
    \hat \,{\bf grad}\, : \hat Y^{0}_h \rightarrow \hat Y^{1}_h \qquad
    \text{and} \qquad   \hat \,{\rm rot}\, : \hat Y^{1}_h \rightarrow \hat
    Y^{2}_h,
  \end{equation}
which, in turn, easily follows from the definitions
\eqref{eq:T-splines-0-forms-on-parametric}--\eqref{eq:T-splines-2-forms-on-parametric}
and from  Proposition~\ref{prop:characterization}.

Exactness of \eqref{eq:Tsplines-complex} means
\begin{equation}
  \label{eq:exactenss-1}
\mathbb{R} =  \mathrm{ker}( \hat \,{\bf grad}\,),
\end{equation}
\begin{equation}
  \label{eq:exactenss-2}
  \mathrm{im}(\hat \,{\bf grad}\,) =  \mathrm{ker}( \hat \,{\rm rot}\,),
\end{equation}
\begin{equation}
  \label{eq:exactenss-3}
  \mathrm{im}(\hat \,{\rm rot}\,)=  Y^{2}_h.
\end{equation}
The first part, i.e., \eqref{eq:exactenss-1}, is obvious. Moreover
\eqref{eq:exactenss-2} is also simple: indeed if $\hat {\bf u} \in \hat
Y^{1}_h$ has null $\hat \,{\rm rot}\, $, then $ \hat {\bf u} = \hat \,{\bf grad}\, \hat
\phi$, where, e.g.,
\begin{equation}
  \label{eq:phi-h-poincare-T-splines}
  \hat \phi(\zeta_1,\zeta_2) = \int _0^{\zeta_1}  \hat
  u_1(\eta,0)\, d \eta + \int _0^{\zeta_2}  \hat
  u_2(\zeta_1,\eta)\, d \eta.
\end{equation}

Since, $ \hat {\bf u} = \hat \,{\bf grad}\, \hat
\phi$, then  $\hat \phi$  has to be element by element (of $\mathcal{M}_{ext}$)
a $p$-degree tensor-product polynomial. Then, $\hat \phi$ inherits  the interelement regularity from $\hat {\bf u}$ and has the one of functions in  $ \hat Y^{0}_h$.  Then, by Proposition~\ref{prop:characterization}\color{black},
 $\hat \phi \in \hat Y^{0}_h$.  The last
point, \eqref{eq:exactenss-3},  follows from the dimension formula
\begin{equation}
  \label{eq:proof-step-2}
  \# (\hat Y^{0}_h ) +   \# (\hat Y^{2}_h ) =   \# (\hat Y^{1}_h ) + 1.
\end{equation}
Indeed, using \eqref{eq:exactenss-2}, \eqref{eq:exactenss-1}, and \eqref{eq:proof-step-2},
\begin{displaymath}
  \begin{aligned}
        \# ( \mathrm{im}(\hat \,{\rm rot}\,)) &=   \# (\hat Y^{1}_h )  - \# (
    \mathrm{ker}( \hat \,{\rm rot}\,))  \\
&=   \# (\hat Y^{1}_h )  - \# ( \mathrm{im}(\hat \,{\bf grad}\,) )) \\
&=   \# (\hat Y^{1}_h )  - \# ( \hat Y^{0}_h) +\# ( \mathrm{ker}(\hat \,{\bf grad}\,) ))\\
&=   \# (\hat Y^{1}_h )  - \# ( \hat Y^{0}_h) +1\\
& = \# (\hat Y^{2}_h ) .
  \end{aligned}
\end{displaymath}
In order to prove \eqref{eq:proof-step-2}, we recall the  Euler's formula for the T-mesh $\mathcal{M}$
\begin{equation}\label{eq:euler}
F_0 + V_0 = E_0 + 1,
\end{equation}
where $F_0$ is the number of faces, $E_0$ the number of edges and $V_0$ the number of vertices of $\mathcal{M}$, including knot repetitions, zero length edges and empty elements. The proof is different for odd and even $p$.

\begin{trivlist}
\item[Case 1)] Let $p$ be odd. We can separate the edges into
  horizontal and vertical ones, and with self-explaining notation we
  have $E_0 = E_0^H + E_0^V$. Similarly, the vertices can be divided
  into horizontal T-junctions, vertical T-junctions and all the other
  vertices (including those on the boundary), in the form $V_0 = V_0^H
  + V_0^V + V_0^+$. For odd $p$ the meshes ${\mathcal M}^{1}_1$ and ${\mathcal M}^{1}_2$ are
  constructed by adding the first-bay face-extension of horizontal and
  vertical T-junctions, respectively. Thus, using the assumption that
  T-junction extensions do not intersect,  the number of horizontal edges in ${\mathcal M}^{1}_1$ is $E_1^H = E_0^H + V_0^H$, and the number of vertical edges in ${\mathcal M}^{1}_2$ is $E_1^V = E_0^V + V_0^V$. Similarly, the mesh ${\mathcal M}^{2}$ is contructed by adding all the first-bay face-extensions, and the number of faces in ${\mathcal M}^{2}$ is equal to $F_2 = F_0 + V_0^H + V_0^V$.

Since $p$ is odd, and from the positions of the anchors in every mesh (see Figure~\ref{fig:Tmesh-p3}), the dimensions of the spaces are
\begin{equation*}
\# (\hat Y^{0}_h) = V_0, \quad \# (\hat Y^{1}_h) = E_1^H + E_1^V = E_0 + V_0^H + V_0^V, \quad \#(\hat Y^{2}_h) = F_2 = F_0 + V_0^H + V_0^V,
\end{equation*}
and using \eqref{eq:euler} the proof is finished.

\item[Case 2)] Let $p$ be even. We denote by $V_0^B$ and $E_0^B$ the number of boundary vertices and boundary edges in $\mathcal{M}$, and we note that $V_0^B = E_0^B$. As before, we distinguish between horizontal and vertical edges, $E_0 = E_0^V + E_0^H$, and also for the boundary edges $E_0^B = E_0^{B,V} + E_0^{B,H}$. For even $p$ the mesh ${\mathcal M}^{1}_1$ (resp. ${\mathcal M}^{1}_2$) is constructed by removing the first and last columns (resp. rows) of elements from $\mathcal{M}$. Hence, the number of vertical edges in ${\mathcal M}^{1}_1$ is $E_1^V = E_0^V - E_0^{B,V}$, and the number of horizontal edges in ${\mathcal M}^{1}_2$ is $E_1^H = E_0^H - E_0^{B,H}$. Similarly, the mesh ${\mathcal M}^{2}$ is constructed by removing the first and last rows and columns of elements from $\mathcal{M}$, thus the number of vertices in ${\mathcal M}^{2}$ is $V_2 = V_0 - V_0^B$.

From the position of the anchors for even $p$ (see Figure~\ref{fig:Tmesh-p2}), the dimensions of the spaces are
\begin{equation*}
\# (\hat Y^{0}_h) = F_0, \quad \# (\hat Y^{1}_h) = E_1^V + E_1^H = E_0 - E_0^B, \quad \#(\hat Y^{2}_h) = V_2 = V_0 - V_0^B.
\end{equation*}
Using \eqref{eq:euler} and that $V_0^B = E_0^B$ the proof is finished. \qed
\end{trivlist}

\end{proof}

\subsection{Three-dimensional De Rham complex based on  T-splines and B-splines} \label{sect:3D-T-splines}

We construct a three-dimensional complex on the parametric domain  by
tensor product of the two-dimensional T-spline complexes
\eqref{eq:Tsplines-complex}--\eqref{eq:Tsplines-complex-*} and the
 one-dimensional complex \eqref{eq:1D-diagram}. Then we define the spaces on the
 parametric domain $\hat \Omega = (0,1)^3$:
\begin{equation}
    \label{eq:Tsplines-3D-discrete-forms-on-omegahat}
  \begin{aligned}
\hat{X}^0_h &:= \hat Y^{0}_h\otimes S_p(\Xi), \\
\hat{X}^1_h &:= [\hat Y^{1}_h\otimes S_p(\Xi)] \times [ \hat Y^{0}_h\otimes S_{p-1}(\Xi')] ,  \\
\hat{X}^2_h &:=  [\hat Y^{1*}_h\otimes S_{p-1}(\Xi')] \times [ \hat Y^{2}_h\otimes S_p(\Xi)], \\
\hat{X}^3_h &:= \hat Y^{2}_h\otimes  S_{p-1}(\Xi') ;
\end{aligned}
\end{equation}
 which form a  complex of the kind \eqref{eq:diag1} (or \eqref{eq:diag2}
 if we also impose homogeneous Dirichlet boundary conditions).

Assume now that the geometry map ${\bf F}$ is
 tensor-product single-patch spline or NURBS, and fulfills Assumption~\ref{ass:F}, now with $\hat{X}^0_h$ defined as in \eqref{eq:Tsplines-3D-discrete-forms-on-omegahat}. Therefore, the push-forwards
 \eqref{eq:diag1}--\eqref{eq:diag2} give the correct \color{black}  complex
 $(X^0_h  $, \ldots,  $X^3_h  )$ \color{black} on $\Omega$:
 this procedure is completely analogous to what we have already
 described in Section~\ref{sec:compatible -splines-phys-domain} and
 is not detailed here.

It is not a difficulty to consider,  more generally, a multi-patch, or a  T-spline geometry
mapping. This is not detailed here, for the sake of brevity, but the
first case will be addressed in the numerical tests of the next section.

\subsection{Concluding remarks on the T-spline complex}

As it appears from our presentation, the understanding of the T-spline complex is much less sound than the one of the spline complex, even in two space dimensions. Moreover some of the properties we have studied for splines do not hold in general for T-splines. For example,
\begin{itemize}
\item the matrices corresponding to the operators are no more the incidence matrices of the mesh $\mathcal{M}$ and a similar fact is true for standard finite elements with hanging nodes, i.e., the T-spline  complex with $p=1$;
\item the definition of control mesh and control fields is not trivial especially when $p$ is even and the analogue of Section \ref{sec:dofs} is not available for T-splines. This deserves further studies.
\end{itemize}

\section{Numerical results}
\label{sec:num}

In this section we present numerical tests showing the behavior of isogeometric methods for electromagnetic problems.
Since numerical tests for B-splines have already been presented in other works, see e.g., \cite{VB10, Buffa_Sangalli_Vazquez, BRSV11}, we will concentrate here on examples involving also T-splines. All our numerical tests have been performed with the Matlab library GeoPDEs \cite{dFRV11}. It should be said though that GeoPDEs does not have full T-splines capability, and in particular does not provide any T-splines adaptivity in the sense of \cite{Scott_Li_Sederberg_Hughes}.

The mappings we use in this section always verify the Assumption in Section~\ref{sect:3D-T-splines}, and can be either single-patch or
multi-patch; the meshes we describe are the ones corresponding to the space $X^0_h$. The meshes for the other spaces are
constructed following the procedure detailed in
Section~\ref{sect:T-compatible-parametric}. In the figures of this
section, repeated lines of the mesh are represented with thicker
lines, independently of the number of repetitions. In all cases,
internal mesh lines have multiplicity one.

 \subsection{Maxwell eigenproblem in the square domain} \label{sect:square}

As a first test we solve the two-dimensional eigenvalue problem: \textit{Find $({\bf u},\omega) \in \,{\bf H}_0({\rm rot};\Omega)\, \times \mathbb{R}$ such that}
\begin{equation} \label{eq:eigenproblem_2D}
\int_\Omega \,{\rm rot}\, {\bf u} \, \,{\rm rot}\, \overline {\bf v} = \omega^2 \int_\Omega {\bf u} \cdot \overline {\bf v} \quad \forall {\bf v} \in \,{\bf H}_0({\rm rot};\Omega)\,,
\end{equation}
in the square domain $\Omega = (0,\pi)^2$, for which the exact
eigenvalues are $ \omega^2 = m^2 + n^2$, with $m,n = 0, 1, \ldots$. The aim of this test is to show that the discretization of the problem with T-splines does not present spurious modes.

The coarsest mesh consists of 8 square (non-empty)
elements in the left half, and 4 rectangular (non-empty) elements in
the right half, thus creating several T-junctions on the vertical line
$\zeta_1 = 0.5$. Finer meshes are created by dividing each element into 4 (see Figure~\ref{fig:mesh_square}).

In Table~\ref{tab:eigv} we present the first non-null eigenvalues for degree~3 and for the sequence of meshes explained above. The results show that there are no spurious eigenvalues, and that a good convergence rate is obtained. In Figure~\ref{fig:square_eigv} we display the first non-null eigenvalues computed with discretizations of degree~4 and~5 in a mesh formed by 768~non-empty elements, and its comparison with the exact eigenvalues. Again, it is seen that the discrete eigenvalues are computed with the right multiplicity.

\begin{figure}[h!]
\centering
\includegraphics[width=.4\textwidth]{T-mesh_square_border_coarse.pdf}
\includegraphics[width=.4\textwidth]{T-mesh_square_border.pdf}
\caption{Coarsest mesh for the square, and mesh after one refinement step.} \label{fig:mesh_square}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=.48\textwidth]{eigv_deg45b.pdf}
\caption{First non-null eigenvalues computed in the square for degrees 4 and 5. } \label{fig:square_eigv}
\end{figure}

\begin{table}[!htbp]
\centering
{\footnotesize
\begin{tabular}{|c|c||c|c|c|c|c|}
\hline
Mode & Exact & \multicolumn{5}{|c|}{Computed} \\
\hline
(1,0) & 1.00000  &1.00001 &  1.00000 &  1.00000  & 1.00000  & 1.00000\\
(0,1) & 1.00000  &1.00005 &  1.00000 &  1.00000  & 1.00000  & 1.00000\\
(1,1) & 2.00000  &2.00016 &  2.00000 &  2.00000  & 2.00000  & 2.00000\\
(2,0) & 4.00000  &4.00396 &  4.00004 &  4.00000  & 4.00000  & 4.00000\\
(0,2) & 4.00000  &4.03882 &  4.00134 &  4.00002  & 4.00000  & 4.00000\\
(2,1) & 5.00000  &5.00395 &  5.00003 &  5.00000  & 5.00000  & 5.00000\\
(1,2) & 5.00000  &5.10164 &  5.00208 &  5.00002  & 5.00000  & 5.00000\\
(2,2) & 8.00000  &8.05454 &  7.99989 &  8.00001  & 8.00000  & 8.00000\\
(3,0) & 9.00000  &9.06255 &  9.00135 &  9.00001  & 9.00000  & 9.00000\\
(0,3) & 9.00000  &9.12399 &  9.02102 &  9.00057  & 9.00001  & 9.00000\\
(3,1) & 10.0000  &10.0614 &  10.0014 &  10.0000  & 10.0000  & 10.0000\\
(1,3) & 10.0000  &10.2361 &  10.0324 &  10.0007  & 10.0000  & 10.0000\\
(3,2) & 13.0000  &12.8159 &  13.0028 &  13.0000  & 13.0000  & 13.0000\\
(2,3) & 13.0000  &13.2002 &  13.0091 &  13.0004  & 13.0000  & 13.0000\\
(4,0) & 16.0000  &17.9413 &  16.0181 &  16.0002  & 16.0000  & 16.0000\\
(0,4) & 16.0000  &19.8934 &  16.2962 &  16.0076  & 16.0001  & 16.0000\\
(4,1) & 17.0000  &19.9586 &  17.0181 &  17.0002  & 17.0000  & 17.0000\\
(1,4) & 17.0000  &20.8937 &  18.0245 &  17.0092  & 17.0001  & 17.0000\\
(3,3) & 18.0000  &21.4707 &  18.7373 &  18.0008  & 18.0000  & 18.0000\\
(4,2) & 20.0000  &24.0689 &  20.0191 &  20.0002  & 20.0000  & 20.0000\\
(2,4) & 20.0000  &26.1844 &  21.6138 &  20.0056  & 20.0001  & 20.0000\\
\hline \multicolumn{2}{|c||}{d.o.f.}  & 74 & 184 & 548 & 1852 & 6764 \\
\hline \multicolumn{2}{|c||}{number of zeros} & 21 & 65 & 225 & 833 & 3201 \\
\hline
\end{tabular}
}
\caption{First non-null eigenvalues computed in the square for $p=3$.}
\label{tab:eigv}
\end{table}

\begin{remark}
We have also solved the previous test by the mixed formulations in
\cite{Boffi07}, that make use of the full two-dimensional De Rham
complex \eqref{eq:Tsplines-complex}. The computed non-null
eigenvalues are the same as for the plain formulation
\eqref{eq:eigenproblem_2D}, while the zero eigenvalues are filtered   with the mixed formulation\color{black}. These results, that we do not present here for the sake of brevity, confirm that the construction of the De Rham complex with T-splines is correct.
\end{remark}

\subsection{Maxwell eigenproblem in the thick L-shaped domain} \label{sect:thickL}
As a second test case, we solve the three-dimensional eigenvalue problem: \textit{Find $({\bf u},\omega) \in \,{\bf H}_0(\mathbf{curl};\Omega)\, \times \mathbb{R}$ such that}
\begin{equation}
  \label{eq:eigenvalues-thickL-3D}
  \int_\Omega \,{\bf curl}\, {\bf u} \cdot\,{\bf curl}\, \overline {\bf v} = \omega^2 \int_\Omega {\bf u} \cdot \overline {\bf v} \quad \forall {\bf v} \in \,{\bf H}_0(\mathbf{curl};\Omega)\,,
\end{equation}
in the thick L-shaped domain $\Omega = \Sigma \times (0,1)$, where $\Sigma = (-1,1)^2 \setminus [-1,0]^2$. From \cite{CD2000}, it is known that the reentrant edge introduces a singularity in the first eigenfunction, which only belongs to the space $H^{2/3-\varepsilon}(\Omega)$ for any $\varepsilon > 0$.

It is well known that in order to recover the optimal convergence
rate we need to suitably refine the mesh toward the reentrant edge,
see e.g., \cite{MR1860445} or \cite{BuCoDa05}. Anisotropic
elements need to be used in this case (see \cite{Beirao_Cho_Sangalli} for some theoretical
background on the topic). We propose here a
dyadic refinement based on T-splines.

For the geometry representation, the thick L-shaped domain is parametrized as
the union of three cubic patches.  Following Section  \ref{sec:multi-patch-complex-continuity}, scalar fields in $X^0_h$ are only continuous at the interfaces between patches, and the fields in $X^1_h$, which are used in the
discretization of \eqref{eq:eigenvalues-thickL-3D}, are  only
tangentially continuous at these interfaces (like for standard
edge finite elements), but at least $C^{p-2}$ within patches.

The refinement is obtained via T-splines by dyadic partitioning of elements which are close to the reentrant edges \cite[Ch.~4]{Apel}.   We perform the refinement first in an L-shaped two-dimensional section, and then propagate to the three-dimensional domain with a uniform mesh in the $z$-direction, as already explained in Section~\ref{sect:3D-T-splines}. The refinement is performed identically for every patch, in such a way that conformity can be ensured at the patch interfaces.

To construct the two-dimensional mesh, at each refinement step, and for every patch, we refine a small square region near the reentrant edge, subdividing each element into 4. Then some T-junction extensions are added, depending on the degree, to make the mesh analysis suitable, as defined in Section~\ref{sect:AS}. For instance, in the example of Figure~\ref{fig:mesh_thickL} we start with an $8\times8$ mesh for each patch, which is drawn in black. At the first step we refine a region of $3 \times 3$ elements on each patch. Since the degree is $p=4$, two-bay extensions must be added to make the mesh analysis suitable. The refined elements at this step are given by the blue lines. At the second refinement step, which is marked in red, we first refine a square region of $2 \times 2$ elements, and again we add the two-bay extensions to make the mesh analysis suitable. Finally we remark that, since the dyadic partition and the definition of analysis suitable T-mesh depend on the degree, different meshes are used for different degrees. \color{black}

\begin{figure}[h!]
\centering
\includegraphics[width=.4\textwidth]{mesh_L_shaped.pdf}
\includegraphics[width=.45\textwidth]{mesh_thickL.pdf}
\caption{Example of a two-dimensional mesh in the L-shaped domain $\Sigma$, and its extension to the three-dimensional domain $\Omega = \Sigma\times (0,1)$ , for $p=4$.} \label{fig:mesh_thickL}
\end{figure}

The problem has been solved for degrees~$4$ and~$5$. In Tables~\ref{tab:thickL_deg4} and~\ref{tab:thickL_deg5} we present the first non-null computed eigenvalues in the three cases, and its comparison with the exact solution. In Figure~\ref{fig:thickL_eigv} we display the convergence rate for the first eigenvalue, comparing the results obtained with T-splines and with a B-spline discretization of the same degree in the corresponding refined tensor-product dyadic mesh. As can be seen, with T-splines we obtain the same error with an important reduction in the number of the degrees of freedom.

\begin{figure}[h!]
\centering
\begin{subfigure}[Degree 4]{\includegraphics[width=.48\textwidth]{thickL_eigv1_deg4.pdf}}
\end{subfigure}
\begin{subfigure}[Degree 5]{\includegraphics[width=.48\textwidth]{thickL_eigv1_deg5.pdf}}
\end{subfigure}
\caption{Convergence of the first eigenvalue in the thick L-shaped domain.} \label{fig:thickL_eigv}
\end{figure}

\begin{table}[!htbp]
\centering
{\footnotesize
\begin{tabular}{|c||c|c|c|c|c|c|}
\hline
Exact & \multicolumn{6}{|c|}{Computed} \\
\hline
9.63972384472	&9.64482260735	&9.64055367165	&9.63986647533	&9.63977706731	&9.63974511214	&9.63972731966 \\
\hline
11.3452262252	&11.3444193267	&11.3450973393	&11.3452056015	&11.3452178503	&11.3452228875	&11.3452256921 \\
\hline
13.4036357679	&13.4036330719	&13.4036359208	&13.4036359870	&13.4036357431	&13.4036357654	&13.4036357699 \\
\hline
15.1972519265	&15.1973643408	&15.1973310163	&15.1973301300	&15.1972556440	&15.1972524223	&15.1972523662 \\
\hline
19.5093282458	&19.5144198480	&19.5101576732	&19.5094708082	&19.5093814308	&19.5093494993	&19.5093317180 \\
\hline
19.7392088022	&19.7392474090	&19.7392464705	&19.7392464473	&19.7392098765	&19.7392090606	&19.7392090522 \\
\hline
19.7392088022	&19.7392474115	&19.7392464714	&19.7392464480	&19.7392098765	&19.7392090617	&19.7392090536 \\
\hline
19.7392088022	&19.7392854949	&19.7392835833	&19.7392835402	&19.7392109156	&19.7392092780	&19.7392092574 \\
\hline
21.2590837990	&21.2591164396	&21.2591199815	&21.2591200740	&21.2590848357	&21.2590840605	&21.2590840611 \\
\hline
d.o.f. &        4042    &    7126   &    11018   &    16162   &    22630    &   34894\\
\hline
\end{tabular}
}
\caption{First non-null eigenvalues computed in the thick L-shaped domain for $p=4$.}
\label{tab:thickL_deg4}
\end{table}

\begin{table}[!htbp]
\centering
{\footnotesize
\begin{tabular}{|c||c|c|c|c|c|}
\hline
Exact & \multicolumn{5}{|c|}{Computed} \\
\hline
9.63972384472	&9.64328299807	&9.64030443177	&9.63981618307	&9.63973902051	&9.63973012738 \\
\hline
11.3452262252	&11.3446611860	&11.3451346157	&11.3452117088	&11.3452238284	&11.3452252153 \\
\hline
13.4036357679	&13.4036342774	&13.4036357233	&13.4036357613	&13.4036357622	&13.4036357630 \\
\hline
15.1972519265	&15.1972704673	&15.1972555872	&15.1972551905	&15.1972551805	&15.19725206245 \\
\hline
19.5093282458	&19.5128817631	&19.5099083654	&19.5094205444	&19.5093433951	&19.5093345086 \\
\hline
19.7392088022	&19.7392095886	&19.7392095763	&19.7392095758	&19.7392095758	&19.7392088273 \\
\hline
19.7392088022	&19.7392095886	&19.7392095763	&19.7392095758	&19.7392095758	&19.7392088273 \\
\hline
19.7392088022	&19.7392103678	&19.7392103456	&19.7392103444	&19.7392103444	&19.7392088514 \\
\hline
21.2590837990	&21.2590824582	&21.2590845213	&21.2590845754	&21.2590845767	&21.2590838179 \\
\hline
d.o.f. &              5891     &   9883   &    14827   &    20723   &    28105\\
\hline
\end{tabular}
}
\caption{First non-null eigenvalues computed in the thick L-shaped domain for $p=5$.}
\label{tab:thickL_deg5}
\end{table}

\subsection{Maxwell source problem in three quarters of the cylinder}
For the third test case we consider the model source problem: \textit{Find ${\bf u} \in \,{\bf H}_{0,\Gamma_D}({\mathbf{curl}};\Omega)\,$ such that}
\begin{equation}
\int_\Omega \,{\bf curl}\, {\bf u} \cdot \,{\bf curl}\, \overline {\bf v} + \int_\Omega {\bf u} \cdot \overline {\bf v} = \int_\Omega {\bf f} \cdot \overline {\bf v} \quad \forall {\bf v} \in \,{\bf H}_{0,\Gamma_D}({\mathbf{curl}};\Omega)\,,
\end{equation}
where $\,{\bf H}_{0,\Gamma_D}({\mathbf{curl}};\Omega)\,$ is the set of functions with null tangential trace on
$\Gamma_D \subset \partial \Omega$, i.e., \[\,{\bf H}_{0,\Gamma_D}({\mathbf{curl}};\Omega)\, := \{{\bf v} \in \,{\bf H}({\mathbf{curl}};\Omega)\, : {\bf v} \times {\bf n} =
\,{\bf 0} \text{ on } \Gamma_D\}. \]

The geometry $\Omega$ is three quarters of a cylinder of radius and length equal to one (see Figure~\ref{fig:mesh_camembert}), that in cylindrical coordinates is given by $\Omega = \{(r, \theta, z): 0<r<1, \, 0 <\theta<\frac{3}{2}\pi, \, 0 <z<1\}$. We impose the null tangential component on $\Gamma_D = \{(r,\theta,z) : \theta \in \{0,\frac{3}{2}\pi\}\}$, and the source term ${\bf f}$ is taken such that the exact solution is ${\bf u} = \,{\bf grad}\, (r^{2/3} \sin (2\theta/3) \sin(\pi z))$, i.e., it is singular in the first two directions, but it is regular in the $z$ direction, for which the local refinement of the previous test is well suited for this case.

As in the previous example, the domain is defined with three patches,
and the discrete fields $X^1_h \subset \,{\bf H}({\mathbf{curl}};\Omega)\, $ are only tangentially continuous between them. The construction of the mesh in the parametric domain is identical to the one in the previous example: for each patch we first create a two-dimensional mesh locally refined towards the corner, and extend it to the three dimensional domain by tensor product. The mesh is then mapped to the physical domain, as can be seen in Figure~\ref{fig:mesh_camembert}.

\begin{figure}[h!]
\centering
\includegraphics[width=.65\textwidth]{mesh_camembert2.pdf}
\vspace*{-2cm}
\caption{Example of a mesh for three quarters of the cylinder.} \label{fig:mesh_camembert}
\end{figure}

The problem is solved with T-splines of degree~3, and also with B-splines of the same degree in the corresponding tensor product mesh. The errors in ${\bf H}(\bf curl)$-norm for the two methods are compared in Figure~\ref{fig:error_camembert}. As in the previous example, with T-splines we are able to obtain results similar to those of B-splines with a reduction of the number of degrees of freedom.

\begin{figure}[h!]
\centering
\includegraphics[width=.4\textwidth]{camembert_error.pdf}
\caption{Comparison of the error for T-splines and B-splines.} \label{fig:error_camembert}
\end{figure}

\subsection{Numerical simulation of a twisted waveguide}

As the last numerical test we use T-splines to simulate the
propagation of a singular mode in a waveguide with a twist. The
configuration, which is presented in Figure~\ref{fig:waveguide}, is
the same given in \cite[Ch.~8]{Jin}, changing the material
discontinuity by a geometric inhomogeneity (the twist). The problem is solved in a waveguide with a twist of 90~degrees, with a section of three quarters of the circle of 2~cm radius, and for which we assume that the walls are perfect electrical conductors. We also assume that the waveguide extends to infinity without other inhomogeneities, and it is truncated by the planes $\Gamma_1$ and $\Gamma_2$ to obtain the computational domain, which consists of three different regions: a middle region where the waveguide is twisted (see Figure~\ref{fig:waveguide}), and two straight regions near the ports, to keep the inhomogeneity far enough from them, in such a way that only the dominant mode TE$_{10}$ can propagate without attenuation. The total length of the computational domain is 24~cm: 4~cm for each straight region, and 16~cm for the region with the twist. The frequency $\omega$ is taken equal to 32~GHz, and it is between the cutoff frequencies for the first mode (21~GHz) and the second mode (33.84~GHz).
\color{black}

Following \cite{Jin}, and working in the time harmonic regime at a given frequency $\omega$, the (complex-valued) electric field ${\bf E} \in \,{\bf H}_{0,\Gamma_D}({\mathbf{curl}};\Omega)\,$ is solution of the problem
\begin{equation}\label{eq:vf_waveguide}
\begin{array}{rcl}
\displaystyle \int_\Omega ( \,{\bf curl}\, {\bf E} \cdot \,{\bf curl}\, \overline {\bf G} - k^2 {\bf E} \cdot \overline {\bf G}) & + & \displaystyle \int_{\Gamma_1 \cup \Gamma_2} {\rm i} \beta_{10} ({\bf n} \times {\bf E}) \cdot ({\bf n} \times \overline {\bf G}) = \\
&& \displaystyle 2{\rm i} \beta_{10} \, \int_{\Gamma_1} {\bf E}^{\rm inc} \cdot \overline {\bf G}, \quad \forall {\bf G} \in \,{\bf H}_{0,\Gamma_D}({\mathbf{curl}};\Omega)\,,
\end{array}
\end{equation}
where
$k = \sqrt{\omega^2 \mu_0 \varepsilon_0}$ with $\mu_0$ and $\varepsilon_0$ the magnetic permeability and electric permittivity of free space. The incident electric field ${\bf E}^{\rm inc}$ at the port $\Gamma_1$, and the wavenumber of the first mode $\beta_{10}$ are defined as
\begin{align*}
{\bf E}^{\rm inc}(x,y,z) = {\bf e}_{10}(x,y) e^{-{\rm i} \beta_{10} z}, &&
\beta_{10} = \sqrt{k^2 - k_{10}^2}.
\end{align*}
In the case of waveguides of rectangular or circular section, the
value of the constant $k_{10}$ and the mode ${\bf e}_{10}$ are known. In
the general case, they can be obtained by solving a 2D eigenvalue
problem on the port $\Gamma_1$, which consists on finding the minimum eigenvalue $k_{10} \in \mathbb{R}$, and its associated eigenvector ${\bf e}_{10} \in {\bf H}_0(\,{\rm rot}\,;\Gamma_1)$, such that
\begin{equation}\label{eq:2d_waveguide}
\int_{\Gamma_1} \,{\rm rot}\, {\bf e}_{10} \, \,{\rm rot}\, \overline {\bf v} = k_{10}^2 \int_{\Gamma_1} {\bf e}_{10} \cdot {\bf v} \quad \forall {\bf v} \in {\bf H}_0(\,{\rm rot}\,;\Gamma_1).
\end{equation}

The electric field ${\bf E}$ in equation \eqref{eq:vf_waveguide} is discretized with T-splines of degree~3, using the approach already explained in Section~\ref{sect:3D-T-splines}. The two-dimensional T-mesh for the section is built as in the previous examples, and in the $z$ direction we use one element for each straight region near the ports, and 4 elements along the twist, for a total of 7936 degrees of freedom. For the solution of the 2D problem \eqref{eq:2d_waveguide}, it is enough to restrict a field in $\,{\bf H}({\mathbf{curl}};\Omega)\,$ to its tangential components on the port $\Gamma_1$, which in practice is equivalent to solve with two-dimensional T-splines.

The magnitude of the real part of the computed solution ${\bf E}$ is shown in Figure~\ref{fig:waveguide_mode}, which shows that the mode is correctly propagated. Finally, we also compute the reflection and transmission coefficients, given by the equations
\begin{equation*}
R = \frac{e^{-{\rm i} \beta z_1} \int_{\Gamma_1} {\bf E} \cdot {\bf e}_{10}}{\int_{\Gamma_1} {\bf e}_{10}\cdot {\bf e}_{10}} - e^{-2 {\rm i} \beta z_1}, \qquad T = \frac{e^{{\rm i} \beta z_2} \int_{\Gamma_1} {\bf E} \cdot {\bf e}_{10}}{\int_{\Gamma_1} {\bf e}_{10}\cdot {\bf e}_{10}} - e^{-2 {\rm i} \beta z_2},
\end{equation*}
and we obtain the values $|R| = 0.0025$ and $|T| = 0.9998$, which confirms that the twist does not affect the propagation of the mode, as expected.

\begin{figure}[h!]
\begin{subfigure}[Geometry of the waveguide]
{\includegraphics[width=.48\textwidth]{waveguide.png} \label{fig:waveguide}}
\end{subfigure}
\begin{subfigure}[Real part of the computed electric field]
{\includegraphics[width=.45\textwidth]{realE.png} \label{fig:waveguide_mode}}
\end{subfigure}
\caption{Geometry of the waveguide and real part of the computed electric field.}
\end{figure}

\bibliographystyle{elsarticle-num}
\def$'$} \def\cprime{$'${$'$} \def$'$} \def\cprime{$'${$'$}