O2 OΞ^Lπ_#1,#2 O2 OΞ^Lπ_#1,#2^per O2 OΞ^Lπ_#1,#2^c,per Oi ^LΠ^#1 Oi LΠ^#1 Oi
^LΠ^#1 Oi LΠ^#1 Oi ^LΠ^#1 hypothesisHypothesis claimClaim Energy conservative
isogeometric methods for the wave equationA.Bressan, A. Buffa, A. Kushova, and
R. Vázquez ex_supplement

# Energy conservative isogeometric methods for the wave equation††thanks:
Submitted to the editors DATE.

Andrea Bressan IMATI-CNR Pavia<EMAIL_ADDRESS>Annalisa Buffa
École Polytechnique Fédérale de Lausanne<EMAIL_ADDRESS><EMAIL_ADDRESS>Alen Kushova Università degli studi di Pavia ().
<EMAIL_ADDRESS>Rafael Vázquez 33footnotemark: 3

###### Abstract

We analyze the wave equation in mixed form, with periodic and/or Dirichlet
homogeneous boundary conditions, and nonconstant coefficients that depend on
the spatial variable. For the discretization, the weak form of the second
equation is replaced by a strong form, written in terms of a projection
operator. The system of equations is discretized with B-splines forming a De
Rham complex along with suitable commutative projectors for the approximation
of the second equation. The discrete scheme is energy conservative when
discretized in time with a conservative method such as Crank-Nicolson. We
propose a convergence analysis of the method to study the dependence with
respect to the mesh size $h$, with focus on the consistency error. Numerical
results show optimal convergence of the error in energy norm, and a relative
error in energy conservation for long-time simulations of the order of machine
precision.

###### keywords:

Isogeometric Analysis, wave equation, splines, quasi-interpolants, energy
conservation.

##  1 Introduction

The purpose of this paper is to provide a numerical analysis of discretization
schemes applied to the anisotropic wave equation in mixed form, in which one
of the variational equations is replaced by a suitably projected equation
within the discrete space, introducing a consistency error. The motivation for
our work stems from analogous discretization techniques recently proposed for
Maxwell’s equations, specifically in the long-term study and simulation of
particle dynamics in plasma, as well as for the electrodynamics of phenomena
resulting from particle motion [23], [21]. The simplified model problem of the
wave equation allows us to analyze the proposed method, paying particular
attention to the introduced consistency error.

The wave equation has been studied extensively in theory and numerical
approximations. Existence and uniqueness results are well known in literature
[11],[18]. The present work considers the wave equation as a first order
hyperbolic system, introducing velocity and pressure fields following
approaches found in prior works [3], [19], [8]. For the resulting variational
formulations in mixed form it is acknowledged that the successful
approximation of the problem requires suitable compatibility conditions within
the involved discrete spaces [7]. In general, mixed methods for the wave
equation consider the discretization of vector fields in some
$\mathbf{H}(\textup{div})$-conforming spaces while scalar fields in some
$L^{2}$ spaces. To meet these conditions, we construct discrete spaces that
adhere to the De Rham complex of exterior calculus [1], within the framework
of isogeometric analysis [22].

Briefly, Isogeometric Analysis (IgA) utilizes spline functions, or their
generalizations, for both representing the computational domain and
approximating solutions to the partial differential equation that models the
relevant problem. This approach aims to streamline the interoperability
between computer-aided design and numerical simulations. Moreover, IgA derives
advantages from the approximation properties of splines, where their high
continuity contributes to enhanced accuracy compared to $C^{0}$ piecewise
polynomials. This characteristic is well-documented in literature, [14, 9,
32]. As regards the De Rham complex for tensor-product B-splines, initially
introduced and analyzed by [10], it has found wide applications in Galerkin
approximation of Maxwell’s equations [31], and divergence-free methods for
incompressible fluid flow [15], [16],[17].

In this paper, together with the discrete B-spline spaces of the De Rham
complex, we build quasi-interpolant projections that commute with the
divergence operator. Various families of quasi-interpolant operators have been
defined and presented in the context of spline approximation, [13],[27]. Here,
we implement a local quasi-interpolant operator, as presented in [25], wherein
explicit formulas are provided for computing the coefficients of the
projection. These explicit formulas entail pointwise evaluation of the
function to be projected. Employing this operator, we project the second
equation of the variational formulation, reducing the problem to a single
equation for the velocity. Subsequently, the pressure field is determined
after the computation of velocity. The resulting semi-discrete problem
preserves the total energy of the system, as for a standard Galerkin method.

Given our interest in preserving energy for long time simulations, the choice
of an energy-preserving method in time is mandatory. We adopt Crank-Nicolson
which is a second order energy conservative time discretization. Other time
discretization methods that conserve energy can be applied with appropriate
modifications to the fully discrete system. The numerical tests show a
relative error in energy conservation for long time simulations of the order
of machine precision.

Finally, error estimates in energy norm for Galerkin discretizations of mixed
formulations are well known in literature [7], and in particular for the wave
equation in mixed form [8]. Here, as the main contribution of this work, we
present an error convergence analysis for our method with a generic family of
projections, focusing on the consistency error. By assuming good approximation
properties of the projections, together with its formal adjoint operator, we
proof high-order convergence with respect to the mesh size $h$. Since the
specific quasi-interpolant that we tested numerically does not fullfil the
requires on the adjoint operator, we relax the assumption and proof that the
numerical scheme converges linearly with respect to the mesh size. Numerical
simulations confirm optimal high-order convergence for the implemented quasi-
interpolant, similar to a standard Galerkin scheme, which suggests an
improvement of our theoretical results can be achieved. The global convergence
of the scheme is of second order, as expected from the use of Crank-Nicolson
method.

The paper’s structure is organized as follows. In section 2 we present the
model problem. In section 3 we provide a brief overview of the isogeometric
framework. Section 4 introduces the discretization spaces along with
commutative projections and presents the discrete problem. In section 5 we
present the a priori error estimates in energy norm. Implementation details
are covered in section 6, with numerical results presented in section 7, and
finally, in section 8 we draw our conclusions.

##  2 Wave equation

We start presenting the wave equation and its formulation as a first order
hyperbolic system.

###  2.1 Strong formulation and boundary conditions

Let the domain $\Omega\subset\mathbb{R}^{d}$, with $d=2,3$, and let the time
domain interval $I=[0,T]$ with $T>0$. We consider our model problem, the
second-order scalar wave equation with space dependent coefficients, that
reads like

(1) $\begin{cases}u_{tt}-\text{div}(c^{2}\gradient{u})=0&\text{in
}\Omega\times I,\\\ u(\cdot,0)=u_{0},\quad u_{t}(\cdot,0)=u_{1}&\text{in
}\Omega,\end{cases}$

with the subscript $\cdot_{t}$ indicating the partial derivative with respect
to time, $u_{0}$ and $u_{1}$ are the initial conditions, and the coefficient
$c$ is a positive, uniformly bounded and smooth scalar field in $\Omega$. The
problem has to be completed with boundary conditions, that we will detail
below. Following the idea in [19], by introducing the new variables,
$\mathbf{v}=c\gradient{u}$ and $\phi=u_{t}$, we rewrite (1) as a first order
hyperbolic system, in the form

(2) $\begin{cases}\mathbf{v}_{t}=c\gradient{\phi}&\text{in
}\Omega\times[0,T],\\\ \phi_{t}=\text{div}(c\mathbf{v})&\text{in
}\Omega\times[0,T],\\\
\mathbf{v}(\cdot,0)=c\gradient{u_{0}},\quad\phi(\cdot,0)=u_{1}&\text{in
}\Omega.\end{cases}$

We will discretize the problem with an isogeometric method, for which we will
assume that the domain is given by a parameterization of the form
$\Omega=\mathcal{F}(\hat{\Omega})$, where $\hat{\Omega}=[0,1]^{d}$ is the
parametric domain, and $\mathcal{F}$ is the isogeometric map to be detailed in
Section 3. We split the boundary of $\Omega$ in two parts, with mutually
disjoint interiors, denoted $\Gamma_{D}=\mathcal{F}(\hat{\Gamma}_{D})$ and
$\Gamma_{P}=\mathcal{F}(\hat{\Gamma}_{P})$, respectively corresponding to
Dirichlet and periodic boundary sides. Moreover, we assume that periodicity
occurs only in the last parametric direction, and split the periodic boundary
into two parts, $\Gamma_{P,1}$ and $\Gamma_{P,2}$, such that
$\Gamma_{P}=\Gamma_{P,1}\cup\Gamma_{P,2}$, where
$\Gamma_{P,i}=\mathcal{F}(\hat{\Gamma}_{P,i})$ for $i=1,2$ with
$\hat{\Gamma}_{P,1}=[0,1]^{d-1}\times\\{0\\}$ and
$\hat{\Gamma}_{P,2}=[0,1]^{d-1}\times\\{1\\}$. An illustration of the
isogeometric map and the split of the boundary is given in Figure 1.

$0$$1$$\hat{x}_{1}$$0$$1$$2$$x_{1}$$1$$\hat{x}_{2}$$1$$2$$x_{2}$$\hat{\Omega}$$\hat{\Gamma}_{P,1}$$\hat{\Gamma}_{P,2}$$\hat{\Gamma}_{D}$$\hat{\Gamma}_{D}$$\Omega$$\mathbf{n}_{1}$$\Gamma_{P,1}$$\Gamma_{D}$$\mathbf{n}_{2}$$\Gamma_{P,2}$$\Gamma_{D}$$\mathcal{F}$
Figure 1: Isogeometric map for the quarter of a bidimensional ring $\Omega$
with mixed Dirichlet and periodic boundaries.

Then, we impose the Dirichlet boundary condition

$u=g\qquad\text{on }\Gamma_{D}\times I,$

where we assume for simplicity that $g$ does not depend on the time $t$. To
impose the periodic boundary conditions, we introduce for scalar fields the
trace operator $\gamma:H^{1}(\Omega)\to H^{\frac{1}{2}}(\partial\Omega)$,
$\gamma:u\mapsto u|_{\partial\Omega}$, and for vector fields the normal trace
operator ${\gamma}_{n}:H(\text{div};\Omega)\to
H^{-\frac{1}{2}}(\partial\Omega)$,
$\gamma_{n}:\mathbf{v}\mapsto\mathbf{v}\cdot\mathbf{n}$, with $\mathbf{n}$ the
outgoing unit normal at $\partial\Omega$. When needed, we will restrict these
operators to a part of the boundary, simply denoting this restriction as
$\cdot|_{\Gamma_{D}}$ or $\cdot|_{\Gamma_{P,1}}$ and $\cdot|_{\Gamma_{P,2}}$.
With this notation, Dirichlet boundary condition reads like
$\gamma(u)|_{\Gamma_{D}}=g$, for every $t\in I$, and for simplicity we will
assume that $g$ does not depend on $t$, while periodic boundary conditions
read like, for every $t\in I$,

$\displaystyle\gamma(u)\big{|}_{\Gamma_{P,1}}=\gamma(u)\big{|}_{\Gamma_{P,2}},\quad$
$\displaystyle\quad\gamma(u_{t})\big{|}_{\Gamma_{P,1}}=\gamma(u_{t})\big{|}_{\Gamma_{P,2}},$
$\displaystyle{\gamma}_{n}(c\nabla u)\big{|}_{\Gamma_{P,1}}$
$\displaystyle={\gamma}_{n}(c\nabla u)\big{|}_{\Gamma_{P,2}}.$

In what follows, we make the assumption that $c$ is periodic, in the sense
that its pull-back $\hat{c}=c\circ\mathcal{F}$, is periodic in the last
parametric direction of $\hat{\Omega}$, that is, it satisfies periodicity
conditions on $\hat{\Gamma}_{P}$. Therefore, the Dirichlet and periodic
boundary conditions for the first order system (2) read as

$\displaystyle\phi=0\quad\text{on }\Gamma_{D}\times I,$
$\displaystyle{\gamma}_{n}(\mathbf{v})\big{|}_{\Gamma_{P,1}}={\gamma}_{n}(\mathbf{v})\big{|}_{\Gamma_{P,2}},\quad\gamma(\phi)\big{|}_{\Gamma_{P,1}}=\gamma(\phi)\big{|}_{\Gamma_{P,2}},\quad\text{
for all }t\in I.$

We remind that existence and uniqueness results for the solution $u$ are well
known, see [18, Section 7.2.2], while the regularity of the solution will
depend on the regularity of the initial conditions [18, Section 7.2.3].

###  2.2 Weak formulation and conservation of energy

Let us define the following Hilbert spaces over $\Omega\subset\mathbb{R}^{d}$.
$L^{2}(\Omega)$ is the usual Hilbert space of square integrable functions over
$\Omega$, endowed with the classical $L^{2}$-norm $\|\cdot\|_{L^{2}(\Omega)}$.
By $\mathbf{L}^{2}(\Omega)$ we denote its vectorial counterpart. The Hilbert
spaces $H^{k}(\Omega)$ denote the functions in $L^{2}(\Omega)$ such that their
$k$th-order derivatives also belong to $L^{2}(\Omega)$, and their vectorial
counterparts will be denoted by $\mathbf{H}^{k}(\Omega)$. We define

$\displaystyle\mathbf{H}(c,\text{div};\Omega)$
$\displaystyle\coloneqq\\{\mathbf{v}\in\mathbf{L}^{2}(\Omega):\text{div
}(c\mathbf{v})\in L^{2}(\Omega)\\},$
$\displaystyle\mathbf{H}_{P}(c,\text{div};\Omega)$
$\displaystyle\coloneqq\\{\mathbf{v}\in\mathbf{H}(c,\text{div};\Omega):{\gamma}_{n}(c\mathbf{v})|_{\Gamma_{P,1}}={\gamma}_{n}(c\mathbf{v})|_{\Gamma_{P,2}}\\},$

We will also make use of the space $\mathbf{H}^{k}(\textup{div};\Omega)$, the
space of functions in $\mathbf{H}^{k}(\Omega)$ such that their divergence
belongs to ${H}^{k}(\Omega)$. We can then write the weak formulation of (2),
that is: find $\mathbf{v}\in\mathbf{H}_{P}(c,\text{div};\Omega)$ and $\phi\in
L^{2}(\Omega)$ such that the following equations hold

(3a) $\displaystyle\int_{\Omega}\mathbf{v}_{t}\cdot\mathbf{w}\
\mathrm{d}\mathbf{x}$ $\displaystyle=-\int_{\Omega}\text{div}(c\mathbf{w})\
\phi\ \mathrm{d}\mathbf{x}\quad$
$\displaystyle\forall\mathbf{w}\in\mathbf{H}_{P}(c,\text{div};\Omega),$ (3b)
$\displaystyle\int_{\Omega}\phi_{t}\psi\ \mathrm{d}\mathbf{x}$
$\displaystyle=\int_{\Omega}\text{div}(c\mathbf{v})\ \psi\
\mathrm{d}\mathbf{x},\quad$ $\displaystyle\forall\psi\in L^{2}(\Omega).$

Notice that adding equation (3a) to equation (3b) and choosing the test
functions as $\mathbf{w}=\mathbf{v}$ and $\psi=\phi$, we obtain

$\int_{\Omega}\mathbf{v}_{t}\cdot\mathbf{v}\mathrm{d}\mathbf{x}+\int_{\Omega}\phi_{t}\phi\mathrm{d}\mathbf{x}=0,$

and defining the quantity of total energy as

(4)
$E(t)\coloneqq\frac{1}{2}\int_{\Omega}|\mathbf{v}|^{2}+|\phi|^{2}\mathrm{d}\mathbf{x},$

it is obviously seen that $\frac{\mathrm{d}E(t)}{\mathrm{d}t}=0$, that is, the
energy is preserved.

Building upon the work conducted by Kraus et al. in the field of
magnetohydrodynamics [23, 21], we aim to introduce and analyze, for this
simplified model problem, a numerical method that preserves the total energy
of the system for long-term simulations, and for which (3a) is solved in a
weak sense, while (3b) is solved in strong form. To achieve this purpose, we
introduce a generic linear operator

$\mathcal{T}:L^{2}(\Omega)\to L^{2}(\Omega),$

and replace the second equation (3b) by

(5) $\phi_{t}=\mathcal{T}(\text{div}(c\mathbf{v})).$

Note that we recover (3b) if $\mathcal{T}$ is the identity operator. This
change in the equation affects the conservation of the total energy of the
system. In order to retrieve this conservation we replace (3a) with

$\int_{\Omega}\mathbf{v}_{t}\cdot\mathbf{w}\
\mathrm{d}\mathbf{x}=-\int_{\Omega}\mathcal{T}(\text{div}(c\mathbf{w}))\phi\
\mathrm{d}\mathbf{x}.$

The modified problem that we obtain, with mixed homogeneous Dirichlet and
periodic boundary conditions, reads: find
$\mathbf{v}\in\mathbf{H}_{P}(c,\text{div};\Omega)$ and $\phi\in L^{2}(\Omega)$
such that the following equations hold

(6a) $\int_{\Omega}\mathbf{v}_{t}\cdot\mathbf{w}\
\mathrm{d}\mathbf{x}=-\int_{\Omega}\mathcal{T}(\text{div}(c\mathbf{w}))\phi\
\mathrm{d}\mathbf{x},\quad\forall\mathbf{w}\in\mathbf{H}_{P}(c,\text{div};\Omega),$
(6b) $\phi_{t}=\mathcal{T}(\text{div}(c\mathbf{v})).$

##  3 B-splines and IGA framework

In this section we recall the definition of B-splines in the univariate and
multivariate context. We also introduce the isogeometric spaces, together with
their regularity assumptions.

###  3.1 Univariate B-splines

For the definition of B-splines we follow the notation and the guidelines of
[5] and [26]. Let us introduce a knot vector
$\Xi=\\{\xi_{1},\dots,\xi_{n+p+1}\\}$, where $\xi_{i}\in\mathbb{R}$ is the
$i$-th knot, $p$ is the polynomial degree and $n$ is the number of basis
functions. We assume that
$\xi_{1}\leq\xi_{2}\leq\dots\leq\xi_{n+p}\leq\xi_{n+p+1},$ and without loss of
generality $\xi_{1}=0,\ \xi_{n+p+1}=1$. We also admit at most $p+1$ repeated
knots. We introduce the B-spline functions $\\{\hat{B}_{i,p}\\}_{i=1}^{n}$ of
degree $p$ over the knot vector $\Xi$ following the Cox-DeBoor recursive
formula [12, Chapter IX]. We obtain a set of $n$ B-splines with the following
properties: non-negativity, partition of unity and local support of each basis
function, see [5]. The univariate B-spline space of degree $p$ over the knot
vector $\Xi$ is:

$\hat{S}_{p}(\Xi)\coloneqq\text{span}\\{\hat{B}_{i,p}:i=1,\dots,n\\}.$

We introduce also the vector $\mathbf{Z}=\\{\zeta_{1},\dots,\zeta_{z}\\}$, of
knots without repetitions, which are also called breakpoints, and denote with
$m_{j}$ the multiplicity of $\zeta_{j}$, such that
$\sum_{j=1}^{z}m_{j}=n+p+1$, and $\zeta_{1}=\xi_{1}=\dots=\xi_{m_{1}}$,
$\zeta_{2}=\xi_{m_{1}+1}=\dots=\xi_{m_{1}+m_{2}}$ and so on. Notice that

$\Xi=\\{\underbrace{\zeta_{1},\dots,\zeta_{1}}_{m_{1}\text{
times}},\underbrace{\zeta_{2},\dots,\zeta_{2}}_{m_{2}\text{
times}},\dots,\underbrace{\zeta_{z},\dots,\zeta_{z}}_{m_{z}\text{ times}}\\},$

with each $\zeta_{j}$ repeated $m_{j}$ times, and $1\leq m_{j}\leq p+1$ for
all internal knots. The breakpoints form a partition of the unit interval, and
we say the partition is locally quasi-uniform if there exists a constant
$\beta\geq 1$ independent of $z$ such that

$\beta^{-1}\leq\dfrac{\zeta_{j+1}-\zeta_{j}}{\zeta_{j}-\zeta_{j-1}}\leq\beta,\quad\forall
j=2,\dots,z-1.$

Assuming that the multiplicity of internal knots is at most $p$, we can write
the derivative of a B-spline as follows

(7)
$\frac{\partial\hat{B}_{i,p}}{\partial\xi}(\xi)=\hat{D}_{i-1,p-1}(\xi)-\hat{D}_{i,p-1}(\xi),$

with the Curry-Schoenberg spline basis:

(8)
$\hat{D}_{i,p-1}(\xi)=\frac{p}{\xi_{i+p+1}-\xi_{i+1}}\hat{B}_{i+1,p-1}(\xi),\quad\text{for
}i=1,\dots,n-1.$

where we assumed that $\hat{D}_{1,p-1}(\xi)=\hat{D}_{n+1,p-1}(\xi)=0$. Note
that the derivative belongs to the spline space $\hat{S}_{p-1}(\Xi^{\prime})$
where $\Xi^{\prime}=\\{\xi_{2},\dots,\xi_{n+p}\\}$.

In order to handle periodic boundary conditions, we follow the construction of
periodic B-spline spaces, see for instance [29]. Consider a uniformly spaced
closed knot vector $\Xi$ of degree $p$, that is $\xi_{1}<\dots<\xi_{p+1}=0$
and $1=\xi_{n+1}<\dots<\xi_{n+p+1}$ with at least $p$ internal knots. We can
construct a periodic basis on such knot vector, by identifying the basis
functions of the B-spline space $\hat{S}_{p}(\Xi)$ in this way:

$\begin{cases}\hat{B}_{i,p}^{Per}\coloneqq\hat{B}_{i,p}+\hat{B}_{n-p+i,p},&\text{for
}i=1,\dots,p;\\\
\hat{B}_{i,p}^{Per}=\hat{B}_{i,p},&\text{otherwise.}\end{cases}$

Figure 2: B-spline basis of degree 3 in the periodic case. In color the first
three B-spline basis. The left hand side matches with regularity $C^{2}$ with
their right hand side.

Here we are gluing together the first $p$ basis functions with the last $p$,
this way we get continuity of derivatives up to order $p-1$. In Figure 2 is
shown the behavior of this periodic basis at the boundary. We can introduce
the periodic B-spline space with highest regularity at the boundary as:

$\hat{S}_{p}^{Per}(\Xi)\coloneqq\text{span}\\{\hat{B}_{i,p}^{Per}:i=1,\dots,n-p\\}.$

It is easy to see that the dimension of the spline space
$\hat{S}_{p}^{Per}(\Xi)$ is $n-p$, which is the same of
$\hat{S}_{p-1}^{Per}(\Xi^{\prime})$. Finally we note that equation (7) holds
for periodic B-splines too. For more information and properties on B-splines,
we refer the reader to [12].

###  3.2 Multivariate B-splines

Multivariate B-splines are introduced and defined by tensor product, starting
from univariate B-splines on each spatial direction. Let $d$ be the space
dimension, usually $d=2,3$, and assume $n_{l}\in\mathbb{N}$ the number of
univariate basis functions in direction $l$, $p_{l}\in\mathbb{N}$ is the
degree, while $\Xi_{l}=\\{\xi_{l,1},\dots,\xi_{l,n_{l}+p_{l}+1}\\}$ and
$\mathbf{Z}_{l}=\\{\zeta_{l,1},\dots,\zeta_{l,z_{l}}\\}$ are respectively the
knots and breakpoints in direction $l$. We also set the polynomial degree
vector $\mathbf{p}=(p_{1},\dots,p_{d})$ and
$\mathbf{\Xi}=\\{\Xi_{1},\dots,\Xi_{d}\\}$. Finally we set
$\mathbf{J}=\\{\mathbf{j}=(j_{1},\dots,j_{d})\subset\mathbb{N}^{d}:1\leq
j_{l}\leq n_{l}\\}$. The set of multivariate B-splines basis functions of
degree vector $\mathbf{p}$ is

$\\{\hat{B}_{\mathbf{i,p}}(\boldsymbol{\xi})=\hat{B}_{i_{1},p_{1}}(\xi_{1})\dots\hat{B}_{i_{d},p_{d}}(\xi_{d}),\text{for
}\mathbf{i}\in\mathbf{J}\\},$

and the B-spline multivariate space of degree vector $\mathbf{p}$ over
$\mathbf{\Xi}$ is

$\hat{S}_{\mathbf{p}}(\mathbf{\Xi})\coloneqq\text{span}\\{\hat{B}_{\mathbf{i,p}}(\boldsymbol{\xi}):\mathbf{i}\in\mathbf{J}\\}=\otimes_{l=1}^{d}\hat{S}_{p_{l}}(\Xi_{l}).$

A similar construction also applies for multivariate periodic B-splines. It is
also possible to combine tensor product of standard B-splines in the first
directions with periodic ones for the last direction, which is the case of the
mixed boundary condition setting for our model problem.

$\mathbf{Q}$$\mathbf{K}$$\mathcal{F}$ Figure 3: Mesh $\widehat{\mathcal{M}}$
in the parametric domain, and its image $\mathcal{M}$ on the physical domain.

###  3.3 Isogeometric analysis framework

The isogeometric map $\mathcal{F}:\hat{\Omega}\to\Omega$, is a
parameterization of the geometry of the physical domain, based on B-splines
(or more often in their rational counterpart of NURBS), usually indicated by

$\mathcal{F}\coloneqq\sum_{\mathbf{i}\in\mathbf{J}}\mathbf{c_{i}}\hat{B}_{\mathbf{i,p}}.$

For the analysis of the discrete problem, we need to introduce some
assumptions on $\mathcal{F}$. First, we introduce the parametric Bézier mesh
$\widehat{\mathcal{M}}$, which is

$\widehat{\mathcal{M}}\coloneqq\\{\mathbf{Q_{j}}=I_{1,j_{1}}\times\dots\times
I_{d,j_{d}}:I_{l,j_{l}}=(\zeta_{l,j_{l}},\zeta_{l,j_{l}+1}),\text{ for }1\leq
j_{l}\leq z_{l}-1\\}.$

Given an element $\mathbf{Q}\in\widehat{\mathcal{M}}$, we set
$h_{\mathbf{Q}}=\textup{diam}(\mathbf{Q}_{\mathbf{j}})$, while
$h=\max\\{h_{\mathbf{Q}},\ \mathbf{Q}\in\widehat{\mathcal{M}}\\}$. Moreover,
given $\mathbf{D}\subset\hat{\Omega}$, we denote by $\widetilde{{\mathbf{D}}}$
its support extension, that is the interior of the union of the supports of
basis functions whose support intersects ${\mathbf{D}}$.

The Bézier mesh is defined as the image of elements in $\widehat{\mathcal{M}}$
through $\mathcal{F}$:

$\mathcal{M}\coloneqq\\{\mathbf{K}\subset\Omega:\mathbf{K}=\mathcal{F}(\mathbf{Q}),\mathbf{Q}\in\widehat{\mathcal{M}}\\},$

see Figure 3. We have
$\widetilde{\mathcal{F}(\mathbf{D})}=\mathcal{F}(\widetilde{\mathbf{D}})$, for
every subset $\mathbf{D}\subset\hat{\Omega}$, denoting the support extension
in the physical domain. The first assumption is the following.

######  Assumption 1.

We assume that $\mathcal{F}$ is a bi-Lipschitz homeomorphism. Moreover,
$\mathcal{F}|_{\overline{\mathbf{Q}}}$ belongs to
$\mathcal{C}^{\infty}(\overline{\mathbf{Q}})$ for all
$\mathbf{Q}\in\widehat{\mathcal{M}}$, where $\overline{\mathbf{Q}}$ denotes
the closure of $\mathbf{Q}$, and $\mathcal{F}^{-1}|_{\overline{\mathbf{K}}}$
belongs to $\mathcal{C}^{\infty}(\overline{\mathbf{K}})$, for all
$\mathbf{K}\in\mathcal{M}$.

This prevents the existence of singularities and self-intersections in the
parameterization $\mathcal{F}$. The second assumption simplifies the dealing
with the boundary.

######  Assumption 2.

The boundary region $\Gamma_{D}\subset\partial\Omega$ is the union of full
faces of the boundary. More precisely,
$\Gamma_{D}=\mathcal{F}(\hat{\Gamma}_{D})$, with $\hat{\Gamma}_{D}$ a
collection of full faces of the parametric domain $\hat{\Omega}$. Finally, if
$\Gamma_{P}\neq\emptyset$, we assume the periodic boundary to be:

$\Gamma_{P}=\mathcal{F}\big{(}[0,1]^{d-1}\times\\{0\\}\big{)}\cup\mathcal{F}\big{(}[0,1]^{d-1}\times\\{1\\}\big{)}.$

Finally, we assume local quasi-uniformity of the univariate partitions, thus
the parametric Bézier mesh is shape regular, that is the ratio between the
smallest edge of $\mathbf{Q}\in\widehat{\mathcal{M}}$ and its diameter
$h_{\mathbf{Q}}$ is bounded uniformly with respect to $h$ and $\mathbf{Q}$.
Analogously, the Bézier mesh is shape regular thanks to Assumption 1.

##  4 Discretization

In this section, our focus is on introducing the discretization of problem
(6). We approximate $\mathcal{T}$ with quasi-interpolant projections that
commute with the divergence operator. The Crank-Nicolson method is employed
for temporal discretization, chosen for its conservativity. We emphasize that
our objective is to study the approximation properties of the proposed method,
both theoretically and numerically. To achieve this, we intentionally maintain
a highly generic notation. In Section 4.1, we present the spatial semi-
discretization through the construction of isogeometric discrete spaces and
the associated projections. In Section 4.2 we discuss a particularization of
such projections, which are quasi-interpolant operators based on point
evaluation of the functions to be projected. This quasi-interpolant operators
will be subjected to numerical study and testing within the proposed
discretization framework. Finally, in Section 4.3 we recall the Crank-Nicolson
semi discretization in time and the fully discrete problem is presented.

###  4.1 Discretization in space with commutative projections

Here we introduce the discrete spaces for the mixed formulation given by
equations (6a) and (6b). We present the case $d=2$, while the case $d=3$ is
completely analogous.

####  4.1.1 Discrete spline spaces

Let us standardize the notation following [4], and start with the case of
Dirichlet boundary conditions, that is $\Gamma_{P}=\emptyset$, we define the
spaces:

${X^{1}\coloneqq\mathbf{H}(c,\text{div};\Omega)},\quad X^{2}\coloneqq
L^{2}(\Omega),\quad{\hat{X}^{1}\coloneqq\mathbf{H}(\hat{c},\text{div};\hat{\Omega})},\quad\hat{X}^{2}\coloneqq
L^{2}(\hat{\Omega}).$

Thanks to Assumption 1, which states that both $\mathcal{F}$ and its inverse
are smooth, we can define the pull-backs that relate these spaces as (see [20,
Section 2.2])

$\displaystyle\iota^{1}(\mathbf{f})$
$\displaystyle\coloneqq\det(D\mathcal{F})(D\mathcal{F})^{-1}(\mathbf{f}\circ\mathcal{F}),$
$\displaystyle\mathbf{f}\in X^{1},$ $\displaystyle\iota^{2}(f)$
$\displaystyle\coloneqq\det(D\mathcal{F})(f\circ\mathcal{F}),$ $\displaystyle
f\in X^{2},$

where $D\mathcal{F}$ is the Jacobian matrix of $\mathcal{F}$. Then, due to the
divergence preserving property of the map $\iota^{1}$, see [28, Section 3.9],
we have $\text{div}\circ\iota^{1}=\iota^{2}\circ\text{div}$. Fixed a
polynomial degree vector $\mathbf{p}=(p_{1},p_{2})$ and
$\mathbf{\Xi}=(\Xi_{1},\Xi_{2})$, we define the discrete spaces on the
parametric domain as:

$\displaystyle\hat{X}^{1}_{h}$
$\displaystyle\coloneqq\hat{S}_{p_{1},p_{2}-1}(\Xi_{1},\Xi_{2}^{\prime})\times\hat{S}_{p_{1}-1,p_{2}}(\Xi_{1}^{\prime},\Xi_{2}),$
$\displaystyle\hat{X}^{2}_{h}$
$\displaystyle\coloneqq\hat{S}_{p_{1}-1,p_{2}-1}(\Xi_{1}^{\prime},\Xi_{2}^{\prime}).$

The choice of the bases follows from [4, Section 5.2] and [31, Section 4],
that is

$\hat{X}^{1}_{h}=span\\{\mathcal{J}_{1}\cup\mathcal{J}_{2}\\},$

where we set

$\displaystyle\mathcal{J}_{1}$
$\displaystyle=\\{\hat{B}_{i_{1},p_{1}}(\xi_{1})\hat{D}_{i_{2},p_{2}-1}(\xi_{2})\mathbf{e_{1}}:\
1\leq i_{1}\leq n_{1},\ 1\leq i_{2}\leq n_{2}-1\\},$
$\displaystyle\mathcal{J}_{2}$
$\displaystyle=\\{\hat{D}_{i_{1},p_{1}-1}(\xi_{1})\hat{B}_{i_{2},p_{2}}(\xi_{2})\mathbf{e_{2}}:\
1\leq i_{1}\leq n_{1}-1,\ 1\leq i_{2}\leq n_{2}\\},$

and $\\{\mathbf{e_{1}},\mathbf{e_{2}}\\}$ is the canonical basis of
$\mathbb{R}^{2}$. As regards the second discrete space, we set

$\hat{X}^{2}_{h}=span\\{\hat{D}_{i_{1},p_{1}-1}(\xi_{1})\hat{D}_{i_{2},p_{2}-1}(\xi_{2}):\
1\leq i_{l}\leq n_{l}-1,\ l=1,2\\}.$

The discrete spaces on the physical domain $\Omega$ can be defined from the
spaces $\hat{X}^{1}_{h}$ and $\hat{X}^{2}_{h}$ by push-forward, that is, the
inverse of the transformations $\iota^{1}$ and $\iota^{2}$, that commute with
the divergence operator. We have the following definitions:

$\displaystyle
X^{1}_{h}\coloneqq\\{\mathbf{f_{h}}:\iota^{1}(\mathbf{f}_{h})\in\hat{X}^{1}_{h}\\},$
$\displaystyle
X^{2}_{h}\coloneqq\\{f_{h}:\iota^{2}(f_{h})\in\hat{X}^{2}_{h}\\}.$

For later use, we introduce a suitable notation for the basis of these
discrete spaces. We denote by $\\{\mathbf{b}_{i,h}\\}_{i=1}^{N}$ and
$\\{\hat{\mathbf{b}}_{i,h}\\}_{i=1}^{{N}}$ the sets of basis functions of
$X^{1}_{h}$ and $\hat{X}^{1}_{h}$ respectively, reordered with lexicographic
ordering, such that $\iota^{1}(\mathbf{b}_{i,h})=\hat{\mathbf{b}}_{i,h}$ for
any $i=1,\dots,{N},$ where $N=n_{1}(n_{2}-1)+(n_{1}-1)n_{2}$. Analogously we
can introduce the notations $\\{b_{i,h}\\}_{i=1}^{M}$ and
$\\{\hat{b}_{i,h}\\}_{i=1}^{M}$ for the basis functions of the spaces
$X^{2}_{h}$ and $\hat{X}^{2}_{h}$ respectively, such that, for $i=1,\dots,M$
with $M=(n_{1}-1)(n_{2}-1)$ we have $\iota^{2}(b_{i,h})=\hat{b}_{i,h}$.

####  4.1.2 Univariate projections

In order to discretize equations (6a) and (6b), We follow the idea in [21] and
in [24], and approximate the operator $\mathcal{T}$ with a projector into the
discrete spline space. We first introduce commutative projectors in the
univariate case, which can be defined by a dual basis, i.e.
$\hat{\pi}_{p,\Xi}:L^{2}(0,1)\longrightarrow\hat{S}_{p}(\Xi)$ such that

(9) $\hat{\pi}_{p,\Xi}(f)=\sum_{i=1}^{n}\lambda_{i,p}(f)\hat{B}_{i,p}$

where $\lambda_{i,p}$ are a set of dual functionals verifying
$\lambda_{i,p}(\hat{B}_{j,p})=\delta_{ij}$, where $\delta_{ij}$ is the
Kronecker delta. Notice that, with this property, the operator is a
projection, that is $\hat{\pi}_{p,\Xi}(s)=s$, for all $s\in\hat{S}_{p}(\Xi)$.

We also recall the construction of univariate projections that commute with
the derivative as in [10, Section 3.1.2]. In order to obtain the one-
dimensional commuting diagram, we want the commutative projection
$\hat{\pi}^{c}_{p-1,\Xi^{\prime}}:L^{2}(0,1)\to S_{p-1}(\Xi^{\prime})$ to
satisfy:

$\hat{\pi}_{p-1,\Xi^{\prime}}^{c}\frac{\partial}{\partial\xi}f=\frac{\partial}{\partial\xi}\hat{\pi}_{p,\Xi}f,$

for all $f$ in $H^{1}(0,1)$. In order to satisfy the previous equation, given
$\hat{\pi}_{p,\Xi}$ as in (9), its commutative projection is defined as

(10)
$\hat{\pi}_{p-1,\Xi^{\prime}}^{c}g(\xi)\coloneqq\frac{\partial}{\partial\xi}\hat{\pi}_{p,\Xi}\int_{0}^{\xi}g(s)\mathrm{d}s,$

for all functions $g$ such that $G(\xi)=\int_{0}^{\xi}g(s)\mathrm{d}s$ is in
the domain of $\hat{\pi}_{p,\Xi}$.

To apply periodic boundary conditions we need to define the commutative
projector in a different way, to ensure that integrating a periodic function
in $L^{2}(0,1)$, we obtain a periodic function in $H^{1}(0,1)$. Hence, given a
periodic projection $\hat{\pi}_{p,\Xi}^{Per}$, we define its commutative one
as follows

(11)
$\hat{\pi}_{p-1,\Xi^{\prime}}^{c,Per}g(\xi)\coloneqq\bigg{(}\frac{\partial}{\partial\xi}\hat{\pi}_{p,\Xi}^{Per}\int_{0}^{\xi}\tilde{g}(s)\mathrm{d}s\bigg{)}\oplus\overline{g},$

for all functions $g=\tilde{g}\oplus\overline{g}\in L^{2}(0,1)$ such that
$\overline{g}\coloneqq\int_{0}^{1}g(s)\mathrm{d}s\in\mathbb{R}$ and
$\tilde{g}\in L^{2}_{0}(0,1)$, the space of functions in $L^{2}(0,1)$ with
zero average. Notice that $G(\xi)=\int_{0}^{\xi}{\tilde{g}}(s)\mathrm{d}s$ is
in the domain of $\hat{\pi}_{p,\Xi}^{Per}$, and again we have

(12)
$\hat{\pi}_{p-1,\Xi^{\prime}}^{c,Per}\frac{\partial}{\partial\xi}f=\frac{\partial}{\partial\xi}\hat{\pi}_{p,\Xi}^{Per}f,$

for all $f$ in the domain of definition of $\hat{\pi}_{p,\Xi}^{Per}$. It is
easy to see that the commutative projections defined in this way preserve
splines.

####  4.1.3 Multivariate construction

The univariate projection operators introduced can be extended to the
multidimensional case by tensor product constructions. Given
$\hat{\Omega}=[0,1]^{d}\subset\mathbb{R}^{d}$, for $i=1,2,\dots,d$, let us
denote with $\hat{\pi}_{p_{i},\Xi_{i}}$ a generic univariate projection as in
(9). We can define a multivariate projection as

(13)
$\hat{\Pi}_{\mathbf{p,\Xi}}\coloneqq\hat{\pi}_{p_{1},\Xi_{1}}\otimes\hat{\pi}_{p_{2},\Xi_{2}}\otimes\dots\otimes\hat{\pi}_{p_{d},\Xi_{d}}.$

It is important to note that it can be expressed as

(14)
$\hat{\Pi}_{\mathbf{p},\mathbf{\Xi}}(f)=\sum_{\mathbf{i}\in\mathbf{J}}\lambda_{\mathbf{i,p}}(f)\hat{B}_{\mathbf{i,p}},$

where $\mathbf{p}$, $\mathbf{\Xi}$ and $\mathbf{J}$ are given as in the
multivariate B-spline construction, and each dual functional is defined from
the univariate dual basis by the expression

(15)
$\lambda_{\mathbf{i,p}}=\lambda_{i_{1},p_{1}}\otimes\lambda_{i_{2},p_{2}}\otimes\dots\otimes\lambda_{i_{d},p_{d}}.$

It is easy to see that the same constructions hold for the multivariate
periodic case, or for combinations with periodic boundary conditions in only
one parametric coordinate.

At this point we have all the ingredients to define projections on the spaces
$\hat{X}^{1}_{h}$ and $\hat{X}^{2}_{h}$. The choice of the projection follows
from the definition of the discrete spaces. More precisely we set

(16) $\displaystyle\hat{\Pi}^{1}$
$\displaystyle=(\hat{\pi}_{p_{1}}\otimes\hat{\pi}_{p_{2}-1}^{c})\times(\hat{\pi}_{p_{1}-1}^{c}\otimes\hat{\pi}_{p_{2}}),$
(17) $\displaystyle\hat{\Pi}^{2}$
$\displaystyle=\hat{\pi}_{p_{1}-1}^{c}\otimes\hat{\pi}_{p_{2}-1}^{c}.$

where for simplicity of notation we have omitted the knot vector from the
subscript of the univariate projectors. These projectors satisfy the spline
preserving properties, see [4, Lemma 5.3].

The projectors into the discrete spline spaces in the physical domain are
defined from the ones in the parametric domain (16) and (17), and the
corresponding pull-backs $\iota^{1}$ and $\iota^{2}$, in such a way they are
uniquely characterized by the equations

(18)
$\displaystyle\iota^{1}(\Pi^{1}\mathbf{f})=\hat{\Pi}^{1}(\iota^{1}(\mathbf{f})),$
$\displaystyle\iota^{2}(\Pi^{2}f)=\hat{\Pi}^{2}(\iota^{2}(f)).$

These projectors satisfy the commutativity property with the divergence
operator, as stated in the following lemma:

######  Lemma 1.

Given $\hat{c}$ and $c$ the coefficients of spaces $\hat{X}^{1}$ and $X^{1}$,
the following equations hold

(19) $\displaystyle\textup{div }(\hat{\Pi}^{1}(\hat{c}\textbf{f}))$
$\displaystyle=\hat{\Pi}^{2}(\textup{div
}(\hat{c}\textbf{f})),\quad\forall\textbf{f}\in\hat{X}^{1},$ (20)
$\displaystyle\textup{div }(\Pi^{1}(c\textbf{f}))$
$\displaystyle=\Pi^{2}(\textup{div }(c\textbf{f})),\quad\forall\textbf{f}\in
X^{1}.$

The proof of (19) is given in [4, Lemma 5.5]. Equation (20) is an immediate
consequence of the definitions, together with the commutativity property (19).
We conclude this paragraph with the following remark on the periodic case.

######  Remark 1.

If $\Gamma_{P}\neq\emptyset$, we proceed analogously introducing the discrete
spaces

(21) $\displaystyle\hat{X}^{1}_{h}$
$\displaystyle\coloneqq\Big{(}\hat{S}_{p_{1}}(\Xi_{1})\otimes\hat{S}_{p_{2}-1}^{Per}(\Xi_{2}^{\prime})\Big{)}\times\Big{(}\hat{S}_{p_{1}-1}(\Xi_{1}^{\prime})\otimes\hat{S}_{p_{2}}^{Per}(\Xi_{2})\Big{)},$
(22) $\displaystyle\hat{X}^{2}_{h}$
$\displaystyle\coloneqq\hat{S}_{p_{1}-1}(\Xi_{1}^{\prime})\otimes\hat{S}_{p_{2}-1}^{Per}(\Xi_{2}^{\prime}),$

and considering the projectors:

(23) $\displaystyle\hat{\Pi}^{1}$
$\displaystyle=(\hat{\pi}_{p_{1}}\otimes\hat{\pi}_{p_{2}-1}^{c,Per})\times(\hat{\pi}_{p_{1}-1}^{c}\otimes\hat{\pi}_{p_{2}}^{Per}),$
(24) $\displaystyle\hat{\Pi}^{2}$
$\displaystyle=\hat{\pi}_{p_{1}-1}^{c}\otimes\hat{\pi}_{p_{2}-1}^{c,Per}.$

####  4.1.4 Semi-discretization in space

We propose to take as a linear operator $\mathcal{T}$ in (6a) and (6b) the
generic tensor product commutative projection $\Pi^{2}$. The semi-discrete
problem in space reads, find $\mathbf{v}_{h}\in X^{1}_{h}$ and $\phi_{h}\in
X^{2}_{h}$, such that the following equations hold:

(25a)
$\displaystyle\int_{\Omega}\left(\mathbf{v}_{h}\right)_{t}\cdot\mathbf{w}_{h}\
\mathrm{d}\mathbf{x}$
$\displaystyle=-\int_{\Omega}\Pi^{2}(\text{div}(c\mathbf{w}_{h}))\ \phi_{h}\
\mathrm{d}\mathbf{x},\quad$
$\displaystyle\forall\mathbf{w}_{h}\in{X}^{1}_{h},$ (25b)
$\displaystyle\left(\phi_{h}\right)_{t}$
$\displaystyle=\Pi^{2}(\text{div}(c\mathbf{v}_{h})).$

In next section we will exploit a particularization of such discretization,
based on quasi-interpolant projections.

###  4.2 Quasi interpolant projections

$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.5$$0$$0.5$$1$$f(\xi_{i+1})$$f\left(\dfrac{\xi_{i+1}+\xi_{i+2}}{2}\right)$$f(\xi_{i+2})$$x$$f(x)=sin(2\pi
x)$$f(x)$ on breakpoints $f(x)$ on midpoints Figure 4: Example of pointwise
evaluation of $f(x)=sin(2\pi x)$ over breakpoints and midpoints for projection
with $\hat{\pi}_{2,\Xi}$ where $\Xi=\\{0,0,0,0.2,0.4,0.6,0.8,1,1,1\\}$. The
three highlighted points are used for the computation of $\lambda_{i,2}(f)$,
for $i=4$, corresponding to the B-spline with support $[0.2,0.8]$.

Here we focus on a particular projection whose construction is given in detail
in [25]. We refer to such projection using the notation $\llm[p]$. We deal
with open knot vectors with non-repeating internal knots, i.e.
$\Xi=\\{\xi_{1},\dots,\xi_{n+p+1}\\}$ such that the vector
$\mathbf{Z}=\\{\zeta_{1},\dots,\zeta_{z}\\}$ of the breakpoints has
multiplicities $m_{1}=m_{z}=p+1$ and $m_{j}=1$ for $j=2,\dots,z-1$. In order
to find an approximation $\llm[p](f)$ of $f:[0,1]\to\mathbb{R}$, we need to
define a dual basis $\lambda_{i,p}$ and compute $\lambda_{i,p}(f)$ for
$i=1,\dots,n$, and the idea is the following:

  1. 1.

Choose $\mathcal{I}=[\xi_{\mu},\xi_{\nu}]\subset[\xi_{p+1},\xi_{n}]$ such that
$\mathcal{I}\cap[\xi_{i},\xi_{i+p+1}]\neq\emptyset$. Notice that
$[\xi_{i},\xi_{i+p+1}]$ is the support of $\hat{B}_{i,p}$, and denote by
$f_{\mathcal{I}}$ the restriction of $f$ to $\mathcal{I}$.

  2. 2.

Choose a local approximation method $\mathcal{P}_{loc}$ for $f_{\mathcal{I}}$,
such that it reproduces splines of degree up to $p$. This approximation is
written analogously to (9) in this way:
$\mathcal{P}_{loc}(f_{\mathcal{I}})=\sum_{j=\mu-p}^{\nu-1}c_{j}\hat{B}_{j,p}$,
and we have $\mu-p\leq i\leq\nu-1$ since the support of $\hat{B}_{i,p}$
intersects $\mathcal{I}$.

  3. 3.

Finally set $\lambda_{i,p}(f)=c_{i}$.

The idea is to use polynomial interpolation as local approximation method. We
exploit the construction for spline spaces of degree $2$ and $3$ given in
[25]. Let us fix now $p=2$, for a given index $i=1,\dots,n$, we have
$[\xi_{i+1},\xi_{i+2}]\subset[\xi_{i},\xi_{i+3}]=supp(\hat{B}_{i,2})$, but
$\xi_{i+1}<\xi_{i+2}$ when $i\neq 1,n$. This means that for $i=2,\dots,n-1$ we
can choose $\mathcal{I}=[\xi_{i+1},\xi_{i+2}]$, and perform three points
interpolation on $\xi_{i+1}<(\xi_{i+1}+\xi_{i+2})/2<\xi_{i+2}$. Notice that we
are interpolating over two consecutive breakpoints and their midpoint in a
uniform partition. Due to this particular choice of the interpolation points,
we can explicitly write the coefficients of $\llm[2]$, that are:

(26)
$\lambda_{i,2}(f)=-\frac{1}{2}f(\xi_{i+1})+2f\Big{(}\frac{\xi_{i+1}+\xi_{i+2}}{2}\Big{)}-\frac{1}{2}f(\xi_{i+2}),\quad\text{for
}i=2,\dots,n-1.$

We remark this expression is valid whenever $\xi_{i+1}<\xi_{i+2}$, which is
not the case for $i=1$ and $i=n$ in the case of an open knot vector, as used
for Dirichlet boundary conditions. In this situation we want the operator to
be interpolant at the boundary and hence we fix $\lambda_{1,2}=f(0)$ and
$\lambda_{n,2}=f(1)$. As regards closed knot vectors for periodic boundary
conditions, the computation is easier, since we can use (26) for all indices
$i=1,\dots,n$. Notice that we need the pointwise evaluation of the function
$f$ we want to approximate, among the breakpoints and their midpoints, as in
the example given in Figure 4.

We recall now the explicit formulae in the case $p=3$, that is the quasi-
interpolant $\llm[3]$, and refer for details to [25]. Here, as local
approximation method, we use five points interpolation over the knots and
midpoints of $\mathcal{I}=[\xi_{i+1},\xi_{i+3}]$. Instead of (26), we can
write now

(27)
$\lambda_{i,3}(f)=\frac{1}{6}f(\xi_{i+1})-\frac{4}{3}f\Big{(}\frac{\xi_{i+1}+\xi_{i+2}}{2}\Big{)}+\frac{10}{3}f(\xi_{i+2})-\frac{4}{3}f\Big{(}\frac{\xi_{i+2}+\xi_{i+3}}{2}\Big{)}+\frac{1}{6}f(\eta_{i+3}).$

Again for open knot vectors, hence Dirichlet boundary conditions, we have the
particular cases for $i=1,2,n-1,n,$ that we recall here:

$\displaystyle\lambda_{1,3}(f)$ $\displaystyle=f(0),$
$\displaystyle\lambda_{2,3}(f)$
$\displaystyle=-\frac{5}{18}f(\xi_{4})+\frac{20}{9}f\Big{(}\frac{\xi_{4}+\xi_{5}}{2}\Big{)}-\frac{4}{3}f(\xi_{5})+\frac{4}{9}f\Big{(}\frac{\xi_{5}+\xi_{6}}{2}\Big{)}-\frac{1}{18}f(\xi_{6}),$
$\displaystyle\lambda_{n-1,3}(f)$
$\displaystyle=-\frac{1}{18}f(\xi_{n-1})+\frac{4}{9}f\Big{(}\frac{\xi_{n-1}+\xi_{n}}{2}\Big{)}-\frac{4}{3}f(\xi_{n})+$
$\displaystyle\quad+\frac{20}{9}f\Big{(}\frac{\xi_{n}+\xi_{n+1}}{2}\Big{)}-\frac{5}{18}f(\xi_{n+1}),$
$\displaystyle\lambda_{n,3}(f)$ $\displaystyle=f(1),$

From the computational point of view, we remember this technique requires
evaluation of $f$ over the breakpoints and midpoints. Notice that for
projections in periodic spline spaces, it is sufficient to use (26) or (27).
Finally, for approximation estimates by quasi-interpolant projections based on
point evaluation of functions, we refer to [27, Section 5].

####  4.2.1 Commutative quasi-interpolant projections

In order to compute a projection $\llm[p]^{c}$ that commutes with $\llm[p]$,
we follow the construction in (10). Since we have to apply $\llm[p]$ to the
integral function $F$, we need its evaluation over all the breakpoints and
midpoints. For this reason, instead of using Gaussian quadrature, we prefer to
use Cavalieri-Simpson composite quadrature formulae. We want to consider both
breakpoints and midpoints in a unified notation. Hence we introduce
$\boldsymbol{\eta}=\\{\eta_{1},\dots,\eta_{2z-1}\\}$ such that:

$\begin{cases}\eta_{2i-1}=\zeta_{i},\quad&\text{for }i=1,\dots,z,\\\
\eta_{2i}=\frac{1}{2}(\zeta_{i}+\zeta_{i+1}),\quad&\text{for
}i=1,\dots,z-1.\end{cases}$

With this notation we have $F(\eta_{1})=0$ and for $i=2,\dots,2z-1$:

$\displaystyle F(\eta_{i})$
$\displaystyle=\int_{0}^{\eta_{i}}f(s)\mathrm{d}s=\sum_{j=1}^{i-1}\int_{\eta_{j}}^{\eta_{j+1}}f(s)\mathrm{d}s$
$\displaystyle\approx\sum_{j=1}^{i-1}\frac{h}{6}\left(f(\eta_{j})+4f\Big{(}\frac{\eta_{j}+\eta_{j+1}}{2}\Big{)}+f(\eta_{j+1})\right),$

where $h=\eta_{j+1}-\eta_{j}$. Notice that this quadrature formula is exact
for polynomials up to degree $3$, but it requires evaluation of the
integrating function over further midpoints, i.e. $(\eta_{j}+\eta_{j+1})/2$,
as it is shown in Figure 5.

$0$$0.2$$0.4$$0.6$$0.8$$1$$-1$$-0.5$$0$$0.5$$1$$f(\eta_{j})$$f(\eta_{j+1})$$f(\eta_{j+2})$$f\left(\dfrac{\eta_{j}+\eta_{j+1}}{2}\right)$$f\left(\dfrac{\eta_{j+1}+\eta_{j+2}}{2}\right)$$x$$f(x)=sin(2\pi
x)$$f(x)$ on breakpoints$f(x)$ on midpointson further midpoints Figure 5:
Example of evaluation of $f(x)=sin(2\pi x)$ in further midpoints, in order to
project with $\llm[p]^{c}$

####  4.2.2 Multivariate quasi-interpolant projections

Finally we define the multivariate quasi-interpolant projections that commute
with the divergence operator in the parametric domain, as

(28) $\displaystyle\LLM[1]$
$\displaystyle=(\llm[p_{1}]\otimes\llm[p_{2}-1]^{c})\times(\llm[p_{1}-1]^{c}\otimes\llm[p_{2}]),$
$\displaystyle\LLM[2]$
$\displaystyle=\llm[p_{1}-1]^{c}\otimes\llm[p_{2}-1]^{c}.$

The projections in the physical domain, namely $\pLLM[1],\pLLM[2]$, are
uniquely determined by (18). We remark that projections
$\LLM[1],\LLM[2],\pLLM[1],\pLLM[2]$ satisfy 1. An analogous constructions
holds for projections in periodic spline spaces, or in spline spaces with
mixed Dirichlet and periodic conditions, as in 1.

###  4.3 Conservative time discretization

Regarding time discretization, the choice of an energy-preserving method is
mandatory since we are interested in preserving the total energy of the
system. We detail the discretization with Crank-Nicolson method, which is
second order and energy conservative, but other conservative methods could be
chosen.

Let us fix $\tau=(t_{0},\dots,t_{N_{T}})$, a partition of the interval
$[0,T]$, such that $t_{0}=0$, $t_{N_{T}}=T$ and $t_{n}<t_{n+1}$. For a
simplified notation we assume the partition to be uniform, that is
$t_{n+1}-t_{n}=k$, with a fixed real number $k>0$, for $n=0,\dots,N_{T}-1$.
Let us denote by $\mathbf{v}_{h}^{n}$ and $\phi_{h}^{n}$ respectively the
value of $\mathbf{v}_{h}$ and $\phi_{h}$ at the time instant $t_{n}$. Applying
Crank-Nicolson method to equations (25), the fully discrete problem is now,
for $n=0,\dots,N_{T}-1$ find $\mathbf{v}_{h}^{n+1}\in X^{1}_{h}$ and
$\phi_{h}^{n+1}\in X^{2}_{h}$, such that

(29a)
$\displaystyle\int_{\Omega}\frac{\mathbf{v}_{h}^{n+1}-\mathbf{v}_{h}^{n}}{k}\cdot\mathbf{w}_{h}\mathrm{d}\mathbf{x}$
$\displaystyle=-\frac{1}{2}\int_{\Omega}\Pi^{2}(\text{div}(c\mathbf{w}_{h}))(\phi_{h}^{n+1}+\phi_{h}^{n})\mathrm{d}\mathbf{x},\quad\forall\mathbf{w}_{h}\in
X^{1}_{h},$ (29b) $\displaystyle\frac{\phi_{h}^{n+1}-\phi_{h}^{n}}{k}$
$\displaystyle=+\frac{1}{2}\big{(}\Pi^{2}(\text{div}(c\mathbf{v}_{h}^{n+1}))+\Pi^{2}(\text{div}(c\mathbf{v}_{h}^{n}))\big{)}.$

Notice that in equation (29b), the left hand side $\phi_{h}^{n+1}$ is written
as a linear combination of terms that depend on the time instant $t_{n}$,
except for $\mathbf{v}_{h}^{n+1}$. Moreover, in view of the definition of the
operator $\Pi^{2}$, (29b) can be seen as a collocation of the equation (6b).
To solve the system, we replace the expression of $\phi_{h}^{n+1}$ from the
second equation into the first one, and bring to the left all the terms in
which the unknown $\mathbf{v}_{h}^{n+1}$ appears, that is

(30)
$\int_{\Omega}\mathbf{v}_{h}^{n+1}\cdot\mathbf{w}_{h}\mathrm{d}\mathbf{x}+\frac{k^{2}}{4}\int_{\Omega}\Pi^{2}(\text{div}(c\mathbf{v}_{h}^{n+1}))\Pi^{2}(\text{div}(c\mathbf{w}_{h}))\mathrm{d}\mathbf{x}=\\\
=\int_{\Omega}\mathbf{v}_{h}^{n}\cdot\mathbf{w}_{h}\mathrm{d}\mathbf{x}-\frac{k^{2}}{4}\int_{\Omega}\Pi^{2}(\text{div}(c\mathbf{v}_{h}^{n}))\Pi^{2}(\text{div}(c\mathbf{w}_{h}))\mathrm{d}\mathbf{x}-k\int_{\Omega}\phi_{h}^{n}\Pi^{2}(\text{div}(c\mathbf{w}_{h}))\mathrm{d}\mathbf{x},$

which must hold for all $\mathbf{w}_{h}\in X^{1}_{h}$. We solve this equation
to compute the unknown $\mathbf{v}_{h}^{n+1}$, which is then used to update
the solution $\phi_{h}^{n+1}$ as indicated by the second equation of (29).

In our numerical tests, we will to compare the proposed method with the
standard Galerkin formulation, that we present for completeness, and which is
given by: find $\mathbf{v}_{h}\in X^{1}_{h}$ and $\phi_{h}\in X^{2}_{h}$, such
that

(31a)
$\displaystyle\int_{\Omega}\frac{\mathbf{v}_{h}^{n+1}-\mathbf{v}_{h}^{n}}{k}\cdot\mathbf{w}_{h}\mathrm{d}\mathbf{x}$
$\displaystyle=-\frac{1}{2}\int_{\Omega}\text{div}(c\mathbf{w}_{h})(\phi_{h}^{n+1}+\phi_{h}^{n})\mathrm{d}\mathbf{x},\quad\forall\mathbf{w}_{h}\in
X^{1}_{h},$ (31b)
$\displaystyle\int_{\Omega}\frac{\phi_{h}^{n+1}-\phi_{h}^{n}}{k}\psi_{h}\mathrm{d}\mathbf{x}$
$\displaystyle=+\frac{1}{2}\int_{\Omega}\big{(}\text{div}(c\mathbf{v}_{h}^{n+1})+\text{div}(c\mathbf{v}_{h}^{n})\big{)}\psi_{h}\mathrm{d}\mathbf{x},\quad\forall\psi_{h}\in
X^{2}_{h},$

######  Remark 2.

Notice that both equations (29) and (31) are energy conservative schemes by
construction. On the other hand, (31) can not be reduced to a single equation.
Thus, the method (29) is intrinsically cheaper than the method (31) and, at
our knowledge, it is the first conservative method that does not require
computation of both primal and dual variables at each time step. REWRITE THIS
SENTENCE? DISCUSS

##  5 Error convergence analysis

In this section we analyze the convergence of the proposed method. In Section
5.1 we first present a convergence study in an abstract setting, for a generic
family of projections as in Section 4.1.4, assuming stability and
approximation properties both for the projectors and for their adjoint
operators. Unfortunately, the adjoints of the quasi-interpolants of Section
4.2 do not satisfy these assumptions. For this reason, we present in Section
5.2 the analysis under some relaxed assumptions, for which only linear
convergence could be proved.

It will be useful for the analysis to write the equations of the problem in
weak form. Let us denote by $(\mathbf{v},\phi)$ the classical solution of the
problem Eq. 2, while $(\mathbf{v}_{h},\phi_{h})$ is the solution of the
semidiscrete problem in space (25), which satisfies also the following
variational equations:

(32a)
$\displaystyle\int_{\Omega}\left(\mathbf{v}_{h}\right)_{t}\cdot\mathbf{w}_{h}\
\mathrm{d}\mathbf{x}+\int_{\Omega}\Pi^{2}(\text{div}(c\mathbf{w}_{h}))\
\phi_{h}\ \mathrm{d}\mathbf{x}$ $\displaystyle=0,\quad$
$\displaystyle\forall\mathbf{w}_{h}\in{X}^{1}_{h},$ (32b)
$\displaystyle\int_{\Omega}\left(\phi_{h}\right)_{t}\psi_{h}\
\mathrm{d}\mathbf{x}-\int_{\Omega}\Pi^{2}(\text{div}(c\mathbf{v}_{h}))\
\psi_{h}\ \mathrm{d}\mathbf{x}$ $\displaystyle=0,\quad$
$\displaystyle\forall\psi_{h}\in{X}^{2}_{h}.$

We introduce the energy norm
$\|\mathbf{v},\phi\|_{E}^{2}\coloneqq\|\mathbf{v}\|_{\mathbf{L}^{2}(\Omega)}^{2}+\|\phi\|_{L^{2}(\Omega)}^{2}$,
the Sobolev norm
$\|\mathbf{v},\phi\|_{\mathcal{H}^{s}}^{2}\coloneqq\|\mathbf{v}\|_{\mathbf{H}^{s}(\text{div};\Omega)}^{2}+\|\phi\|_{H^{s}(\Omega)}^{2}$,
for a positive integer $s\geq 0$, and the following norms in the space time
domain $\Omega\times[0,T]$

$\displaystyle\|\mathbf{v},\phi\|_{\infty,E}$
$\displaystyle\coloneqq\sup_{t\in[0,T]}\|\mathbf{v},\phi\|_{E},$
$\displaystyle\|\mathbf{v},\phi\|_{W^{1,\infty},\mathcal{H}^{s}}$
$\displaystyle\coloneqq\max\\{\sup_{t\in[0,T]}\|\mathbf{v},\phi\|_{\mathcal{H}^{s}},\sup_{t\in[0,T]}\|\mathbf{v}_{t},\phi_{t}\|_{\mathcal{H}^{s}}\\}.$

In order to alleviate notation, we indicate by
$\|\phi,\psi\|_{E}^{2}=\|\phi\|_{L^{2}(\Omega)}^{2}+\|\psi\|_{L^{2}(\Omega)}^{2}$
the equivalent of the energy norm when $\phi,\psi\in L^{2}({\Omega})$ are both
scalar fields. We are interested in bounding

(33)
$\|\mathbf{v}-\mathbf{v}_{h},\phi-\phi_{h}\|_{\infty,E}=\sup_{t\in[0,T]}\|\mathbf{v}-\mathbf{v}_{h},\phi-\phi_{h}\|_{E},$

and we will follow the ideas used in [2] for elasticity problems in mixed form
with weak symmetry.

###  5.1 Convergence analysis based on the adjoint projections

Let us introduce the adjoint projection $\Pi^{2,*}\in\textit{End}(X^{2})$,
which is defined by the following fundamental relation

(34)
$\prescript{}{X^{2}}{\langle}\Pi^{2}(\phi),\psi\rangle_{X^{2,*}}=\prescript{}{X^{2}}{\langle}\phi,\Pi^{2,*}(\psi)\rangle_{X^{2,*}},\quad\forall\
\phi\in X^{2},\ \psi\in X^{2,*},$

where $X^{2,*}$ is the dual space of $X^{2}$, and $\langle\cdot,\cdot\rangle$
denotes the duality pairing. We identify $X^{2}=L^{2}(\Omega)$ with its dual
space, as it is usually done, and in particular (34) can be expressed in terms
of the $L^{2}(\Omega)$ scalar product. Throughout this section, we require
that $\Pi^{2}$ and $\Pi^{2,*}$ are $L^{2}$-stable, and $\Pi^{1},\Pi^{2}$ and
$\Pi^{2*}$ have good approximation properties, summarizing our request in the
following assumption.

######  Assumption 3.

The projections $\Pi^{2}$ and $\Pi^{2,*}$ are $L^{2}$-stable, that is

$\displaystyle\|\Pi^{2}(\phi)\|_{L^{2}(\Omega)}\leq\|\phi\|_{L^{2}(\Omega)},\quad\text{and}\quad\|\Pi^{2,*}(\psi)\|_{L^{2}(\Omega)}\leq\|\psi\|_{L^{2}(\Omega)},\quad\forall\phi,\,\psi\in
L^{2}(\Omega).$

Moreover, given $\mathbf{w}\in\mathbf{H}^{m}(\mathrm{div};\Omega)$ and
$\psi\in H^{m}(\Omega)$, for $0\leq l\leq m\leq p$, and
$p=\min_{i=1}^{d}\\{p_{i}\\}$, there exists a constant $C$ such that it holds:

(35a)
$\displaystyle\|(Id-\Pi^{1})\mathbf{w}\|_{\mathbf{H}^{l}(\mathrm{div};\Omega)}$
$\displaystyle\leq
Ch^{m-l}\|\mathbf{w}\|_{\mathbf{H}^{m}(\mathrm{div};\Omega)},$ (35b)
$\displaystyle\|(Id-\Pi^{2})\psi\|_{{H}^{l}(\Omega)}$ $\displaystyle\leq
Ch^{m-l}\|\psi\|_{{H}^{m}(\Omega)},$ (35c)
$\displaystyle\|(Id-\Pi^{2,*})\psi\|_{{H}^{l}(\Omega)}$ $\displaystyle\leq
Ch^{m-l}\|\psi\|_{{H}^{m}(\Omega)},$

where $Id$ is the identity operator, and $h$ is the mesh size.

Projections that commute with the divergence operator and satisfy
$L^{2}$-stability, together with (35a) and (35b), are well known in
literature, see [10, Remark 5.1]. The usefulness of this working hypothesis is
thus to ensure (35c). One option is to explore quasi-interpolant operators
constructed with biorthogonal dual bases, the latter having the capability of
reproducing polynomials, see i.e. [30].

Let us assume from now on that the classical solution is sufficiently regular,
i.e. $\mathbf{v}\in\mathbf{H}^{m}(\text{div};\Omega)$ and $\phi\in
H^{m}(\Omega)$, for $0\leq m\leq p$ and $p$ as in Assumption 3. We can start
bounding (33) by using the triangular inequality, and separate the error in
the following way:

(36)
$\|\mathbf{v}-\mathbf{v}_{h},\phi-\phi_{h}\|_{\infty,E}\leq\|\mathbf{v}-\mathbf{v}_{h}^{P},\phi-\phi_{h}^{P}\|_{\infty,E}+\|\mathbf{v}_{h}^{P}-\mathbf{v}_{h},\phi_{h}^{P}-\phi_{h}\|_{\infty,E},$

where $\mathbf{v}_{h}^{P}$ and $\phi_{h}^{P}$ are suitable approximations of
the fields $\mathbf{v}$ and $\phi$. Here we fix
$\mathbf{v}_{h}^{P}\coloneqq\Pi^{1}(\mathbf{v})$ and
$\phi_{h}^{P}\coloneqq\Pi^{2}(\phi)$, but other approximation choices that
satisfy (35a) and (35b) can also be considered. We recall the following
approximation result.

######  Lemma 2.

There exists a constant $C>0$ such that, for $0\leq m\leq p$ and $p$ as in
Assumption 3, it holds

(37a)
$\displaystyle\|{\mathbf{v}-\mathbf{v}_{h}^{P},\phi-\phi_{h}^{P}}\|_{\infty,E}$
$\displaystyle\leq Ch^{m}\|{\mathbf{v},\phi}\|_{\infty,\mathcal{H}^{m}},$
(37b) $\displaystyle\|{\phi-\phi_{h}^{P},\
\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})}\|_{\infty,E}$
$\displaystyle\leq Ch^{m}\|{\mathbf{v},\phi}\|_{\infty,\mathcal{H}^{m}},$
(37c)
$\displaystyle\|{(\mathbf{v}-\mathbf{v}_{h}^{P})_{t},(\phi-\phi_{h}^{P})_{t}}\|_{\infty,E}$
$\displaystyle\leq
Ch^{m}\|{\mathbf{v}_{t},\phi_{t}}\|_{\infty,\mathcal{H}^{m}}.$

###### Proof.

In order to bound (37a) we argue like this:

$\displaystyle\|{\mathbf{v}-\mathbf{v}_{h}^{P},\phi-\phi_{h}^{P}}\|_{\infty,E}^{2}$
$\displaystyle\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny
def}}}}{{=}}}\sup_{t\in[0,T]}\left(\|\mathbf{v}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}^{2}+\|\phi-\phi_{h}^{P}\|_{L^{2}(\Omega)}^{2}\right)$
$\displaystyle\leq\sup_{t\in[0,T]}\left(Ch^{2m}\|\mathbf{v}\|_{\mathbf{H}^{m}(\text{div};\Omega)}^{2}+Ch^{2m}\|\phi\|_{H^{m}(\Omega)}^{2}\right)$
$\displaystyle=Ch^{2m}\|{\mathbf{v},\phi}\|_{\infty,\mathcal{H}^{m}}^{2},$

where we used only (35a) and (35b), and taking the square roots we proved
(37a). The proof of (37c) is analogous. As regards (37b) we have that

$\|{\phi-\phi_{h}^{P},\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})}\|_{\infty,E}^{2}=\sup_{t\in[0,T]}\left(\|\phi-\phi_{h}^{P}\|_{L^{2}(\Omega)}^{2}+\|\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}^{2}\right).$

By applying the Leibnitz rule and Young’s inequality, $(a+b)^{2}\leq
2a^{2}+2b^{2}$ for each $a,b>0$, on the last term of the right hand side we
have

$\displaystyle\|\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\left(\|c\
\textup{div}(\mathbf{v}-\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}+\|\gradient{c}\cdot(\mathbf{v}-\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}\right)^{2}$
$\displaystyle\leq 2\|c\
\textup{div}(\mathbf{v}-\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}^{2}+2\|\gradient{c}\cdot(\mathbf{v}-\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq
C\|\mathbf{v}-\mathbf{v}_{h}^{P}\|_{\mathbf{H}^{0}(\text{div};\Omega)}^{2},$

because $c$ is smooth and uniformly bounded. Now again we can apply (35a) and
(35b), which gives us

$\displaystyle\|{\phi-\phi_{h}^{P},\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})}\|_{\infty,E}^{2}$
$\displaystyle\leq
C\sup_{t\in[0,T]}\left(\|\phi-\phi_{h}^{P}\|_{L^{2}(\Omega)}^{2}+\|\mathbf{v}-\mathbf{v}_{h}^{P}\|_{\mathbf{H}^{0}(\text{div};\Omega)}^{2}\right)$
$\displaystyle\leq
C\sup_{t\in[0,T]}\left(h^{2m}\|\phi\|_{H^{m}(\Omega)}^{2}+h^{2m}\|\mathbf{v}\|_{\mathbf{H}^{m}(\text{div};\Omega)}^{2}\right)$
$\displaystyle\leq Ch^{2m}\|{\mathbf{v},\phi}\|_{\infty,\mathcal{H}^{m}}^{2}.$

and again taking the square roots we end the proof of this lemma.

Next we state the main result of this paper.

######  Theorem 1.

Under Assumption 3, together with the commutativity of the projectors,
$\Pi^{2}\textup{div}=\textup{div}\Pi^{1}$, given
$\mathbf{v}\in\mathbf{H}^{m}(\mathrm{div};\Omega)$ and $\phi\in
H^{m}(\Omega)$, it holds that

(38) $\|\mathbf{v}-\mathbf{v}_{h},\phi-\phi_{h}\|_{\infty,E}\leq
Ch^{m-1}\|\mathbf{v},\phi\|_{W^{1,\infty},\mathcal{H}^{m}}.$

###### Proof.

Let us consider both the weak formulation (3) and the semi-discrete problem
(32). With the notation $(\cdot,\cdot)$ indicating the $L^{2}$ scalar product
over the physical domain $\Omega$ we can write:

$\displaystyle\left((\mathbf{v}_{h})_{t},\mathbf{w}_{h}\right)+\left(\Pi^{2}\text{div}(c\mathbf{w}_{h}),\phi_{h}\right)$
$\displaystyle=\left(\mathbf{v}_{t},\mathbf{w}_{h}\right)+\left(\text{div}(c\mathbf{w}_{h}),\phi\right),$
$\displaystyle\left((\phi_{h})_{t},\psi_{h}\right)-\left(\Pi^{2}\text{div}(c\mathbf{v}_{h}),\psi_{h}\right)$
$\displaystyle=\left(\phi_{t},\psi_{h}\right)-\left(\text{div}(c\mathbf{v}),\psi_{h}\right),$

for all $\mathbf{w}_{h}\in X^{1}_{h}$ and $\psi_{h}\in X^{2}_{h}$. We can
respectively subtract from the previous equations the following quantities:

$\left((\mathbf{v}_{h}^{P})_{t},\mathbf{w}_{h}\right)+\left(\Pi^{2}\textup{div}(c\mathbf{w}_{h}),\phi_{h}^{P}\right),\quad\text{and}\quad\left((\phi_{h}^{P})_{t},\psi_{h}\right)-\left(\Pi^{2}\textup{div}(c\mathbf{v}_{h}^{P}),\psi_{h}\right),$

and rearrange the two equations such that it holds that

(39)
$\left((\mathbf{v}_{h}-\mathbf{v}_{h}^{P})_{t},\mathbf{w}_{h}\right)+\left(\Pi^{2}\textup{div}(c\mathbf{w}_{h}),(\phi_{h}-\phi_{h}^{P})\right)=\\\
=\left((\mathbf{v}-\mathbf{v}_{h}^{P})_{t},\mathbf{w}_{h}\right)+\left(\Pi^{2}\textup{div}(c\mathbf{w}_{h}),(\phi-\phi_{h}^{P})\right)+\left((Id-\Pi^{2})\textup{div}(c\mathbf{w}_{h}),\phi\right),$

and

(40)
$\left((\phi_{h}-\phi_{h}^{P})_{t},\psi_{h}\right)-\left(\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P}),\psi_{h}\right)=\\\
=\left((\phi-\phi_{h}^{P})_{t},\psi_{h}\right)-\left(\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P}),\psi_{h}\right)-\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}),\psi_{h}\right).$

By making the choice $\mathbf{w}_{h}=\mathbf{v}_{h}-\mathbf{v}_{h}^{P}$ and
$\psi_{h}=\phi_{h}-\phi_{h}^{P}$, and adding together the equations we get

(41)
$\displaystyle\begin{split}\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}&-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}^{2}=\\\
&=\left((\mathbf{v}-\mathbf{v}_{h}^{P})_{t},\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\right)+\left((\phi-\phi_{h}^{P})_{t},\phi_{h}-\phi_{h}^{P}\right)+\\\
&+\left(\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P}),\phi-\phi_{h}^{P}\right)-\left(\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P}),\phi_{h}-\phi_{h}^{P}\right)+\\\
&+\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P}),\phi\right)-\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}),\phi_{h}-\phi_{h}^{P}\right).\end{split}$

Here, by construction $\Pi^{2}\textup{div}=\textup{div}\,\Pi^{1}$. Combining
this with the definition of the formal adjoint operator, we can rewrite the
last two terms on the right hand side as

(42)
$\displaystyle\begin{split}\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P}),\phi\right)&=\left(\textup{div}(c\mathbf{v}_{h}^{P}-c\mathbf{v}_{h}),(Id-\Pi^{2,*})\phi\right),\\\
\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}),\phi_{h}-\phi_{h}^{P}\right)&=\left(\textup{div}(Id-\Pi^{1})(c\mathbf{v}),\phi_{h}-\phi_{h}^{P}\right).\end{split}$

Now on the right hand side we can apply the Cauchy-Schwartz inequality, which
gives

$\displaystyle\begin{split}\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}^{2}\leq\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t}\|_{\mathbf{L}^{2}(\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+\|(\phi-\phi_{h}^{P})_{t}\|_{L^{2}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+\|(\phi-\phi_{h}^{P})\|_{L^{2}(\Omega)}&\|\Pi^{2}\textup{div}\left(c\mathbf{v}_{h}^{P}-c\mathbf{v}_{h}\right)\|_{L^{2}(\Omega)}+\\\
+\|\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+\|(Id-\Pi^{2,*})\phi\|_{L^{2}(\Omega)}&\|\textup{div}(c\mathbf{v}_{h}^{P}-c\mathbf{v}_{h})\|_{L^{2}(\Omega)}+\\\
+\|\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{L^{2}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}.\end{split}$

By boundedness of the projections as in Assumption 3 and using boundedness of
the coefficient $c$, together with inverse inequalities, there exists a
constant $C>0$ such that:

$\displaystyle\begin{split}\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}^{2}\leq\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t}\|_{\mathbf{L}^{2}(\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+\|(\phi-\phi_{h}^{P})_{t}\|_{L^{2}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+Ch^{-1}\|(\phi-\phi_{h}^{P})\|_{L^{2}(\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+C\|\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+Ch^{-1}\|(Id-\Pi^{2,*})\phi\|_{L^{2}(\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+\|\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{L^{2}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}.\end{split}$

Again by Cauchy-Schwartz inequality, applied to ordered couples of terms in
the right hand side, we have:

$\displaystyle\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{E}^{2}\leq\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t},\
(\phi-\phi_{h}^{P})_{t}\|_{E}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{E}$ $\displaystyle+$
$\displaystyle+Ch^{-1}\|(\phi-\phi_{h}^{P}),\
\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{E}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{E}$ $\displaystyle+$
$\displaystyle+Ch^{-1}\|(Id-\Pi^{2,*})\phi,\
\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{E}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{E}$ .

Dividing both sides by $\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{E}$ and integrating in time in $[0,t]$, we obtain

$\displaystyle\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{E}$
$\displaystyle\leq\|\mathbf{v}_{h}(\mathbf{x},0)-\mathbf{v}_{h}^{P}(\mathbf{x},0),\
\phi_{h}(\mathbf{x},0)-\phi_{h}^{P}(\mathbf{x},0)\|_{E}+$
$\displaystyle+\int_{0}^{t}\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t},\
(\phi-\phi_{h}^{P})_{t}\|_{E}\mathrm{d}\tau+$
$\displaystyle+\int_{0}^{t}Ch^{-1}\|(\phi-\phi_{h}^{P}),\
\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{E}\mathrm{d}\tau+$
$\displaystyle+\int_{0}^{t}Ch^{-1}\|(Id-\Pi^{2,*})\phi,\
\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{E}\mathrm{d}\tau.$

Specifying the initial condition for problem (32) to be
$\left(\mathbf{v}_{h}^{P}(\mathbf{x},0),\phi_{h}^{P}(\mathbf{x},0)\right)$ it
holds $\|\mathbf{v}_{h}(\mathbf{x},0)-\mathbf{v}_{h}^{P}(\mathbf{x},0),\
\phi_{h}(\mathbf{x},0)-\phi_{h}^{P}(\mathbf{x},0)\|_{E}=0$. Since the previous
estimate is valid for almost every $t\in[0,T]$ we can take the maximum, and
there exists a new constant $C>0$ such that

$\displaystyle\begin{split}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\
\phi_{h}-\phi_{h}^{P}\|_{\infty,E}\leq
C&\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t},\
(\phi-\phi_{h}^{P})_{t}\|_{\infty,E}+\\\ +Ch^{-1}&\|(\phi-\phi_{h}^{P}),\
\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{\infty,E}+\\\
+Ch^{-1}&\|(Id-\Pi^{2,*})\phi,\
\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{\infty,E}.\end{split}$

Now the left hand side in (38) splits as in (36), therefore we have

$\displaystyle\begin{split}\|\mathbf{v}-\mathbf{v}_{h},\
\phi-\phi_{h}\|_{\infty,E}\leq\ &\|\mathbf{v}-\mathbf{v}_{h}^{P},\
\phi-\phi_{h}^{P}\|_{\infty,E}+\\\ +C&\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t},\
(\phi-\phi_{h}^{P})_{t}\|_{\infty,E}+\\\ +Ch^{-1}&\|(\phi-\phi_{h}^{P}),\
\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{\infty,E}+\\\
+Ch^{-1}&\|(Id-\Pi^{2,*})\phi,\
\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{\infty,E}.\end{split}$

Finally, from Lemma 2 and Assumption 3 we obtain

$\displaystyle\begin{split}\|\mathbf{v}-\mathbf{v}_{h},\
\phi-\phi_{h}\|_{\infty,E}\leq C&\Big{(}h^{m-1}\|\mathbf{v},\
\phi\|_{\infty,\mathcal{H}^{m}}+h^{m}\|\mathbf{v}_{t},\
\phi_{t}\|_{\infty,\mathcal{H}^{m}}\Big{)},\end{split}$

which concludes the proof of this theorem.

######  Remark 3.

The suboptimality of the estimate comes from using inverse inequalities to
control the divergence of the error. Optimality can be recovered if we define
suitable approximations $\mathbf{v}_{h}^{P}$ and $\phi_{h}^{P}$ that control
the divergence, such as in [8]. However, the projections $\pLLM[1]$ and
$\pLLM[2]$, that we implement and test numerically, do not satisfy Assumption
3, therefore we directly move to next section, dedicated to error estimates
for these projections.

###  5.2 Weaker approximation assumptions

In this section we proof a convergence result for weaker assumptions than the
ones in Assumption 3. The reason is the following. Of considerable practical
interest are projections in discrete spaces that have good approximation
properties and are computationally fast and efficient to compute, i.e. [23,
21]. A very general family of such projections are the quasi-interpolators
proposed in [25]. Since we do not have good approximation properties of the
adjoint projection of such quasi-interpolants, we will require here to relax
the assumptions we made in the previous section.

######  Assumption 4.

The projections $\Pi^{1}$ and $\Pi^{2}$ satisfy

$\|\Pi^{1}\|_{\mathbf{L}^{\infty}(\widetilde{K})\to\mathbf{L}^{\infty}(K)}=C_{1},\quad\text{and}\quad\|\Pi^{2}\|_{{L}^{\infty}(\widetilde{K})\to{L}^{\infty}(K)}=C_{2}.$

where $K\in\mathcal{M}$ and $\widetilde{K}$ is its support extension.
Moreover, given $\mathbf{w}\in\mathbf{H}^{m}(\mathrm{div};\Omega)$ and
$\psi\in H^{m}(\Omega)$, with at least $m\geq 2$, for $0\leq l<m\leq p$, and
$p=\min_{i=1}^{d}\\{p_{i}\\}$, there exists a constant $C$ such that it holds:

(43a)
$\displaystyle\|(Id-\Pi^{1})\mathbf{w}\|_{\mathbf{H}^{l}(\mathrm{div};\Omega)}$
$\displaystyle\leq
Ch^{m-l}\|\mathbf{w}\|_{\mathbf{H}^{m}(\mathrm{div};\Omega)},$ (43b)
$\displaystyle\|(Id-\Pi^{2})\psi\|_{{H}^{l}(\Omega)}$ $\displaystyle\leq
Ch^{m-l}\|\psi\|_{{H}^{m}(\Omega)},$

where $Id$ is the identity operator, and $h$ is the mesh size.

Notice that, the projectors $\pLLM[1]$ and $\pLLM[2]$ satisfy Assumption 4,
but not Assumption 3. We consider again
$\mathbf{v}_{h}^{P}=\Pi^{1}(\mathbf{v})$ and $\phi_{h}^{P}=\Pi^{2}(\phi)$, and
we notice that under Assumption 4 it still holds 2. Therefore we have the
following theorem.

######  Theorem 2.

Under Assumption 4, together with the commutativity of the projectors,
$\Pi^{2}\textup{div}=\textup{div}\Pi^{1}$, given
$\mathbf{v}\in\mathbf{H}^{m}(\mathrm{div};\Omega)$ and $\phi\in
H^{m}(\Omega)$, with at least $m\geq 2$, and
$c\in\mathcal{C}^{\infty}(\Omega)$, it holds that

(44) $\|\mathbf{v}-\mathbf{v}_{h},\phi-\phi_{h}\|_{\infty,E}\leq
Ch\|\mathbf{v},\phi\|_{W^{1,\infty},\mathcal{H}^{m}}.$

###### Proof.

The first steps of the proof are as in Theorem 1 until the choice of the test
functions, and in particular (41) is valid. Then, by construction we have
$\Pi^{2}\textup{div}=\textup{div}\Pi^{1}$, and by integrating by parts we
obtain

$\displaystyle\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P}),\phi\right)=$
$\displaystyle-\left((Id-\Pi^{1})(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P}),\gradient{\phi}\right)$
$\displaystyle\left((Id-\Pi^{2})\textup{div}(c\mathbf{v}),\phi_{h}-\phi_{h}^{P}\right)=$
$\displaystyle\left(\textup{div}(Id-\Pi^{1})(c\mathbf{v}),\phi_{h}-\phi_{h}^{P}\right),$

note the difference in the first equation with respect to (42). By using
Cauchy-Schwartz inequality as in the proof of Theorem 1, we have

$\displaystyle\begin{split}\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}^{2}\leq\|(\mathbf{v}-\mathbf{v}_{h}^{P})_{t}\|_{\mathbf{L}^{2}(\Omega)}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+\|(\phi-\phi_{h}^{P})_{t}\|_{L^{2}(\Omega)}\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+\|\phi-\phi_{h}^{P}\|_{L^{2}(\Omega)}{\|\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(\Omega)}}+\\\
+{\|\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}(\Omega)}}\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+\|\phi\|_{H^{1}(\Omega)}\|(Id-\Pi^{1})(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(\Omega)}+\\\
+\|\textup{div}(Id-\Pi^{1})(c\mathbf{v})\|_{L^{2}(\Omega)}\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)},\end{split}$

and the only differences with respect to the proof of the previous theorem are
in the stability assumptions of $\Pi^{2}$ and in the fifth term of the right-
hand side. In order to bound this last one, we can apply the superconvergence
results stated in [6, Theorem 2.2], that is

(45)
$\|(Id-\Pi^{1})(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(\Omega)}\leq
Ch\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}.$

Next, we bound
$\|\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(\Omega)}$
by local arguments. Given $K\in\mathcal{M}$, since $|K|\leq h^{d}$ where $d$
is the dimension of $\Omega$, we first apply Hölder inequality, then stability
of Assumption 4 together with uniform boundedness of $c$ and finally inverse
inequalities for splines, that is

$\displaystyle\begin{split}\|\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(K)}^{2}&\leq|K|\,\|\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{\infty}(K)}^{2}\\\
&\leq h^{d}C_{2}^{2}\|\nabla
c\cdot(\mathbf{v}_{h}-\mathbf{v}_{h}^{P})+c\,\textup{div}(\mathbf{v}_{h}-\mathbf{v}_{h}^{P})\|_{{L}^{\infty}(\widetilde{K})}^{2}\\\
&\leq
h^{d}C\left(\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{\infty}(\widetilde{K})}+\|\textup{div}(\mathbf{v}_{h}-\mathbf{v}_{h}^{P})\|_{{L}^{\infty}(\widetilde{K})}\right)^{2}\\\
&\leq
h^{d}C\left(\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{\infty}(\widetilde{K})}+h^{-1}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{\infty}(\widetilde{K})}\right)^{2}\\\
&\leq
h^{d}Ch^{-2}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{\infty}(\widetilde{K})}^{2}\\\
&\leq
Ch^{-2}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\widetilde{K})}^{2}.\end{split}$

Notice that in the last inequality it is essential to have a shape regularity
assumption of the mesh. Taking the square root and by standard arguments we
have immediately the global inequality

$\|\Pi^{2}\textup{div}(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(\Omega)}\leq
C{h^{-1}}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}.$

Recall that $\mathbf{v}_{h}^{P}=\Pi^{1}(\mathbf{v})$ and
$\phi_{h}^{P}=\Pi^{2}(\phi)$, together with
$\textup{div}\Pi^{1}=\Pi^{2}\textup{div}$. Again, with the same arguments, and
with $|K|\leq h^{d}$, we have

$\displaystyle\begin{split}\|\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}({K})}^{2}&\leq
h^{d}\|\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{\infty}(K)}^{2}\\\
&\leq h^{d}C_{2}^{2}\|\nabla
c\cdot(\mathbf{v}-\mathbf{v}_{h}^{P})+c\,\textup{div}(\mathbf{v}-\mathbf{v}_{h}^{P})\|_{L^{\infty}(\widetilde{K})}^{2}\\\
&\leq
h^{d}C\left(\|\mathbf{v}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{\infty}(\widetilde{K})}+\|\textup{div}(\mathbf{v}-\mathbf{v}_{h}^{P})\|_{{L}^{\infty}(\widetilde{K})}\right)^{2}\\\
&\leq
h^{d}C\left(\|(Id-\Pi^{1})\mathbf{v}\|_{\mathbf{L}^{\infty}(\widetilde{K})}^{2}+\|(Id-\Pi^{2})\textup{div}(\mathbf{v})\|_{{L}^{\infty}(\widetilde{K})}^{2}\right)\\\
&\leq
h^{d}C\left(h^{2m-d}\|\mathbf{v}\|_{\mathbf{H}^{m}(\widetilde{\widetilde{K}})}^{2}+h^{2m-d}\|\textup{div}(\mathbf{v})\|_{{H}^{m}(\widetilde{\widetilde{K}})}^{2}\right)\\\
&\leq
Ch^{2m}\|\mathbf{v}\|_{H^{m}(\textup{div};\widetilde{\widetilde{K}})}^{2},\end{split}$

where $\widetilde{\widetilde{K}}$ is the support extension of $\widetilde{K}$.
Notice that we used Cauchy-Schwartz inequality on fourth row, and the local
approximation estimates of the kind [27, Theorem 10.2] on fifth row. Taking
the square root and by standard arguments this global inequality is
straightforward

$\|\Pi^{2}\textup{div}(c\mathbf{v}-c\mathbf{v}_{h}^{P})\|_{L^{2}({\Omega})}\leq
Ch^{m}\|\mathbf{v}\|_{H^{m}(\textup{div};\Omega)}.$

We bound the remaining of approximation errors as it is done for Lemma 2, and
together with the above inequalities, we have

$\displaystyle\begin{split}\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}^{2}\leq
Ch^{m}\|\mathbf{v}_{t}\|_{\mathbf{H}^{m}(\textup{div};\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+Ch^{m}\|\phi_{t}\|_{H^{m}(\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+Ch^{m-1}\|\phi\|_{H^{m}(\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+Ch^{m}\|\mathbf{v}\|_{\mathbf{H}^{m}(\textup{div};\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}+\\\
+Ch\|\phi\|_{H^{1}(\Omega)}&\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P}\|_{\mathbf{L}^{2}(\Omega)}+\\\
+Ch^{m}\|\mathbf{v}\|_{\mathbf{H}^{m}(\textup{div};\Omega)}&\|\phi_{h}-\phi_{h}^{P}\|_{L^{2}(\Omega)}.\end{split}$

By using Cauchy-Schwartz inequality and adding together the similar terms, we
have

$\displaystyle\dfrac{1}{2}\partialderivative{t}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}^{2}$
$\displaystyle\leq Ch^{m}\|\mathbf{v}_{t},\
\phi_{t}\|_{\mathcal{H}^{m}}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}$
$\displaystyle+Ch\|\mathbf{v},\
\phi\|_{\mathcal{H}^{m}}\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}.$

Dividing by $\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}$
and integrating in time, we get

$\displaystyle\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{E}$
$\displaystyle\leq\|\mathbf{v}_{h}(\boldsymbol{x},0)-\mathbf{v}_{h}^{P}(\boldsymbol{x},0),\phi_{h}(\boldsymbol{x},0)-\phi_{h}^{P}(\boldsymbol{x},0)\|_{E}+$
$\displaystyle+Ch^{m}\int_{0}^{t}\|\mathbf{v}_{t},\
\phi_{t}\|_{\mathcal{H}^{m}}\mathrm{d}\tau+Ch\int_{0}^{t}\|\mathbf{v},\
\phi\|_{\mathcal{H}^{m}}\mathrm{d}\tau.$

Specifying the initial condition for problem (32) to be
$\left(\mathbf{v}_{h}^{P}(\boldsymbol{x},0),\phi_{h}^{P}(\boldsymbol{x},0)\right)$,
it holds
$\|\mathbf{v}_{h}(\boldsymbol{x},0)-\mathbf{v}_{h}^{P}(\boldsymbol{x},0),\
\phi_{h}(\boldsymbol{x},0)-\phi_{h}^{P}(\boldsymbol{x},0)\|_{E}=0$. Taking the
maximum over $[0,T]$, we have

(46)
$\|\mathbf{v}_{h}-\mathbf{v}_{h}^{P},\phi_{h}-\phi_{h}^{P}\|_{\infty,E}\leq
Ch\|\mathbf{v},\ \phi\|_{W^{1,\infty},\mathcal{H}^{m}}.$

By putting together (36), Lemma 2 and (46) we end up proving (44).

######  Remark 4.

We have proved linear convergence under $h$ refinement for the semi-
discretization (25) using projections as in Assumption 4. This is the case of
the projections $\pLLM[1]$ and $\pLLM[2]$. However, in Section 7 we
investigate numerically this error bound, and show high order rates of
convergence. There is numerical evidence that
$\left((Id-\Pi^{1})(c\mathbf{v}_{h}-c\mathbf{v_{h}^{P}}),\nabla\phi\right)\approx
h^{p}$, while
$\|(Id-\Pi^{1})(c\mathbf{v}_{h}-c\mathbf{v}_{h}^{P})\|_{\mathbf{L}^{2}(\Omega)}$
depends linearly on the mesh size $h$.

##  6 Implementation

In this section we introduce the matrix form associated to (30), giving
further details about the computation of the projections. Let us start from
equation (30), that must hold for all $\mathbf{w}_{h}\in X_{h}^{1}$. Consider
as test functions $\\{\mathbf{b}_{i,h}\\}_{i=1}^{N}$, the basis functions of
$X_{h}^{1}$, which are the push-forward with the Piola transformation map
$\iota^{1}$ of the B-splines on the parametric domain. To assemble the
matrices involved in (30) we compute the projections of the basis functions

$\Pi^{1}(c\mathbf{b}_{i,h})=\sum_{l=1}^{N}\theta^{i}_{l}\mathbf{b}_{l,h},$

and we denote by
$\boldsymbol{\theta}^{i}=(\theta^{i}_{1},\dots,\theta^{i}_{N})^{T}$, the
column vectors with the coefficients of the projections
$\Pi^{1}(c\mathbf{b}_{i,h})$ into the space $X^{1}_{h}$. We then define the
matrix $\tilde{\mathbf{A}}\in\mathbb{R}^{N\times N}$ that encapsulates the
second term of (30), and with the previous notation each entry of the matrix
can be computed as

$\begin{split}[\tilde{\mathbf{A}}]_{i,j}&\coloneqq\int_{\Omega}\text{div}\left(\Pi^{1}(c\mathbf{b}_{i,h})\right)\text{div}\left(\Pi^{1}(c\mathbf{b}_{j,h})\right)\mathrm{d}\mathbf{x}\\\
&=\sum_{l=1}^{N}\sum_{m=1}^{N}\theta^{i}_{l}\theta^{j}_{m}\int_{\Omega}\text{div}(\mathbf{b}_{l,h})\text{div}(\mathbf{b}_{m,h})\mathrm{d}\mathbf{x}=(\boldsymbol{\theta}^{i})^{T}\mathbf{A}\boldsymbol{\theta}^{j},\end{split}$

for $i,j=1,\ldots,N$, with $\mathbf{A}\in\mathbb{R}^{N\times N}$ defined as
$[\mathbf{A}]_{l,m}\coloneqq\int_{\Omega}\text{div}(\mathbf{b}_{l,h})\text{div}(\mathbf{b}_{m,h})\mathrm{d}\mathbf{x}$.
It is therefore convenient to store the coefficients of the projectors in the
matrix
$\mathbf{\Theta}=[\boldsymbol{\theta}^{1}|\dots|\boldsymbol{\theta}^{N}]\in\mathbb{R}^{N\times
N}$, from which we obtain
$\tilde{\mathbf{A}}=\mathbf{\Theta}^{\top}\mathbf{A\Theta}$.

For the computation of the projections, we make use of the commutativity
property (18), from which we obtain the two equivalent expressions

$\displaystyle\iota^{1}(\Pi^{1}(c\mathbf{b}_{i,h}))$
$\displaystyle=\hat{\Pi}^{1}(\iota^{1}(c\mathbf{b}_{i,h}))=\hat{\Pi}^{1}(\hat{c}\hat{\mathbf{b}}_{i,h})=\sum_{l=1}^{N}\hat{\theta^{i}_{l}}\hat{\mathbf{b}}_{l,h},$
$\displaystyle\iota^{1}(\Pi^{1}(c\mathbf{b}_{i,h}))$
$\displaystyle=\iota^{1}(\sum_{l=1}^{N}{\theta^{i}_{l}}{\mathbf{b}}_{l,h})=\sum_{l=1}^{N}{\theta^{i}_{l}}\iota^{1}({\mathbf{b}}_{l,h})=\sum_{l=1}^{N}{\theta^{i}_{l}}\hat{\mathbf{b}}_{l,h},$

where $\hat{\theta}^{i}_{l}$ are the coefficients of the projection
$\hat{\Pi}^{1}(\hat{c}\hat{\mathbf{b}}_{i,h})$, and
$\hat{c}=c\circ\mathcal{F}$. Notice that
$\hat{\theta}^{i}_{l}=\theta^{i}_{l}$, and therefore the projections can be
computed in the parametric domain, exploiting the tensor-product structure.
The details on the computation of these projections are given in Appendix A.
We finally recall that, in order to commute with the divergence operator, the
computation of the projection involves (10), or (11) for the periodic case,
both of them requiring the application of a quadrature formula.

Finally, by introducing the notation $\underline{\mathbf{v}}^{n}$ and
$\underline{\phi}^{n}$ for the coefficients of the unknown fields
$\mathbf{v}_{h}^{n}$ and $\phi_{h}^{n}$ respectively, equation (30) can be
written in the following matrix form:

$\left(\mathbf{M}+\frac{k^{2}}{4}\mathbf{\Theta}^{T}\mathbf{A}\mathbf{\Theta}\right)\underline{\mathbf{v}}^{n+1}=\left(\mathbf{M}-\frac{k^{2}}{4}\mathbf{\Theta}^{T}\mathbf{A}\mathbf{\Theta}\right)\underline{\mathbf{v}}^{n}-k\mathbf{\Theta}^{T}{\mathbf{B}}\underline{\phi}^{n},\quad\text{for
}n=0,\dots,N-1,$

where the matrix $\mathbf{M}\in\mathbb{R}^{N\times N}$ denotes the mass matrix
for the space $X^{1}_{h}$, and $\mathbf{B}\in\mathbb{R}^{N\times M}$ is
defined as
$[\mathbf{B}]_{i,j}=\int_{\Omega}b_{j,h}\text{div}(\mathbf{b}_{i,h})\mathrm{d}\mathbf{x}$.

######  Remark 5.

Notice that the computation of the involved matrices is independent of time,
that is, $\mathbf{M},\mathbf{A},\mathbf{B}$ and $\mathbf{\Theta}$ can be
computed only once at the beginning of our method. If the coefficient $c$ is
time dependent, it is necessary to recompute only $\mathbf{\Theta}$ at every
time step.

##  7 Numerical results

In this section, we have numerically investigated the approximation properties
of the method by conducting academic tests to determine the convergence order
under $h$-refinements. Furthermore, concerning the employment of the quasi-
interpolant in Section 4.2 within the described numerical method, even if it
is tested on a simplified model problem, we have supplemented the array of
numerical tests conducted by Holderied et al. [21]. Indeed, their studies
observed the favorable approximation properties of the quasi-interpolant but
not of the overall global method. Lastly, we have verified the conservation of
the total energy for this method. All the numerical tests have been performed
in Matlab, with the isogeometric analysis open source package GeoPDEs [33].

(a)  (b)

Figure 6:  (a) Initial data and numerical solutions of (29) at $T=1$ with
${c}=sin(2\pi{x}_{1})sin(2\pi{x}_{2})+2$ for Dirichlet homogeneous boundary
conditions. (b) $h$-convergence rate estimation for both quasi-interpolant
(Q.I.) and Galerkin (G) methods.

###  7.1 Dirichlet boundary conditions

We consider the wave equation as presented in Section 2. The domain is
$\Omega=\mathcal{F}(\widehat{\Omega})$, where $\widehat{\Omega}=[0,1]^{2}$ and
$\mathcal{F}$ is the NURBS map describing the quarter of a ring, as in Figure
1. We assume a full Dirichlet boundary with homogeneous conditions, that is
$\Gamma_{D}=\partial\Omega$ and $g=0$ in (1). The time interval is $[0,T]$,
where we fix $T=1$. We consider a space dependent coefficient
$c(x_{1},x_{2})\coloneqq sin(2\pi{x}_{1})sin(2\pi{x}_{2})+2$, which is smooth,
with bounded derivatives and bounded away from zero.

We assume that the solution is of the form
$u(\boldsymbol{x},t)=\chi(\boldsymbol{x})Re(e^{i\omega t})$, with $\chi\in
H^{1}_{0}(\Omega)$ solution of

(47)
$\int_{\Omega}c^{2}\nabla\chi\cdot\nabla\Phi\mathrm{d}\boldsymbol{x}=\omega^{2}\int_{\Omega}\chi\Phi\mathrm{d}\boldsymbol{x},\quad\forall\Phi\in
H^{1}_{0}(\Omega).$

Since $\chi$ is not known explicitly, we use as reference solution its
approximation with splines of degree $\mathbf{p}=(6,6)$ in a uniform mesh with
mesh width $h=1/128$. We compute the reference solution, $\chi_{h}$, relative
to the fourth smallest eigenvalue, that is $\omega\approx 4\pi$. Our reference
solution is $u(\boldsymbol{x},t)=\chi_{h}(\boldsymbol{x})Re(e^{i\omega t})$,
and its velocity and pressure fields are respectively
$\mathbf{v}(\boldsymbol{x},t)=c(\boldsymbol{x})\nabla\chi_{h}(\boldsymbol{x})Re(e^{i\omega
t})$ and $\phi(\boldsymbol{x},t)=\chi_{h}(\boldsymbol{x})Re(i\omega e^{i\omega
t})$. In order to study the convergence of our method, we estimate the error
between the numerical solution and the reference one for different uniform
mesh sizes. We let the space mesh width $h$ vary in
$\\{1/8,1/16,1/32,1/64\\}$. We choose a uniform partition $\tau$ of the
interval $[0,T]$ with $k=5\times 10^{-4}$. We fix $p=3$, and we set the
initial data from the reference solution, and they can be seen in Figure
LABEL:fig:5.3a, for $h=1/64$. The plot on the right of Figure LABEL:fig:5.3a
shows the solutions obtained at the final time. For each value of $h$ that we
are considering, we record
$\|\mathbf{v}-\mathbf{v}_{h},\phi-\phi_{h}\|_{\infty,E}$. The plot on the left
of Figure LABEL:fig:5.3b shows the two components of the computed errors, that
are
${\|\mathbf{v}-\mathbf{v}_{h}\|_{\infty,2}}\coloneqq\sup_{t\in\tau}\|\mathbf{v}-\mathbf{v}_{h}\|_{\mathbf{L}^{2}(\Omega)}$
and
${\|\phi-\phi_{h}\|_{\infty,2}}\coloneqq\sup_{t\in\tau}\|\phi-\phi_{h}\|_{L^{2}(\Omega)}$.
We compare these errors with the ones obtained from solving (31), i.e., the
standard Galerkin method. The plot shows third order convergence rates both
for our discretization and for Galerkin. In order to investigate also
convergence with $L^{2}$ norm in time we will measure the errors in the
discrete norm
${\|f(\boldsymbol{x},t)\|_{2,2}}\coloneqq\frac{k}{2}\sum_{t\in\tau\setminus\\{T\\}}(\|f(\boldsymbol{x},t)\|_{L^{2}(\Omega)}^{2}+\|f(\boldsymbol{x},t+k)\|_{L^{2}(\Omega)}^{2})^{\frac{1}{2}}$.
We define the analogous norms for the vector fields. The plot on the right of
Figure LABEL:fig:5.3b shows the same error study under $h$-refinement using
the norm $\|\cdot\|_{2,2}$. We obtain the same error convergence rates, that
are of order $h^{p}$, better then what predicted in the Theorem 2.

(a)  (b)

Figure 7:  (a) Initial data and numerical solutions of (29) at $T=1$ with
${c}=sin(2\pi{x}_{1})sin(2\pi{x}_{2})+2$ for mixed boundary conditions. (b)
$h$-convergence rate estimation for both quasi-interpolant (Q.I.) and Galerkin
(G) methods.

###  7.2 Mixed boundary conditions

In this second example we consider the domain as above. The boundary of
$\Omega$ is defined as in Figure 1, that is $\Gamma_{D}$ is the Dirichlet
boundary, while $\Gamma_{P}$ is the periodic side. The Dirichlet conditions
are homogeneous, that is $g=0$. The time interval is $[0,T]$, with $T=1$. Here
we used the same coefficient for the example in Section 7.1, since it is
periodic in $\Gamma_{P}$. We use separation of variables to construct the
reference solution, as it was done in the example of Section 7.1. The
reference solution is
$u(\boldsymbol{x},t)=\chi_{h}(\boldsymbol{x})Re(e^{i\omega t})$, where
$\omega\approx 4\pi$ and $\chi_{h}$ is the approximation with splines of the
solution of (47) with mixed Dirichlet and periodic boundary condition. In
order to have a fine approximation, we used $\mathbf{p}=(6,6)$ and a uniform
mesh with mesh width $h=1/128$.

As it was done for the previous example, we study the convergence of the
numerical scheme by estimating the error between the numerical solution and
the reference one. We discretize the problem as in (29), choosing $p=3$ and a
uniform mesh. We let the space mesh width $h$ vary in
$\\{1/8,1/16,1/32,1/64\\}$. We choose a uniform partition $\tau$ of the
interval $[0,T]$ with time step size $k=5\times 10^{-4}$. We compute the
initial data from the reference solution, and they can be seen in Figure
LABEL:fig:5.4a, for $h=1/64$. The plot on the right of Figure LABEL:fig:5.4a
shows the solutions obtained at the final time.

As it was done for the example in Section 7.1, for each value of $h$ that we
are considering we compute the errors in the energy norm. We also measured the
errors with the discrete norm $\|\cdot\|_{2,2}$. Figure LABEL:fig:5.4b shows
the two components of the computed errors, compared with the ones obtained
from solving with the standard Galerkin method. Here it seems the error
convergence rate is close to a third order of convergence. Although in the
last refinement step the convergence is reduced, the same behavior is observed
for the solution with the standard Galerkin method.

(a)  (b)

Figure 8:  Global convergence rate estimation, for quasi-interpolant (Q.I.)
and Galerkin (G) methods, for homogeneous Dirichlet boundary conditions (a),
while mixed boundary conditions in (b)

###  7.3 Convergence study under time refinement

So long we discussed the $h$-convergence of the quasi-interpolant method. Here
we check the global convergence rate when refining both the mesh size $h$ and
the time step $k$. We let vary $h$ as before from $1/8$ to $1/64$, this time
taking the time step $k=h$. We show in Figure 8 the convergence of both fields
$\mathbf{v}_{h}$ and $\phi_{h}$ in the same norms introduced above, for the
example of Section 7.1 in Figure LABEL:fig:5.5a, and for the example of
Section 7.2 in Figure LABEL:fig:5.5b. The obtained convergence rate is of
second order, as expected from the use of Crank-Nicolson method.

(a)  (b)

Figure 9:  Energy conservation plots for Dirichlet homogeneous boundary
conditions (a) and mixed boundary conditions (b)

###  7.4 Energy conservation

In Remark 2 we point out that our numerical scheme is preserving the total
energy of the system, as defined in (4). Here we check the energy conservation
for long time simulations, hence we fix $T=300$ and choose two different
uniform partitions $\tau_{1}$ and $\tau_{2}$ of the interval $[0,T]$, the
first with step $k_{1}=0.2$, and the second with step size $k_{2}=0.01$. We
also fix the meshsize $h=1/32$. Recall that we consider as reference solution
a stationary wave with a fixed time-frequency $\omega=4\pi$, and the solution
evolves in time as $\Psi(t)\approx Re(e^{i\omega t})$, therefore it performs
around two complete oscillations for each unit interval. We solve problem (29)
and compute the energy as in (4) for every $t_{n}\in\tau\cap\mathbb{Z}$. In
Figure 9 we show the evolution of the relative errors
$\frac{|E_{0}-E_{t_{n}}|}{|E_{0}|}$, where $E_{t_{n}}$ is the energy evaluated
at time $t_{n}$. Figure LABEL:fig:5.6a shows the evolution of the relative
energy errors in semi-logarithmic scale for the solution of the problem with
homogeneous Dirichlet boundary conditions, as the example in Section 7.1. We
can see that for both time steps the error remains at the level of round-off
errors, with lower numbers for the finer time grid. The same plots are
reproduced in Figure LABEL:fig:5.6b for the solutions of the problem with
mixed boundary conditions of Section 7.2, and we observe a similar behavior.
Notice that on the finer grid, for both Dirichlet and mixed boundary
conditions, we performed $30000$ steps in time without losing energy. The
increasing behavior of the errors for the coarse grid solutions seems only due
to accumulation of round-off errors. We can conclude that the method preserves
the total energy of the system as we expected.

##  8 Conclusions

In this paper we proposed an isogeometric discretization of the anisotropic
wave equation in mixed form, with mixed Dirichlet and periodic boundary
conditions. Our method relies on tensor product projections into spline spaces
with good approximation properties and that commute with the divergence
operator, according to the De Rham complex for splines. The conservation of
the total energy of the system, is imposed weakly by modifying the
differential equations according to the energy constraint and to projections
operators that we introduced. Regarding time discretization we employed Crank-
Nicolson method, which is of second order and energy conservative, though
alternative conservative methods could be selected.

The method employed for energy preservation finds precedent in the field of
magnetohydrodynamics, see for example the discretization proposed by[23] for
the Vlasov-Maxwell equations. What distinguishes our work, except for the
simplified model problem, is the a priori error estimate analysis for the full
method, not just for the projection approximation. In cases where the
introduced projections, in addition to commutativity with the divergence
operator, exhibit $L^{2}$-stability and preserve splines, we assume good
approximation properties of the adjoint projection operator to establish that
the method is of high order - specifically, $h^{m-1}$ assuming the solutions
reside in $H^{m}(\Omega)$ and its vectorial counterpart. However, the
stability requirement in $L^{2}$-norm, is not always guaranteed, as
exemplified by the projections introduced in section 4.2. In light of this,
under the less stringent assumption that the projections are locally stable in
the uniform norm, without requiring approximation properties of the adjoint
projection operators, we prove that the method converges, at least linearly in
$h$.

Numerical tests conducted indicate that, in scenarios where we expect linear
convergence, the method exhibits high order convergence in $h$, matching the
order of convergence achieved through a Galerkin approximation without
projections. This observation suggests that the proposed error estimate might
be improved, representing a potential avenue for future development.

Finally, the conservation of energy is also confirmed by numerical results,
where the relative error remains of the order of machine precision at the
final time.

##  Appendix A Quasi-interpolant projections

(a)  (b)

Figure 10: Graphical visualization of intermediate steps in the computation of
$\mu_{i_{1},i_{2}}$.

###  A.1 Computation of quasi-interpolants

In this section we want to discuss the implementation of $\hat{\Pi}^{1}$, as
given in (16). The construction of the projection in (23) is analogous. Let us
recall that

$\hat{\Pi}^{1}=\pi_{p_{1}}\otimes\pi_{p_{2}-1}^{c}\times\pi_{p_{1}-1}^{c}\otimes\pi_{p_{2}}.$

This is applied to project bivariate vector functions
$\mathbf{f}=\left(f_{1},f_{2}\right)$, that in our practical case are the
functions $\hat{c}\mathbf{\hat{b}}_{i,h}$ as in Section 6. This means that we
need to compute $\pi_{p_{1}}\otimes\pi^{c}_{p_{2}-1}(f_{1})$ and
$\pi^{c}_{p_{1}-1}\otimes\pi_{p_{2}}(f_{2})$. In order to alleviate notation
we fix $p=p_{1}=p_{2}$, and omit the subscript $p$ from the univariate
projections. With some abuse of notation, these tensorizations must be
interpreted as commuting compositions of the kind

$\pi\otimes\pi^{c}(f_{1})=\pi\big{(}\pi^{c}(f_{1})\big{)},$

Notice that $\pi^{c}$ projection is applied to $f_{1}$ seen as a function of
$x_{2}$ for a fixed $x_{1}$, while the $\pi$ projection is applied at the
result of the previous computation by seeing it as a function of $x_{1}$ for
any fixed $x_{2}$, as it is specified in [4, Section 2.2.2]. We focus on
applying $\pi\otimes\pi^{c}(f_{1})$, that is:

$\displaystyle\pi\Big{(}\pi^{c}\big{(}f_{1}(x_{1},x_{2})\big{)}\Big{)}$
$\displaystyle=\pi\left(\sum_{i_{2}=1}^{n_{2}-1}\mu_{i_{2}}(x_{1})\hat{D}_{i_{2}}(x_{2})\right)$
$\displaystyle=\sum_{i_{2}=1}^{n_{2}-1}\pi\big{(}\mu_{i_{2}}(x_{1})\big{)}\hat{D}_{i_{2}}(x_{2})$
$\displaystyle=\sum_{i_{2}=1}^{n_{2}-1}\sum_{i_{1}=1}^{n_{1}}\mu_{i_{1},i_{2}}\hat{B}_{i_{1}}(x_{1})\hat{D}_{i_{2}}(x_{2}),$

and the coefficients $\mu_{i_{1},i_{2}}$ will depend only on $f_{1}$. In order
to compute these coefficients $\mu_{i_{1},i_{2}}$ we may perform the following
steps:

  1. 1.

First evaluate $f_{1}$ on the Cartesian grid given by breakpoints and
midpoints of the first univariate direction, and breakpoints, midpoints and
further midpoints of the second univariate direction, see Figure
LABEL:fig:5.1.A.

  2. 2.

Use these pointwise evaluations to project $f_{1}(x_{1},\cdot)$ with
$\pi^{c}$. This means to perform $\pi^{c}(f_{1}(x_{1},\cdot))$ for every
$x_{1}$ in the set of breakpoints and midpoint of the first univariate
direction. We recall that this requires to use Cavalieri-Simpson composite
quadrature formula. We end up with a set of coefficients
$\\{\mu_{i_{2}}(x_{1})\\}_{i_{2}=1}^{n_{2}-1}$, for every fixed $x_{1}$. In
Figure LABEL:fig:5.1.B, these coefficients are highlighted in blue and green
circles, respectively for breakpoints and midpoints of $x_{1}$ direction, for
the case $p=2$.

  3. 3.

Due to the choice of step 1, for a fixed index $i_{2}$, the function
$\mu_{i_{2}}(x_{1})$ is already evaluated over the breakpoints and midpoints
of $x_{1}$ direction. Therefore we can perform a projection $\pi(\mu_{i_{2}})$
for every $i_{2}=1,\dots,n_{2}-1$. We end up with a set of coefficients
$\\{\mu_{i_{1},i_{2}}\\}_{i_{1}=1,i_{2}=1}^{n_{1},n_{2}-1}$ over the
parametric domain that uniquely identify the bivariate spline that
approximates $f_{1}$.

The procedure to compute $\pi^{c}\otimes\pi(f_{2})$ is analogous, we just need
to swap the evaluation points required per univariate direction. After
projecting the two scalar components $f_{1}$ and $f_{2}$, we have two sets of
coefficients whose joint corresponds to the degrees of freedom of the spline
function that approximates $\mathbf{f}$ in the discrete space
$\hat{X}^{1}_{h}$.

## Acknowledgments

## References

  * [1] D. N. Arnold, R. S. Falk, and R. Winther, Finite element exterior calculus, homological techniques, and applications, Acta numerica, 15 (2006), pp. 1–155.
  * [2] D. N. Arnold and J. J. Lee, Mixed methods for elastodynamics with weak symmetry, SIAM Journal on Numerical Analysis, 52 (2014), pp. 2743–2769.
  * [3] E. Bécache, P. Joly, and C. Tsogka, An analysis of new mixed finite elements for the approximation of wave propagation problems, SIAM Journal on Numerical Analysis, 37 (2000), pp. 1053–1084.
  * [4] L. Beirão da Veiga, A. Buffa, G. Sangalli, and R. Vázquez, Mathematical analysis of variational isogeometric methods, Acta Numerica, 23 (2014), p. 157.
  * [5] L. Beirão da Veiga, A. Buffa, G. Sangalli, and R. Vázquez, An Introduction to the Numerical Analysis of Isogeometric Methods, vol. 9, SEMA SIMAI Springer Series, 07 2016, pp. 3–69, https://doi.org/10.1007/978-3-319-32146-2_1.
  * [6] S. Bertoluzza, The discrete commutator property of approximation spaces, Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, 329 (1999), pp. 1097–1102.
  * [7] D. Boffi, F. Brezzi, M. Fortin, et al., Mixed finite element methods and applications, vol. 44, Springer, 2013.
  * [8] D. Boffi, A. Buffa, and L. Gastaldi, Convergence analysis for hyperbolic evolution problems in mixed form, Numerical Linear Algebra with Applications, 20 (2013), pp. 541–556.
  * [9] A. Bressan and E. Sande, Approximation in FEM, DG and IGA: a theoretical comparison, Numerische Mathematik, (2019).
  * [10] A. Buffa, J. Rivas, G. Sangalli, and R. Vázquez, Isogeometric discrete differential forms in three dimensions, SIAM Journal on Numerical Analysis, 49 (2011), pp. 818–844.
  * [11] R. Dautray and J. L. Lions, Evolution Problems 1, Mathematical Analysis and Numerical Methods for Science and Technology, vol. 5, Springer, Berlin, 2000.
  * [12] C. De Boor, A practical guide to splines, vol. 27, springer-verlag New York, 1978.
  * [13] C. De Boor and G. Fix, Spline approximation by quasiinterpolants, Journal of Approximation Theory, 8 (1973), pp. 19–45.
  * [14] J. A. Evans, Y. Bazilevs, I. Babuška, and T. J. R. Hughes, $n$-widths, sup-infs, and optimality ratios for the $k$-version of the isogeometic finite element method, Computer Methods in Applied Mechanics and Engineering, 198 (2009), pp. 1726–1741.
  * [15] J. A. Evans and T. J. Hughes, Isogeometric divergence-conforming b-splines for the darcy–stokes–brinkman equations, Mathematical Models and Methods in Applied Sciences, 23 (2013), pp. 671–741.
  * [16] J. A. Evans and T. J. Hughes, Isogeometric divergence-conforming b-splines for the steady navier–stokes equations, Mathematical Models and Methods in Applied Sciences, 23 (2013), pp. 1421–1478.
  * [17] J. A. Evans and T. J. Hughes, Isogeometric divergence-conforming b-splines for the unsteady navier–stokes equations, Journal of Computational Physics, 241 (2013), pp. 141–167.
  * [18] L. C. Evans, Partial differential equations, American Mathematical Society, 2010.
  * [19] R. Glowinski and S. Lapin, Solution of a wave equation by a mixed finite element-fictitious domain method, Computational Methods in Applied Mathematics, 4 (2004), pp. 431–444.
  * [20] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numerica, 11 (2002), pp. 237–339.
  * [21] F. Holderied, S. Possanner, and X. Wang, Mhd-kinetic hybrid code based on structure-preserving finite elements with particles-in-cell, Journal of Computational Physics, 433 (2021), p. 110143.
  * [22] T. J. Hughes, J. A. Cottrell, and Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering, 194 (2005), pp. 4135–4195.
  * [23] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, Gempic: Geometric electromagnetic particle-in-cell methods, Journal of Plasma Physics, 83 (2017).
  * [24] A. Kushova, Energy conservative isogeometric techniques for the wave equation, master’s thesis, Università degli studi di pavia, September 2020.
  * [25] B.-G. Lee, T. Lyche, and K. Mørken, Some examples of quasi-interpolants constructed from local spline projectors, Mathematical Methods for Curves and Surfaces: Oslo, (2000), pp. 243–252.
  * [26] T. Lyche and K. Mørken, Spline methods. Draft, Department of Informatics, Center of Mathematics for Applications, University of Oslo, Oslo, (2008).
  * [27] T. Lyche and L. L. Schumaker, Local spline approximation methods, Journal of Approximation Theory, 15 (1975), pp. 294–325.
  * [28] P. Monk et al., Finite element methods for Maxwell’s equations, Oxford University Press, 2003.
  * [29] M. Montardini, F. Remonato, and G. Sangalli, Isogeometric methods for free boundary problems, in Conference on Isogeometric Analysis and Applications, Springer, 2018, pp. 131–155.
  * [30] P. Oswald and B. Wohlmuth, On polynomial reproduction of dual fe bases, in Thirteenth international conference on domain decomposition methods, 2001, pp. 85–96.
  * [31] A. Ratnani and E. Sonnendrücker, An arbitrary high-order spline finite element solver for the time domain maxwell equations, Journal of Scientific Computing, 51 (2012), pp. 87–106.
  * [32] G. Sangalli and M. Tani, Matrix-free weighted quadrature for a computationally efficient isogeometric $k$-method, Computer Methods in Applied Mechanics and Engineering, 338 (2018), pp. 117 – 133.
  * [33] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Computers and Mathematics with Applications, 72 (2016), pp. 523–554.