Léo Morin

SPECTRAL ASYMPTOTICS FOR THE SEMICLASSICAL BOCHNER LAPLACIAN OF A LINE BUNDLE WITH CONSTANT RANK CURVATURE

The goal of this paper is manyfold. Firstly, we want to give a short introduction to the Bochner Laplacian and explain why it acts locally as a magnetic Laplacian. Secondly, given a confining magnetic field, we use Agmon-like estimates to reduce its spectral study to magnetic Laplacians, in the semiclassical limit. Finally, we use this to translate already-known spectral asymptotics for the magnetic Laplacian to the Bochner Laplacian.

-Introduction

-Motivations and context

The spectral theory of the magnetic Laplacian, and the Bochner Laplacian, has given rise to many interesting questions. First motivated by the Ginzburg-Landau theory, bound states of the magnetic Laplacian (ihd+A) * (ihd+A) on a Riemannian manifold in the semiclassical limit h → 0 were studied in many works (see the books [START_REF] Fournais | Spectral methods in surface superconductivity[END_REF][START_REF] Raymond | Bound States of the Magnetic Schrödinger Operator[END_REF]), and appears to have very various behaviours according to the variations of the magnetic field B = dA and the boundary conditions. The first main technique consisted in the construction of approximated eigenfunctions (see for instance the works of Helffer-Mohamed [START_REF] Helffer | Semiclassical Analysis for the Ground State Energy of a Schrödinger Operator with Magnetic Wells[END_REF] and Helffer-Kordyukov [START_REF] Helffer | Semiclassical spectral asymptotics for a magnetic Schrödinger operator with non-vanishing magnetic field[END_REF][START_REF] Helffer | Semiclassical spectral asymptotics for a twodimensional magnetic Schrödinger operator. II The case of degenerate wells[END_REF][START_REF] Helffer | Semiclassical spectral asymptotics for a twodimensional magnetic Schrödinger operator: The case of discrete wells[END_REF]). More recently, an other approach was developped, which consists in an approximation of the operator itself, using semiclassical tools such as microlocalisation estimates and Birkhoff normal forms (As in Raymond-Vu Ngoc [START_REF] Raymond | Geometry and Spectrum in 2D Magnetic wells[END_REF] and Helffer-Kordyukov-Raymond-Vu Ngoc [START_REF] Helffer | Magnetic Wells in Dimension Three[END_REF]). In the semiclassical limit h → 0, we recover the classical behaviour of a particle exposed to the magnetic field B, since the magnetic Laplacian is the quantification of the classical energy.

If we are given a magnetic field B which is not exact, there is no potential A and we cannot define the magnetic Laplacian. However, the Bochner Laplacian 1 p 2 ∆ L p appears to be the suitable generalization in this case, since it acts locally as a magnetic Laplacian. In this context the semiclassical parameter is p = h -1 . Its spectral theory appears to be deeply related to holomorphic structures and to the Kodaira Laplacian (or the renormalized Bochner Laplacian more generaly). For instance, this is exploited in the works of Marinescu-Savale [START_REF] Marinescu | Bochner Laplacian and Bergman kernel expansion of semi-positive line bundles on a Riemann surface[END_REF], and Kordyukov [START_REF] Kordyukov | Semiclassical eigenvalue asymptotics for the Bochner Laplacian of a positive line bundle on a symplectic manifold[END_REF]. In this last paper, the case of non-degenerate magnetic wells with full-rank magnetic field is studied, and expansions of the ground states energies are given, using quasimodes. In [START_REF] Morin | A Semiclassical Birkhoff Normal Form for Symplectic Magnetic Wells[END_REF], similar results were obtained in the special case of the magnetic Laplacian, using a Birkhoff normal form, but also giving a description of semi-excited states, and a Weyl law. In [START_REF] Morin | A semiclassical Birkhoff normal form for constant-rank magnetic fields[END_REF], these result are generalized to constant-rank magnetic fields.

The spectral theory of the Bochner Laplacian is also deeply related to the global geometry of complex manifolds. For instance, in [START_REF] Demailly | Champs magnétiques et inégalités de Morse pour la d"-cohomologie[END_REF], a Weyl law was proven for ∆ L p , and used to get Morse inequalities and Riemann-Roch formulas.

In this paper, we use Agmon-like estimates to reduce the spectral theory of the Bochner Laplacian with magnetic wells to magnetic Laplacians. Then, we deduce spectral asymptotics for the Bochner Laplacian using the results of [START_REF] Morin | A Semiclassical Birkhoff Normal Form for Symplectic Magnetic Wells[END_REF][START_REF] Morin | A semiclassical Birkhoff normal form for constant-rank magnetic fields[END_REF].

-The Bochner Laplacian on a line bundle

Let (M, g) be a compact oriented manifold of dimension d > 1. We consider a complex line bundle L → M over M, endowed with a Hermitian metric h. In other words, we associate to each x ∈ M a 1-dimensional complex vector space L x , and a Hermitian product

h x on L x . L is a d + 1-dimensional manifold such that L = x∈M L x . A smooth section of L (or L-valued function) is a smooth function s : M → L such that s(x) ∈ L x .
It is the generalisation of the notion of function f : M → C, but here the target space can vary with x ∈ M. Similarily, L-valued k-forms are sections of ∧ k T * M ⊗ L. We denote by C ∞ (M, L) the set of smooth sections of L, and Ω k (M, L) the set of smooth L-valued k-forms.

We take ∇ L a Hermitian connexion on (L, h). It is the generalisation of the exterior derivative d. The underlying idea is that the "derivative" of a L-valued function should be L-valued too.

∇ L : Ω k (M, L) → Ω k+1 (M, L) satisfies: (1.1) ∇ L (sα) = ∇ L s ∧ α + sdα, ∀s ∈ C ∞ (M, L), α ∈ Ω k (M, C), (1.2) dh(s 1 , s 2 ) = h(∇ L s 1 , s 2 ) + h(s 1 , ∇ L s 2 ), ∀s 1 , s 2 ∈ C ∞ (M, L).
One can prove that (∇ L ) 2 : Ω 0 (M, L) → Ω 2 (M, L) acts as a multiplication. There exists a real closed 2-form B on M such that:

(1.3) (∇ L ) 2 s = iBs, ∀s ∈ C ∞ (M, L).
Example : The trivial line bundle. The line bundle L = M × C, such that L x = {x} × C is called the trivial line bundle. We identify sections s ∈ C ∞ (M, L) with functions f ∈ C ∞ (M) by s(x) = (x, f (x)). Similarily, L-valued k-forms are identified with C-valued k-forms, and we recover the usual differential objects on M. If L is endowed with the Hermitian product h x ((x, z 1 ), (x, z 2 )) = z 1 z 2 , we call (L, h) the trivial Hermitian line bundle. We write h(z 1 , z 2 ) for short. Hermitian connexions on the trivial line bundle are given by ∇ α = d + iα where α ∈ Ω 1 (M, R) and d is the exterior derivative. The curvature of ∇ α is ∇ 2 α = idα, as shown by the easy but enlightening calculation:

(1.4) ∇ 2 α f = (d + iα)(df + if α) = d 2 f + iα ∧ df + id(f α) + if α ∧ α = iα ∧ df + idf ∧ α + if dα = if dα.
Let us describe now the Bochner Laplacian ∆ L associated to a Hermitian connexion ∇ L on a Hermitian complex line bundle (L, h). First note that the spaces C ∞ (M, L) = Ω 0 (M, L) and Ω 1 (M, L) are endowed with L 2 -norms. The norm of a section s ∈ C ∞ (M, L) is:

(1.5) s 2 = M h x (s(x), s(x))dν g (x),
where dν g denotes the volume form of the oriented Riemannian manifold (M, g). We denote by L 2 (M, L) the completion of C ∞ (M, L) for this norm. The definition of the norm of a L-valued 1-form α is a little more involved. First, using a partition of unity, it is enough to define it for α ∈ Ω 1 (U, L) where U is a small open subset of M. If U is small enough, there exists a section e ∈ C ∞ (U, L) such that h x (e(x), e(x)) = 1. Then for any α ∈ Ω 1 (U, L), there exists a unique X ∈ T M such that α x (•) = g x (X x , •)e x (we identify 1-forms with tangent vectors using the metric g). We define:

(1.6)

α 2 = M g x (X x , X x )dν g (x)
.

The completion of Ω 1 (M, L) for this norm is denoted by L 2 Ω 1 (M, L): it is the space of square-integrable L-valued 1-forms. These norms are associated with scalar products, denoted by brackets ., . .

The formal adjoint of ∇

L : Ω 0 (M, L) → Ω 1 (M, L) for these scalar products is denoted by (∇ L ) * : Ω 1 (M, L) → Ω 0 (M, L). The Bochner Laplacian ∆ L is the self- adjoint extension of (∇ L ) * ∇ L .
It is the operator associated with the quadratic form:

(1.7) Q(s 1 , s 2 ) = ∇ L s 1 , ∇ L s 2 .
We denote by Dom(∆ L ) its domain. C ∞ (M, L) is a dense subspace of Dom(∆ L ) and:

(1.8)

∆ L s 1 , s 2 = ∇ L s 1 , ∇ L s 2 , ∀s 1 , s 2 ∈ Dom(∆ L ).
Since M is compact, one can prove that ∆ L has compact resolvent, and we denote by

(1.9) λ 1 (∆ L ) ≤ λ 2 (∆ L ) ≤ ...
the non-decreasing sequence of its eigenvalues. We will use the following notation for the Weil counting function:

N(∆ L , λ) := ♯{j; λ j (∆ L ) ≤ λ}.
In this paper, we are interested in the semiclassical limit, i.e. the high curvature limit "B → +∞". We can increase the curvature B using tensor products of L. For any p ∈ N, we denote by L p = L ⊗ ... ⊗ L the p-th tensor power of L. L p is still a complex line bundle of M, with

L p x = L x ⊗... ⊗L x . It is endowed with the Hermitian product h p x (s 1 ⊗ ... ⊗ s p , s 1 ⊗ ... ⊗ s p ) = Π p i=1 h x (s i , s i ).
The connexion ∇ L induces a Hermitian connexion ∇ L p on L p by the formula:

∇ L p (s 1 ⊗ ... ⊗ s p ) = (∇ L s 1 ) ⊗ ... ⊗ s p + ... + s 1 ⊗ ... ⊗ (∇ L s p ). The curvature of ∇ L p is (1.10) (∇ L p ) 2 = ipB.
Hence, the high curvature limit is p → +∞. We want to invastigate the behaviour of λ j (∆ L p ) and the corresponding eigensections in the limit p → +∞.

-Main results

One can measure the "intensity" of the curvature B (the magnetic field) in the following way. We denote by B x : T x M → T x M the linear operator defined by

(1.11) g x (B x U, V ) = B x (U, V ), ∀U, V ∈ T x M.
B x is real and skew-symmetric with respect to g x , and thus its eigenvalues lie on the imaginary axis and are symmetric with respect to the origin. We denote its non-zero eigenvalues by:

(1.12)

±iβ 1 (x), • • • , ±iβ s (x),
with β j (x) > 0. Hence, the rank of B x is 2s and might depend on x, but we will soon assume that s is constant, at least locally. The magnetic intensity is the function b : M → R + defined by:

(1.13) b(x) = s(x) j=1 β j (x).
This function is continuous on M, but not smooth in general. However, note that it is smooth on a neighborhood of any point x 0 where the (β j (x 0 )) 1≤j≤s are simple (if s is locally constant near x 0 ).

One of the purposes of this article is to show that the eigensections of ∆ L p are localized near the minimum points of b, and to deduce that the low-lying spectrum of ∆ L p is given by magnetic Laplacians on neighborhoods of the minimal points of b.

We will do the following assumptions.

Assumptions. (A1)

The minimal value of b is only reached at non degenerate points

x 1 , • • • , x N ∈ M. We denote by b 0 = b(x j ) = min x∈M b. (A2) The rank of B is constant on small neighborhoods U 1 , • • • , U N of x 1 , • • • , x N .
We denote by 2s j the rank of B x j .

(A3) We assume b 0 > 0, which is equivalent to say that s j > 0 for any j.

As noticed in several papers, one can prove using Agmon-like estimates that the eigensections of ∆ L p associated to low-lying eigenvalues are exponentially localized near {x ∈ M, b(x) = b 0 }, in the limit p → +∞. Now let us present the local model operators on U j .

Recall that the 2-form B is closed: dB = 0. Hence, if the open sets U j are small enough, B is exact on U j : there exists A j ∈ Ω 1 (U j ) such that B = dA j on U j . We denote by L (j) p the Dirichlet realization of (d + ipA j ) * (d + ipA j ) on L 2 (U j ). It is the self-adjoint operator associated to the following sesquilinear form on C ∞ 0 (U j ):

(1.14) Q j (u, v) = M (du + ipA j u)(dv + ipA j v)dν g .
We prove the following Theorem.

Theorem 1. Let α ∈ (0, 1/2). Under assumptions (A1) and (A3), if η, ε > 0 are small enough, then:

(1.15) λ k (∆ L p ) = λ k L (1) p ⊕ ... ⊕ L (N ) p + O(exp(-εp α )), uniformly with respect to k ∈ [1, K p ],
where

K p = min N(∆ L p , (b 0 + η)p), N(L (1) p ⊕ ... ⊕ L (N ) p , (b 0 + η)p) ,
and N(A, λ) denotes the number of eigenvalues of an operator A below λ, counted with multiplicities.

As a corollary, we can deduce spectral asymptotics for ∆ L p from already-known results for L (j) p . Let us recall some of these results here.

-The full-rank case

Under the assumptions (A1) -(A2) -(A3), we fix a j ∈ {1, • • • , N}, and we denote by B j = dA j . Hence, B j is just the restriction of B to the small open set U j , where it admits a primitive A j . L (j) h is the magnetic Laplacian with Dirichlet boundary conditions on U j , with magnetic field B j . We first focus on the full-rank case, when the rank of B j is maximal: 2s j = d. We define r j ∈ N by the condition

(1.16) ∀n ∈ Z s j , 0 < s j ℓ=1 |n ℓ | < r j ⇒ s j ℓ=1 n ℓ β ℓ (x j ) = 0.
Note that, if the β ℓ (x j ) are pairwise distinct, we can choose r j ≥ 3. Moreover, if the open set U j is small enough we have, for all x ∈ U j and n ∈ Z s j , (1.17) 0 <

s j ℓ=1 |n ℓ | < r j =⇒ s j ℓ=1 n ℓ β ℓ (x) = 0.
The following Theorem was proved in [START_REF] Morin | A Semiclassical Birkhoff Normal Form for Symplectic Magnetic Wells[END_REF].

Theorem 2. We assume (A1) -(A2) -(A3) with 2s j = d and r ≥ 3 in (1.16). Let η, ε > 0 small enough. Then there exists a symplectomorphism ψ : U j → T * R d/2 such that:

(1.18) 1 p 2 λ k (L (j) p ) = λ k n∈N d N [j,n] p + O(p -r j /2+ε ),
uniformly with respect to k ∈ [1, Kp ], where N

[j,n] p is a pseudo-differential operator with principal symbol:

σ(N [j,n] p ) = 1 p s j ℓ=1 (2n ℓ + 1)β ℓ • ψ -1 (x, ξ),
and

Kp = min N(L (j) p , (b 0 + η)p), N(⊕ n N [j,n] p , (b 0 + η)p -1 )
.

Hence, we have a description of the semi-excited states of L (j)

p . In the same paper, Weyl estimates are proven for L (j) p . We directly deduce the following Weyl estimates for ∆ L p . A Similar formula was proven in [START_REF] Demailly | Champs magnétiques et inégalités de Morse pour la d"-cohomologie[END_REF] using a local approximation of the magnetic field by a constant field. Corollary 3. Assume (A 1 ) -(A 2 ) -(A 3 ), and for any j ∈ {1, • • • , N} that s j = d/2, and that (β ℓ (x j )) 1≤ℓ≤N are pairwise distinct. Then, for η > 0 small enough,

(1.19) N(∆ L p , (b 0 + η)p) ∼ p 2π d/2 n∈N d/2 b [n] (x)≤b 0 +η B d/2 (d/2)! , in the limit p → +∞, where b [n] (x) = d/2 ℓ=1 n ℓ β ℓ (x).
Finally, we also deduce asymptotic expansions of the first eigenvalues. and that (β ℓ (x j )) 1≤ℓ≤N are pairwise distinct, and r := min j r j ≥ 5. Then, for any k ∈ N and ε > 0,

λ k (∆ L p ) = b 0 p + r-5 i=0 α i,k p -i/2 + O(p 2-r/2+ε ),
for some coefficients α i,k ∈ R.

Remark. We also have geometric interpretations of the coefficients. First, the full expansion comes from the effective operator N [j,0] p , which is the reduction of L (j) p to the lowest energy of the Harmonic oscillator describing the classical cyclotron motion. Moreover, α 0,k is given by an eigenvalue of an other Harmonic oscillator whose symbol is the Hessian of b at x j (for some 1 ≤ ≤ N): it describes a slow drift of the classical particle arround x j . If the eigenvalues of this oscillator are simple, then a Birkhoff normal form can be used to show that α i,k = 0 if i is odd.

-The constant-rank case

In the non-full-rank case, the kernel of B (which corresponds to the directions of the field lines), has a great influence on the spectrum of ∆ L p . Fix 1 ≤ j ≤ N. If the rank of B j is constant, equal to 2s j , then its kernel as dimension k j = d -2s j . The partial Hessian of b at x j , in the directions of the Kernel of B j , is non-degenerate. we denote by

(1.20) ν 2 j,1 , • • • , ν 2 j,k j
its eigenvalues. For simplicity, we will make the following non-resonance assumptions (however, we can deal with resonances using a resonance order r as in the full-rank case).

Assumptions. (A4) For every j, (β ℓ (x j )) 1≤ℓ≤s j are non-resonnant in the following sense:

∀n ∈ Z s j , n = 0 =⇒ s j ℓ=1 n ℓ β ℓ (x j ) = 0.
(A5) For every j such that k j > 0, (ν j,ℓ ) 1≤ℓ≤k j are non resonnant in the following sense:

∀n ∈ Z k j , n = 0 =⇒ k j ℓ=1 n ℓ ν j,ℓ = 0.
Applying the results of [START_REF] Morin | A semiclassical Birkhoff normal form for constant-rank magnetic fields[END_REF] to get spectral asymptotics for L (j) h , we deduce from Theorem 1 the following corollary.

λ n (∆ L p ) = b 0 p + κp 1/2 + i≥0 α i,n p -i/2 + O(p -∞ ).

Moreover:

• If there is at least one j such that k j = 0, then κ = 0. with ∇e = eiA. Moreover,

• If ∀j ∈ {1, • • • , N}, k j > 0, then κ = min j=1,••• ,N k j ℓ=1 ν j,ℓ .

-Some Remarks

∇ 2 s = ∇e ∧ [(d + iA)u] + ed[(d + iA)u] = e(iA ∧ du) + e(iA ∧ iA)u + ed 2 u + ieudA + edu ∧ iA = ieudA = (idA)s,
and thus B = dA. Hence ∇ acts locally as d+iA, and ∆ L as the magnetic Laplacian (d + iA) * (d + iA). This is the core of Theorem 1.

-On the quantization of a magnetic field

If we are given a closed 2-form B (the magnetic field), the quantization question constist in finding a quantum operator associated to B. If B is exact, this question is answered by the semiclassical magnetic Laplacian ( d+iA) * ( d+iA), with B = dA. Here, > 0 is the semiclassical parameter (Planck's constant) and the semiclassical limit is → 0.

If B is not exact, but if there exists an Hermitian line bundle with Hermitian connexion such that ∇ 2 = iB, then the Bochner Laplacian ∇ * ∇ acts locally as the magnetic Laplacian and hence it is a good candidate. Moreover, we have locally

∆ L p = (d + ipA) * (d + ipA) = p 2 ( 1 p d + iA) * ( 1 p d + iA),
so that the semiclassical parameter is now = 1 p (Also notice the p 2 factor which is important for the eigenvalue asymptotics). The limit → 0 is equivalent to p → +∞ exept that the semiclassical parameter becomes discrete (p ∈ N).

A new question arises : When does such an Hermitian line bundle exists ? Weil's Theorem states that it exists if and only if B satisfies the prequantization condition:

(2.1) [B] ∈ 2πZ,
where [B] denotes the cohomology class of B. This condition also enlightens the discreteness of the semiclassical parameter. Indeed, if one wants to quantize the magnetic field 1 B, then one must have 1 B ∈ 2πZ, and thus 1 p ∈ Z, unless [B] = 0 which means that B is exact (and thus we can use the magnetic Laplacian !).

-Proof of Theorem 1

-Agmon-like estimates

In this section, we prove exponential decay on the eigensections of ∆ L p , away from the set {x 1 , • • • x N }. We need the following Lemma. Lemma 6. There exist p 0 > 0 and C 0 > 0 such that, for p ≥ p 0 and s ∈ C ∞ (M, L),

(1 + C 0 p 1/4 ) ∇ L s 2 ≥ p M (b(x) - C 0 p 1/4 )|s(x)| 2 dx.
Proof. Take a partition of unity (χ α ) α on M, such that suppχ α ⊂ U α with U α a small open subset of M, and

(3.1) 1 = α χ 2 α , α |dχ α | 2 ≤ C.
Then, for any s ∈ C ∞ (M, L) we have

(3.2) ∇ L p s 2 = α ∇ L p (χ α s) 2 - α |dχ α |s 2 ,
and χ α s ∈ C ∞ 0 (U α , L). If every U α are small enough, we can find a non-vanishing section e α on U α which we can use to trivialize the line bundle. Writting χ α s = u α e α for some u α ∈ C ∞ 0 (U α ), we have (3.3)

∇ L p (χ α s) = [(d + ipA α )u α ] e α ,
where A α ∈ Ω 1 (U α ) is such that B = dA α . For the magnetic Laplacian (d + ipA α ) * (d + ipA α ), the desired inequality is well-known: There exist p α ∈ N, C α > 0 such that, for every u α ∈ C ∞ (U α ), and p ≥ p α :

(3.4)

(1 + C α p 1/4 ) (d + ipA α )u α 2 ≥ p Uα (b(x) - C α p 1/4 )|u α (x)| 2 dx.
Using (3.2), (3.3), and (3.4), we deduce that

(3.5) (1 + C 0 p 1/4 ) ∇ L p s 2 ≥ p M (b(x) - C 0 p 1/4 )|s(x)| 2 dx -(1 + C 0 p 1/4 ) α |dχ α |s 2 ,
with C 0 = max α C α , and for p ≥ max α p α . Finally, (3.1) yields

(3.6) (1 + C 0 p 1/4 ) α |dχ α |s 2 ≤ C(1 + C 0 p 1/4 ) s 2 ≤ Cp 3/4 s 2 .
Hence, up to changing C 0 into C 0 + C, Lemma 6 is proved. Now we can use Lemma 6 to prove Agmon-like decay estimates.

Theorem 7. Let α ∈ (0, 1/2), η > 0, and K η = {b(x) ≤ b 0 + 2η}. There exist C > 0 and p 0 > 0 such that, for all p ≥ p 0 and all eigenpair (λ, ψ) of

∆ L p with λ ≤ (b 0 + η)p, M |e d(x,Kη)p α ψ| 2 dq ≤ C ψ 2 .
Proof. Let Φ : M → R be a Lipschitz function. The Agmon formula is:

(3.7) ∆ L p e Φ ψ, e Φ ψ = λ e Φ ψ 2 + dΦe Φ ψ 2 .
Using Lemma 6, we deduce that:

(pb(x) -C 0 p 3/4 -(1 + C 0 p -1/4 )(λ + |dΦ| 2 ))|e Φ ψ| 2 dx ≤ 0.
We split this integral into two parts.

K c η (pb(x) -C 0 p 3/4 (1 + C 0 p -1/4 )(λ + |dΦ| 2 ))|e Φ ψ| 2 dx ≤ Kη (-pb(x) + C 0 p 3/4 + (1 + C 0 p -1/4 )(λ + |dΦ| 2 ))|e Φ ψ| 2 dx
We choose Φ: Φ m (x) = χ m (d(x, K η ))p α , for m > 0, where χ m (t) = t for t < m, χ m (t) = 0 for t > 2m, and χ ′ m uniformly bounded with respect to m. Since Φ m = 0 on K η and pb(x) -C 0 p 3/4 > 0 for p large enough, we have:

K c η (pb(x) -C 0 p 3/4 -(1 + C 0 p -1/4 )(λ + |dΦ m | 2 ))|e Φm ψ| 2 dx ≤ (b 0 + η)p Kη (1 + C 0 p -1/4 )|ψ| 2 dx ≤ Cp ψ 2 . Moreover, since λ ≤ (b 0 + η)p and |dΦ m | 2 ≤ Cp 2α , K c η (pb(x) -C 0 p 3/4 -(1 + C 0 p -1/4 )(b 0 p + ηp + Cp 2α )|e Φm ψ| 2 dx ≤ Cp ψ 2 p K c η (b(x) -(b 0 + η) -C 0 p -1/4 -Cp 2α-1 )|e Φm ψ| 2 dx ≤ Cp ψ 2 ,
for p large enough. But b(x) > b 0 + 2η on K c η , so there is a δ > 0 and p 0 > 0 such that, for p ≥ p 0 :

δ K c |e Φm ψ| 2 dq ≤ C ψ 2 .
Since Φ m = 0 on K, we get a new C > 0 such that:

e Φm ψ 2 ≤ C ψ 2 ,
and we can use Fatou's lemma in the limit m → +∞ to get the desired inequality. Corollary 8. Let ε > 0 and χ : M → [0, 1] be a smooth cutoff function, being 1 on a small neighborhood of

K η + ε = {x; d(x, K η ) < ε}.
Then, for any eigenpair (λ, ψ) of ∆ L p , with λ ≤ (b 0 + η)p we have: We can use these Agmon estimates on ∇ L p ψ to get our second result.

ψ = χψ + O(e -εp α ) ψ ,
(3.10)

∇ L p ((1 -χ)ψ) ≤ (∇ L p χ)ψ + (1 -χ)∇ L p ψ
The first term is dominated by

(3.11) (∇ L p χ)ψ ≤ C (1 -χ)ψ
where χ is a cutoff function such that χ = 1 on K η + ε and χ = 0 on supp(1χ).

We can apply (3.8) to χ to get:

(3.12) (∇ L p χ)ψ ≤ Ce -εp α ψ .

The second term of (3.10) is dominated as in (3.8), using (3.9):

(3.13) (1χ)∇ L p ψ ≤ Cp 1/2 e -εp α ψ .

Finally, (3.10) with (3.12) and (3.13) yields ∇ L p ((1χ)ψ) ≤ Cp 1/2 e -εp α ψ .

3.2 -Comparison of the spectrum of ∆ L p and L (j) p

We recall that the minimum b 0 of b is reached at x 1 , • • • , x N in a non-degenerate way. For η > 0 small enough, the compact set K η = {b(x) ≤ b 0 + η} has N disjoint connected components K (j) η such that x j ∈ K (j) η . We fix the value of η, and we take U j a neighborhood of K (j) η . For ε > 0 sufficiently small, K (j) η + 2ε ⊂ U j .

We denote by B j the restriction of B to U j . L (j) p is the Dirichlet realisation of (d + ipA j ) * (d + ipA j ), with A j ∈ Ω 1 (U j , L) such that B j = dA j . It is the self adjoint operator associated to the quadratic form: We split the proof of Theorem 1 into two Lemmas.

Corollary 4 .

 4 Assume (A1) -(A2) -(A3), and for any j ∈ {1, • • • N} that s j = d/2

Corollary 5 .

 5 Assume (A1) -(A2) -(A3) -(A4) -(A5), and let n ∈ N. Then λ n (∆ L p ) admits a full asymptotic expansion in powers of p -1/2 :

2. 1 -

 1 The Bochner Laplacian and the magnetic Laplacian are locally the same If U is any open subset of M such that there exists a non-vanishing section e ∈ C ∞ (U, L), then any s ∈ C ∞ (U, L) can be written s = ue for u ∈ C ∞ (M). Hence, ∇s = (∇e)u + e(du) = e[(d + iA)u],

  and∇ L p (χψ) = ∇ L p ψ + O(p 1/2 e -εp α ) ψ , uniformly with respect to (λ, ψ).Proof. By Theorem 7, we have:(3.8) (1χ)ψ 2 ≤ (Kη+ε) c |ψ| 2 dq ≤ M e -2εp α |e d(x,Kη)p α ψ| 2 dx ≤ Ce -2εp α ψ 2 ,which gives the first estimates. Moreover, we have with Φ(x) = d(x, K η ),e Φp α ∇ L p ψ ≤ ∇ L p (e Φp α ψ) + p α dΦe Φp α ψ ,and using Agmon's formula 3.7 Theorem 7:∇ L p (e Φp α ψ) 2 = λ e Φp α ψ 2 + p 2α dΦe Φp α ψ 2 ≤ C 2 p ψ 2 . Φp α ∇ L p ψ ≤ Cp 1/2 ψ 2 .

( 3 .

 3 14) Q j (u, v) = U j (d + ipA j )u(d + ipA j )vdx, ∀u, v ∈ H 1 0 (U j ).Let us denote by (3.15) K p = min N(∆ L p , (b 0 + η)p); N N j=1 L (j) p , (0, b 0 + η)p .

 Lemma 9. Let α ∈ (0, 1/2). We have:

uniformily with respect to k ∈ [1, K p ].

Proof. We prove this using the min-max principle. For k ∈ [1, J p ], let ψ k be the normalized eigenfunction associated to λ k (∆ L p ). We will define the quasimode u j,k ∈ C ∞ 0 (U j ) using a trivialisation of L p on U j . Let e j ∈ C ∞ (U j , L) be the nonvanishing local section of L such that, for any u ∈ C ∞ (U j ),

Let χ j ∈ C ∞ 0 (U j ) be a smooth cutoff function, such that χ j = 1 on K

Moreover, by (3.16),

Now, χ = N j=1 χ j satisfies the assumptions of Corollary 8 (with 2ε instead of ε). Thus,

uniformly with respect to k. ψ k being the eigensection associated to λ k (∆ L p ), it remains:

This is true for every

and the Lemma follows from the min-max principle, because the vector space ranged by (u i ) 1≤i≤k is k-dimensional (and p 1/2 e -2εp α = O(e -εp α )).

The reverse inequality is proven similarily.

Lemma 10. Let α ∈ (0, 1/2). We have:

Proof. The k-th eigenvalue of

p is given by an eigenpair

uniformly with respect to k. We define s k = χ k u k e j k , where e j k satisfies (3.16), and we extend s k by 0 outside U j k . Then,

Hence the min-max principle implies

which is the desired inequality.


Maodong Pan
email: mdpan@mail.ustc.edu.cn

Bert Jüttler
email: bert.juettler@jku.at

Angelos Mantzaflaris
email: angelos.mantzaflaris@inria.fr

Efficient Matrix Assembly in Isogeometric Analysis with Hierarchical B-splines

Keywords: Matrix formation, Isogeometric analysis, Hierarchical B-splines, Quasi-interpolation, look-up and sum-factorization, Computational cost

. The new algorithm consists of three stages: approximating the integrals by quasi-interpolation, building three compact look-up tables and assembling the matrices via sum-factorization. A detailed analysis shows that the complexity of our method has the order O(N p 3 ) under a mild assumption about mesh admissibility, where N and p denote the number of degrees of freedom and spline degree respectively. Finally, several experimental results are demonstrated to verify the theoretical results and to show the performance of the proposed method.

Introduction

Introduced to improve the interoperability across the design and analysis pipeline, isogeometric analysis (IgA) [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF][START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURB-S, exact geometry and mesh refinement[END_REF] has been recognized as a promising tool in several disciplines such as nonlinear solid mechanics [START_REF] Elguedj | B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements[END_REF], electromagnetic problems [START_REF] Buffa | Isogeometric methods for computational electromagnetics: B-spline and T-spline discretizations[END_REF], shape optimization [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF] and plate analysis [START_REF] Beirão Da Veiga | An isogeometric method for the Reissner-Mindlin plate bending problem[END_REF]. The core idea of IgA is to employ the same basis functions (e.g., nonuniform rational B-splines (NURBS)) both for the representation of geometry in computer aided design (CAD) and for the approximation of solution fields in finite element method (FEM), allowing from the beginning the elimination of geometry errors [START_REF] Nguyen | Isogeometric analysis: an overview and computer implementation aspects[END_REF]. Another important feature of IgA is that it allows for higher global regularity of the basis functions, up to C p-1 across mesh elements for splines of degree p, which leads to an increased accuracy per degree of freedom [START_REF] Bazilevs | Isogeometric analysis: approximation, stability and error estimates for h-refined meshes[END_REF][START_REF] Beirão Da Veiga | Some estimates for h-p-k-refinement in isogeometric analysis[END_REF] and improved spectrum properties [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF][START_REF] Hughes | Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS[END_REF]. We refer the interested reader to the reference [START_REF] Beirão Da Veiga | Mathematical analysis of variational isogeometric methods[END_REF] for a detailed mathematical analysis of the IgA approach.

However, the use of high-degree basis functions in IgA raises the issue of expensive computational effort required to assemble the system matrices arising in isogeometric Galerkin discretizations [START_REF] Drzisga | The surrogate matrix methodology: Low-cost assembly for isogeometric analysis[END_REF]. For a d-dimensional spline space with spline degree p and N degrees of freedom, standard assemblers based on element-wise Gauss quadrature require O(N p 3d ) operations to compute the system matrix, which is far from being optimal. For this reason, several important contributions have been made towards the development of efficient quadrature rules for tensor-product-based isogeometric discretizations. The first line of research addressing this problem aims at developing reduced or specialized quadrature rules [2, 4-6, 15, 31, 32, 36, 46] that are more efficient than the standard Gauss quadrature with O(p) nodes per knot span. The second category of quadrature rules are the isogeometric collocation methods [START_REF] Auricchio | Isogeometric collocation methods[END_REF][START_REF] Gomez | The variational collocation method[END_REF][START_REF] Schillinger | Isogeometric collocation: Cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations[END_REF], which can be interpreted as one-point quadrature schemes in the IgA context. Another approach aims at exploiting the tensor-product structure, which includes the low-rank tensor approximation [START_REF] Hofreither | A black-box low-rank approximation algorithm for fast matrix assembly in isogeometric analysis[END_REF][START_REF] Mantzaflaris | Low rank tensor methods in Galerkinbased isogeometric analysis[END_REF][START_REF] Scholz | Partial tensor decomposition for decoupling isogeometric Galerkin discretizations[END_REF] and sum-factorization [START_REF] Antolin | Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization[END_REF][START_REF] Bressan | Sum factorization techniques in isogeometric analysis[END_REF][START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF][START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF] techniques. Finally, let us mention two quadrature-free approaches [START_REF] Mantzaflaris | Integration by interpolation and look-up for Galerkin-based isogeometric analysis[END_REF][START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF] that are the starting points of this paper. In the interpolation and look-up approach [START_REF] Mantzaflaris | Integration by interpolation and look-up for Galerkin-based isogeometric analysis[END_REF], the common factor occurring in the integrand is firstly transformed into tensor-product splines, and subsequently the resulting integrals are evaluated exactly using the precomputed look-up tables. Based on this method, the recent work [START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF] further reduced the computational cost of matrix assembly by exploiting the tensor-product structure via the sum-factorization technique. Instead of repeating the description of these prior works, we refer to the papers [START_REF] Calabrò | Quadrature Rules in the Isogeometric Galerkin Method: State of the Art and an Introduction to Weighted Quadrature[END_REF][START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF], whose introductory sections contain an overview of the state of the art.

In the framework of IgA, hierarchical B-splines have attracted increasing attention [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF][START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF][START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF][START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF] since they possess local refinement abilities, which is impossible for tensor-product splines. Unfortunately, similar to the tensor-product case, isogeometric simulations based on hierarchical B-splines face a great computational burden at the point of matrix assembly, especially for high polynomial degrees. However, to the best of our knowledge, so far few works on this topic have been proposed [START_REF] Giannelli | Efficient matrix assembly for hierarchical B-splines[END_REF].

In this paper, we extend the previous work [START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF] to bivariate hierarchical B-splines. Although both methods consist of three stages: spline projection, building look-up tables and assembling the matrices via sum-factorization, performing each step of the new proposed method is more technical and challenging. A detailed complexity analysis demonstrates that the computational cost for matrix assembly has order O(N p 3 ) under a mesh admissibility assumption. Several numerical tests in terms of the N -dependence and p-dependence behavior are provided, which confirm the theoretical result of the new algorithm. In addition, the convergence behavior and accuracy of our method are also verified by applying it to solving a Poisson problem with an adaptive refinement strategy. It should be noted that the extension of the new algorithm to trivariate hierarchical B-splines is not a trivial task, since the last stage -matrix assembly via sum factorization presented in current paper -is not suitable for the trivariate case. This will be discussed in another paper.

The remainder of this paper is structured as follows. Section 2 recalls the basic concepts of hierarchical B-splines and admissible meshes followed by introducing the isogeometric discretizations of an elliptic problem. Section 3-6 describe the outline and three stages of the proposed method. A detailed theoretical analysis of the computational cost of the new algorithm is given in Section 7. Section 8 presents several numerical examples to demonstrate the performance and to verify the theoretical results of our method. Finally, we conclude the paper with a summary and future work in Section 9.

Preliminaries

This section includes three parts. The first part recalls the basic concepts and defines the global index space of hierarchical B-splines (HB-splines). Then the definition and some important properties of the admissible HB-splines are presented. Finally, we use the introduced hierarchical spline space to derive the isogeometric Galerkin discretizations of an elliptic boundary value problem.

Hierarchical B-splines

We consider two sequences (one per coordinate direction, i.e., d = 1 or d = 2) of spaces B d composed of univariate splines of degree p on [0, 1]. The spaces are indexed by the level = 0, . . . , L.

For each level , the two spaces B 1 and B 2 of level are both defined by 2 knot spans. Therefore, they are spanned by bases

B d = {β d,i d (ξ d ) : i d ∈ I } (1)
that consist of (p + 2 ) B-splines

β d,i d (ξ d ).
In order to simplify the notation, we shall employ the global index space Z, where Bsplines of level have indices i d in the subspace

I = {2 + p, . . . , 2 +1 + ( + 1)p -1}.
(

) 2
Note that the index subspaces associated with different levels are mutually disjoint. We denote with level(i d ) the level of an index i d , i.e., i d ∈ I level(i d ) . Additionally, we assume that the knot spans of any level + 1 are created by performing non-uniform dyadic refinement to the ones of the previous level. More precisely, the knot spans of level are split into two knot spans of level + 1, where non-uniform splits are allowed. This assumption implies

• the nestedness of the spline spaces, B d ⊂ B +1 d , • the maximum smoothness C p-1 ,

and

• a level-wise bound on the number of B-splines that possess a support intersection with a given B-spline.

The latter fact is made precise in the following lemma.

Lemma 1. For any given B-spline β d,i d of level , there are at most 2 -(p + 1) + p + 1 B-splines of level that possess a support intersection1 with it.

The tensor-product spline bases B = B 1 ⊗ B 2 consist of all products of univariate B-splines,

B = {β i 1 ,i 2 (ξ 1 , ξ 2 ) = β 1,i 1 (ξ 1 )β 2,i 2 (ξ 2 ) : i 1 ∈ I , i 2 ∈ I }. (3)
They span nested tensor-product spline spaces

B = span B , B ⊂ B +1 . ( 4
)
The arguments of the tensor-product splines are ξ = (ξ 1 , ξ 2 ). After introducing multi-indices i = (i 1 , i 2 ) which is an element of the global index space Z 2 , we can simply write

β i (ξ) = β 1,i 1 (ξ 1 )β 2,i 2 (ξ 2 ), (5)
where B-splines of level have indices i in the subspace I × I . The two consequences of the assumption about dyadic refinement of the univariate spline spaces carry over to the bivariate setting, implying again

• the nestedness of the spline spaces, B ⊂ B +1 , and

• a level-wise bound on the number of B-splines that possess a support intersection with a given B-spline, as described in the following lemma.

Lemma 2. For any given B-spline β i of level , there are at most ( 2 -(p + 1) + p + 1) 2 B-splines of level that possess a support intersection with it.

In order to define hierarchical B-splines (HB-splines), we need to consider a sequence of nested subdomains

[0, 1] 2 = Ω 0 ⊇ Ω 1 ⊇ . . . ⊇ Ω L , (6)
where we assume that the boundary of Ω is aligned to the knot lines of B -1 . The HB-splines

H = {β i : i ∈ H}, (7)
which are defined by the index set

H = {i ∈ L =0 I × I : Ω ⊇ supp β i ⊆ Ω +1 } (8)
with Ω L+1 = ∅, span the hierarchical spline space

H = span H. (9)
Under the definition (8), the indices of the HB-splines H possess the block structure shown in Figure . 1, since only univariate B-splines of the same level are considered in the individual tensor-product bases.

Ensuring admissibility

We recall the notion of admissible mesh configurations from [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF]: The HB-splines H, which are defined by a given subdomain hierarchy, are said to be admissible of class r ≥ 2 if the following condition holds:

supp β i ∩ Ω +r = ∅, ∀i ∈ (I × I ) ∩ H, = 0, 1, . . . , L -r. ( 10
)
Figure 2 illustrates this definition through a simple example of biquadratic hierarchical splines of class 2. 3 and Ω 3 are highlighted in black, red, green and blue respectively.

0 \ Ω 1 , Ω 1 \ Ω 2 , Ω 2 \ Ω
Admissible HB-splines possess three important properties:

• the basis functions in H that take non-zero values over any element of the hierarchical mesh belong to at most r successive levels;

• the number of basis functions acting on each element of the hierarchical mesh does not exceed r(p + 1) 2 , independently of the overall number of hierarchical levels;

• under the assumption of dyadic refinement strategy for creating the sequence of nested spaces, the number of basis functions in H that possess a support intersection with any given HB-spline

β i in H does not exceed r-1 k=1-r p + 1 2 k + 1 + p 2 = O(p 2 ).
These properties play an important role in the analysis of adaptive isogeometric analysis [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF][START_REF] Buffa | Complexity of hierarchical refinement for a class of admissible mesh configurations[END_REF]. In particular, the last one is directly related to the sparsity properties of the isogeometric Galerkin matrices.

In the remainder of this paper we shall assume admissibility of class r = 2. Adaptive mesh refinement subject to this constraint can be performed using the algorithms presented in [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF]. Given an input mesh, the admissibility can be restored by suitably enlarging the nested subdomains Ω . An example is shown in Figure 3.

Isogeometric Galerkin discretizations

To describe the proposed method for assembling the system matrices arising in isogeometric Galerkin discretizations of partial differential equations (PDEs), we consider a general second-order linear elliptic boundary value problem

-div(K∇u) + λ T ∇u + σu = f in Ω, u = 0 on ∂Ω (11)
on the bounded domain Ω ⊂ R 2 . The matrix-valued function K ∈ L ∞ (Ω, R 2×2 ) is symmetric and uniformly positive definite, the scalar function σ ∈ L ∞ (Ω) is non-negative, the vectorvalued function λ and the source function

f satisfy λ ∈ L ∞ (Ω, R 2 ), f ∈ L 2 (Ω).
Considering the Sobolev space V = H 1 0 (Ω), we recall the weak formulation of the problem [START_REF] Buffa | Complexity of hierarchical refinement for a class of admissible mesh configurations[END_REF]

: Find u ∈ V satisfying a(u, v) = F (v) for all v ∈ V, (12)
where the bilinear and linear forms are

a(u, v) = Ω ∇u T (x)K(x)∇v(x) + λ T (x)∇u(x)v(x) + σ(x)u(x)v(x) dx and F (v) = Ω f (x)v(x) dx respectively.
In the framework of isogeometric analysis, the physical domain Ω is parameterized using spline functions. More precisely, we have a diffeomorphism G from an axis-aligned parametric domain Ω

:= [0, 1] 2 satisfying Ω = G( Ω).
For the sake of simplicity, we assume that the parameterization G is given by single-patch hierarchical splines

G(ξ) = i∈H c i β i (ξ), (13)
which are represented as a linear combination of HB-splines β i (ξ) with associated global indices i and coefficients c i ∈ R 2 . See [START_REF] Buchegger | Planar multi-patch domain parameterization via patch adjacency graphs[END_REF][START_REF] Falini | THB-splines multi-patch parameterization for multiply-connected planar domains via Template Segmentation[END_REF] and the references therein for the extension to the multi-patch case. When using the isogeometric framework, the discretization of (12) utilizes the given parameterization of the domain Ω and is obtained by employing the same basis used for defining the parameterization [START_REF] Calabrò | Efficient quadrature rules based on spline quasiinterpolation for application to IGA-BEMs[END_REF]. Specifically, the discretization space V h is spanned by the isogeometric functions

ϕ i = β i • G -1 , β i ∈ H.
Applying the Galerkin approach to the variational problem [START_REF] Buffa | Isogeometric methods for computational electromagnetics: B-spline and T-spline discretizations[END_REF], we arrive at the following discrete formulation: Find

u h in V h satisfying a(u h , v h ) = F (v h ) for all v h ∈ V h . (14)
Solving the discrete variational problem ( 14) requires the formation of several matrices and vectors. More precisely, we need to assemble the stiffness matrix S with the elements

S i,i = Ω ∇ϕ i T K∇ϕ i dx = Ω ∇ ξ β i T ∇ ξ G -1 K∇ ξ G -T ∇ ξ β i | det ∇ ξ G|dξ, (15)
the advection matrix A with the elements

A i,i = Ω λ T ∇ϕ i ϕ i dx = Ω λT ∇ ξ G -T ∇ ξ β i β i | det ∇ ξ G|dξ, (16)
the mass matrix M with the elements

M i,i = Ω σϕ i ϕ i dx = Ω σβ i β i | det ∇ ξ G|dξ, (17)
and the load vector b with the elements

b i = Ω f ϕ i dx = Ω f β i | det ∇ ξ G|dξ. (18)
Note that the hat symbolˆdenotes the pull-back of a function from the physical domain Ω to the parametric domain Ω, i.e.,

K = K • G, λ = λ • G, σ = σ • G, f = f • G.
In order to distinguish it from the gradient on the physical domain, we use the symbol ∇ ξ to denote the gradient in the parametric domain Ω.

Framework overview

In order to keep the presentation concise, we focus on the evaluation of the elements of the stiffness matrix, which can be rewritten as

S i,i = Ω ∇ ξ β i T W ∇ ξ β i dξ, W = | det ∇ ξ G | ∇ ξ G -1 K∇ ξ G -T . ( 19
)
The approach can be easily adapted to other relevant matrices in ( 16), ( 17) and [START_REF] Beirão Da Veiga | Some estimates for h-p-k-refinement in isogeometric analysis[END_REF]. The elements of the kernel W in [START_REF] Beirão Da Veiga | Mathematical analysis of variational isogeometric methods[END_REF] are generally rational functions possessing high degrees. To make the calculation more efficient, we propose a quadrature-free method via spline projection, look-up and sum-factorization techniques. The overall framework of our approach is outlined as follows:

• Step 1. Approximating the integrals via spline projection. We first approximate the elements of W by hierarchical spline functions in H,

W (ξ) ≈ i ∈H W i β i (ξ), ( 20
)
where W i = (w µ,ν i ) µ,ν=1,2 are the associated coefficient matrices. This transforms the stiffness matrix elements [START_REF] Beirão Da Veiga | Mathematical analysis of variational isogeometric methods[END_REF] into sums of integrals of tensor-product B-splines

S i,i ≈ 2 µ=1 2 ν=1 i ∈H w µ,ν i Ω ∂ µ β i ∂ ν β i β i dξ, (21)
where the partial derivatives are defined as ∂ d β = ∂β/∂ξ d . This is the only approximation step in the entire method.

• Step 2. Building compact look-up tables. In order to perform an exact evaluation of the sums [START_REF] Elguedj | B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements[END_REF], we define the following tri-product integrals

T θ,ϑ i d ,i d ,i d = 1 0 β θ d,i d β ϑ d,i d β d,i d dξ d , (22)
where β d,i d , β d,i d and β d,i d are univariate B-splines of degree p (maybe from different levels), and θ, ϑ ∈ {0, 1} represent the order of derivatives. We collect all these values in a look-up table. Consequently, the quantities ( 21) can be rewritten as

S i,i ≈ 2 µ=1 2 ν=1 i ∈H w µ,ν i 2 d=1 T δ(µ,d),δ(ν,d) i d ,i d ,i d , (23)
where δ(•, •) is the Kronecker delta function. The second step of our approach is to generate three compact look-up tables for the exact values of the resulting tri-product integrals [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF].

•

Step 3: Assembling the matrices. In order to obtain a computationally efficient method, we employ sum-factorization to assemble the matrix, based on the previously established look-up tables.

The three steps will be described in more detail in the following three sections.

Approximating the integrals via spline projection

The first step of our approach is to project the matrix-valued function W that appears in the integrals [START_REF] Beirão Da Veiga | Mathematical analysis of variational isogeometric methods[END_REF] into a suitable hierarchical spline space. This can be accomplished by the global fitting or quasi-interpolation methods. We assume that the spline space used for the projector [START_REF] Drzisga | The surrogate matrix methodology: Low-cost assembly for isogeometric analysis[END_REF] is the same as the one for isogeometric discretizations described in Section 2.3. In this case, the size of the required look-up tables in the subsequent steps is relatively small and independent of the mesh size of the HB-splines. Also in the tensor-product case, it was noted that this choice gives optimal performance [START_REF] Mantzaflaris | Integration by interpolation and look-up for Galerkin-based isogeometric analysis[END_REF].

Global fitting

The most straightforward method for approximating a function g defined on Ω with a hierarchical spline s in H is based on global L 2 fitting. This is done by solving the minimization problem

min s∈H Ω(s -g) 2 dξ, ( 24
)
where the unknown spline function

s = j∈H c j β j
depends on the unknown coefficients c j of the HB-splines β j ∈ H. In order to deal with the above problem, one may proceed by using Gaussian quadrature, arriving at the discrete optimization problem min

c j i γ i j∈H c j β j (ξ i ) -g(ξ i ) 2 , (25)
where the points ξ i and positive weights γ i are derived from the Gaussian nodes and weights, respectively. The unknowns c j are determined by the normal equations

X T Xc = X T v (26)
with

X ij = √ γ i β j (ξ i ), v i = √ γ i g(ξ i ), c = (c j ) j∈H . (27)
It is known that the linear system (26) always has a unique solution. The solutions of the two problems ( 24) and ( 25) are identical if the Gauss quadrature exactly integrates the products β i β j and β i g. Otherwise, the solution of the discrete problem is an approximation of the exact one.

Although global fitting is a simple method for spline projection, it requires evaluating the given function and the HB-splines at numerous Gauss points and solving a linear system. In fact, the matrix X T X is equal to the mass matrix [START_REF] Beirão Da Veiga | An isogeometric method for the Reissner-Mindlin plate bending problem[END_REF] of the HB-splines in the case σ| det ∇ ξ G| = 1. We conclude that this simple method is computationally too expensive.

Quasi-interpolation

Quasi-interpolation (QI) is another well-established methodology to construct approximations of a given function by splines. Typically it leads to lower computational costs than global fitting. Recently, the QI method has also been used to deal with the integrals arising in Isogeometric Boundary Element Methods (IGABEM), see [START_REF] Calabrò | Efficient quadrature rules based on spline quasiinterpolation for application to IGA-BEMs[END_REF][START_REF] Falini | An adaptive IgA-BEM with hierarchical B-splines based on quasi-interpolation quadrature schemes[END_REF][START_REF] Falini | A study on spline quasi-interpolation based quadrature rules for the isogeometric Galerkin BEM[END_REF].

Given a function g with domain Ω and a set of blending functions Φ = {φ j : j ∈ J }, a quasi-interpolant of g takes the general form Qg = j∈J λ j (g)φ j [START_REF] Giannelli | Efficient matrix assembly for hierarchical B-splines[END_REF] with certain linear coefficient functionals

λ j (g) = supp φ j g(ξ)ψ j (ξ)dξ (29)
that are determined by suitably chosen weight functions ψ j . The linear operator Q is called the quasi-interpolation operator. This operator is required to reproduce polynomials of degree p, in order to possess good approximation properties. Moreover, certain types of QI operators even reproduce the functions in span Φ, which makes them projectors.

A popular choice for the coefficient functionals is to use linear combinations of sampled values, corresponding to weight functions that are linear combinations of translates of Dirac delta functions. Here we consider only coefficient functionals of this type. QI operators in hierarchical spline space were studied already by Kraft [START_REF] Kraft | Adaptive und linear unabhängige Multilevel B-Splines und ihre Anwendungen[END_REF]. He introduced an HB-spline projector that is constructed in a recursive way. Later, these results have been extended to truncated hierarchical B-splines (THB-splines, see [START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF][START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF]) by Speleers and Manni in [START_REF] Speleers | Effortless quasi-interpolation in hierarchical spaces[END_REF]. In a subsequent paper [START_REF] Speleers | Hierarchical spline spaces: quasi-interpolants and local approximation estimates[END_REF], Speleers showed how to construct general QI operators for hierarchical spline spaces, which are not projectors but still possess good approximation properties. Recently, Giust et al. [START_REF] Giust | Local (T)HB-spline projectors via restricted hierarchical spline fitting[END_REF] established a new quasi-interpolation scheme for both HB-splines and THB-splines by employing local least-square fitting in restricted hierarchical spline spaces. This provides more efficient spline projectors, while obtaining the same accuracy as the previous quasi-interpolants.

The computational costs of QI operators depend on the required number of sampled values (#smpl) and on the number of floating point operations (#flops) needed to evaluate the coefficient functionals. Note that the evaluations of these functionals is perfectly suited for parallel computing. Table 1 summarizes the properties of various quasi-interpolants for bivariate hierarchical splines. It should be noted that these properties are based on certain assumptions the admissibility of the hierarchical meshes.

Building compact look-up tables

After replacing the elements of W in [START_REF] Beirão Da Veiga | Mathematical analysis of variational isogeometric methods[END_REF] by hierarchical spline functions via HB-spline projection, the integrals are transformed into sums of integrals of tensor-product B-splines [START_REF] Elguedj | B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements[END_REF]. The next step of our approach is to evaluate those tri-product integrals, which are defined in the look-up tables [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF], thereby preparing the accurate evaluation of [START_REF] Elguedj | B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements[END_REF].

Owing to the small supports of the B-splines, the majority of the entries in the look-up tables vanish. Non-zero entries are created only by B-spline-triplets ( β d,i d ,β d,i d ,β d,i d ) with overlapping supports. As noted before, we assume that the HB-splines used for isogeometric discretizations and spline projection are the same and that they are admissible of class 2. In this situation, the number of non-zero elements in [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF] is relatively small and the B-splines β d,i d , β d,i d and β d,i d belong to the same level or to two successive levels (see Figure 4).

For uniform B-splines, the tri-product integrals [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF] that do not involve boundary Bsplines2 can be classified into the following three cases:

Λ σ,σ ,σ ij = R βp,σ [0,1,...,p+1] ( ξ) βp,σ [i,i+1,...,i+p+1] ( ξ) βp,σ [j,j+1,...,j+p+1] ( ξ) d ξ (30)
with σ, σ , σ ∈ {0, 1}, σσ σ = 0 and 0 ≤ i, j ≤ p, where βp,σ Ξ denotes the σth order derivative of the B-spline β of degree p and with knot vector Ξ. Note that the case σ = σ = σ = 1 is not considered. These normalized integrals (30) can be precomputed and stored in a look-up table with 7(p + 1) 2 elements. As an example, Table 2 reports these values for B-splines of degree 2.

The integrals [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF] required for assembling the system matrices can be evaluated via

T µ,ν i d ,i d ,i d = 1 h µ+ν-1 d        Λ µ,ν,0 i d -i d ,i d -i d if i d = min{i d , i d , i d } Λ ν,µ,0 i d -i d ,i d -i d if i d = min{i d , i d , i d } Λ 0,µ,ν i d -i d ,i d -i d if i d = min{i d , i d , i d } . ( 31
)
the overall computations. They can be evaluated via Gauss quadrature without compromising the efficiency of the approach.

i, j i, j ΓΞ 0,0,0 We shift them such that the one at level is transformed to the B-spline with global index 2 + ( + 1)p + (p + 1)/2 . We then consider the standardized tri-product integrals

Ξ 0,0,0 ij Ξ 0,0,1 ij Ξ 0,1,0 ij Ξ 0,1,1 ij Ξ 1,0,0 ij Ξ 1,0,1 ij Ξ 1,
i ΓΞ 0,0,0 i3 ΓΞ 0,0,1 i3 ΓΞ 0,1,0 i3 ΓΞ 0,1,1 i3 ΓΞ 1,0,0 i3 ΓΞ 1,0,1 i3 ΓΞ
ij ΓΞ 0,0,1 ij ΓΞ 0,1,0 ij ΓΞ 0,1,1 ij ΓΞ 1,0,0 ij ΓΞ 1,0,1 ij ΓΞ
Λ σ,σ ,σ ij = R βp,σ [ i 2 , i+1 2 ,..., i+p+1 2
] ( ξ) βp,σ

[ p+1 2 , p+1 2 +1,..., p+1 2
+p+1] ( ξ)

• βp,σ

[ j 2 , j+1 2 ,..., j+p+1 2
] ( ξ) d ξ ( 32
)
with σ, σ , σ ∈ {0, 1}, σσ σ = 0 and 0 ≤ i, j ≤ 2 3(p + 1)/2 -1, which can be precalculated and constitute a look-up table with 7(2 3(p + 1)/2 ) 2 entries (see Table 3 for a slice of it). We then evaluate integrals ( 22) via

T µ,ν i d ,i d ,i d = 2 hd µ+ν-1        Λ ν,µ,0 i d -2i d +ρ,i d -2i d +ρ if i d = min{i d , i d , i d } Λ µ,ν,0 i d -2i d +ρ,i d -2i d +ρ if i d = min{i d , i d , i d } Λ 0,µ,ν i d -2i d +ρ,i d -2i d +ρ if i d = min{i d , i d , i d } (33)
with ρ = 2 (p + 1)/2 + p, where hd represents the length of the knot span of B-splines at level + 1 in the d-th coordinate direction.

• Case 3 -Figure 4(c): Two of the B-splines β d,i d , β d,i d and β d,i d belong to the coarser level, the other one belongs to the finer level. Again we assume these B-splines are at level and + 1 respectively, and shift them such that the first B-spline of level (i.e., the one with the lower index) has the global index ( + 1)p + 2 + (p + 1)/2 . Consequently, the integrals ( 22) can be derived from the tri-product integrals

Λ σ,σ ,σ ij = R βp,σ [ i 2 , i+1 2 ,..., i+p+1 2
] ( ξ) βp,σ

[ p+1 2 , p+1 2 +1,..., p+1 2
+p+1] ( ξ)

• βp,σ

[j+ p+1 2 ,j+ p+1 2 +1,...,j+ p+1 2 +p+1] ( ξ) d ξ ( 34
)
via the formula

T µ,ν i d ,i d ,i d = 2 hd µ+ν-1                        Λ µ,ν,0 i d -2i d +ρ,i d -i d if i d ≤ i d < i d Λ µ,0,ν i d -2i d +ρ,i d -i d if i d ≤ i d < i d Λ ν,µ,0 i d -2i d +ρ,i d -i d if i d ≤ i d < i d Λ ν,0,µ i d -2i d +ρ,i d -i d if i d ≤ i d < i d Λ 0,µ,ν i d -2i d +ρ,i d -i d if i d ≤ i d < i d Λ 0,ν,µ i d -2i d +ρ,i d -i d if i d ≤ i d < i d , (35)
where hd and ρ are as in the previous case. The indices in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURB-S, exact geometry and mesh refinement[END_REF] and ( 35) satisfy σ, σ , σ ∈ {0, 1}, σσ σ = 0, 0 ≤ i ≤ 2 3(p + 1)/2 -1, 0 ≤ j ≤ p and ρ = 2 (p + 1)/2 + p respectively, and all the integrals (34) can be precomputed and stored in a look-up table with 7(p + 1)(2 3(p + 1)/2 ) entries. Table 4 presents three slices of this table.

By exploiting the symmetry and compact supports of the B-splines, the look-up tables can be further compressed. Consequently, the duplicated entries and most of the zero entries in Tables 2, 3 and 4 can be omitted. On the other hand, the symmetric structure of these tables simplifies the subsequent look-up procedure.

Our method is also applicable to non-uniform B-splines. In this case, we evaluate all the non-zero tri-product integrals [START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF] in each of the 2 coordinate directions via standard Gauss quadrature, forming 6 look-up tables. To estimate the sizes of these look-up tables, we take the number of levels of the HB-splines into account and use the quantity (2) and Lemma 1 introduced in the preliminary section. Consequently the numbers of nonzero entries in Tables 2, 3 and 4 in each coordinate direction are O(Lp 3 + 2 L+1 p 2 ), O(Lp 3 + 2 L p 2 ) and O(Lp 3 + 2 L p 2 ) respectively. As a difference to the case of uniform B-splines, the sizes of the look-up tables for non-uniform case depend on both the polynomial degree and on the number of levels of the HB-splines.

Efficient matrix assembly

With the established look-up tables at hand, the remaining step of our approach is to perform the sum and product operations in ( 23) more efficiently. To this end, we further employ the sum-factorization technique to accelerate the matrix assembly procedure.

When used in isogeometric analysis, the method of sum factorization [START_REF] Antolin | Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization[END_REF][START_REF] Bressan | Sum factorization techniques in isogeometric analysis[END_REF][START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF][START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF] achieves its efficiency by exploiting the tensor-product structure of multivariate splines. Essentially, a univariate quadrature is performed in each stage of the algorithm realizing the multivariate quadrature. Applying this technique in our framework, we rewrite the sums in [START_REF] Falini | An adaptive IgA-BEM with hierarchical B-splines based on quasi-interpolation quadrature schemes[END_REF] as

S i,i ≈ 2 µ=1 2 ν=1 i 2 ∈H 2 T δ(µ,2),δ(ν,2) i 2 ,i 2 ,i 2 i 1 ∈H 1 [i 2 ] w µ,ν i 1 i 2 T δ(µ,1),δ(ν,1) i 1 ,i 1 ,i 1 = S µ,ν,1 i 1 i 1 i 2 = S i 1 i 2 i 1 i 2 , ( 36
)
where we use the index sets

H 2 = {i 2 ∈ Z | ∃i 1 ∈ Z : (i 1 , i 2 ) ∈ H} and H 1 [i 2 ] = {i 1 ∈ Z | (i 1 , i 2 ) ∈ H}.
Here we introduce two auxiliary sparse tensors S (1)

:= S µ,ν,1 i 1 i 1 i 2 , S (2) := S i 1 i 2 i 1 i 2 .
The sums for all pairs (i, i ) are evaluated in a recursive way: The innermost loop computes and stores the first tensor S (1) with at most O(N p 2 ) non-zero elements, which is subsequently used for evaluating the tensor S (2) with O(N p 2 ) non-zero entries in the second stage of the recursion.

Note that we visit the elements of the index set H in a particular way: The outer loop runs through all indices i 2 that have at least one associated multi-index (i 1 , i 2 ) ∈ H, and the inner loop visits all indices i 1 that are relevant for the given value of i 2 . Consequently, referring to Figure 1 (right), the outer loop runs through all the rows with at least one index in H, while the inner loop visits all elements of the row under consideration. To implement this efficiently, a row-wise ordering of H would be needed.

We will show that it is possible to perform the assembly without such an ordering. In fact, an arbitrary enumeration of the index set H can be used instead. The assembly is realized by the following algorithm:

Algorithm StiffnessMatrixAssembly

Stiffness matrix assembly via sum-factorization

for µ = 1, 2 do for ν = 1, 2 do for i ∈ H do for i 1 ∈ I 1 (i 1 ) do indices of B-splines that overlap β 1,i 1 for i 1 ∈ I 1 (i 1 ) ∩ I 1 (i 1 ) do indices of B-splines that overlap β 1,i1 and β 1,i 1 S µ,ν,1 i1i 1 i 2 = 0
initialization of the first tensor

for i ∈ H do for i 1 ∈ I 1 (i 1 ) do for i 1 ∈ I 1 (i 1 ) ∩ I 1 (i 1 ) do S µ,ν,1 i1i 1 i 2 + = w µ,ν i 1 i 2 T δ(µ,1),δ(ν,1) i1,i 1 ,i 1
evaluation of the first tensor

for i ∈ H do for i ∈ N (i) do indices of B-splines in H that overlap β i S i1i2i 1 i 2 = 0 initialization of the stiffness matrix for µ = 1, 2 do for ν = 1, 2 do for i 2 ∈ I 2 (i 2 ) ∩ I 2 (i 2 ) do indices of B-splines that overlap β 2,i2 and β 2,i 2 S i1i2i 1 i 2 + = T δ(µ,2),δ(ν,2) i2,i 2 ,i 2 S µ,ν,1 i1i 1 i 2
evaluation of the stiffness matrix return S

The index sets I d (i d ) and N (i) are defined by respectively.

I d (i d ) = {i d ∈ I level(i d )-1 ∪ I level(i d ) ∪ I level(i d )+1 : supp(β d,i d β d,i d ) = ∅} and N (i) = {i ∈ H : supp(β i β i ) = 0} (
For simplifying the implementation and reducing the complexity of the above algorithm, we evaluate some extra elements S µ,ν,1 i 1 i 1 i 2 that are not needed for computing the tensor S (2) . The sub-indices i 1 , i 1 of these elements do not belong to the set

H 1 = {i 1 ∈ Z | ∃i 2 ∈ Z : (i 1 , i 2 ) ∈ H}.
Otherwise, we have to check if the indices i 1 , i 1 exist in H 1 in each loop, which results in a non-constant cost per function. Also in the second stage, some elements S µ,ν,1 i 1 i 1 i 2 required for evaluating the final stiffness matrix are not pre-computed in the first stage and are zero by default. The sub-indices i 2 of these elements do not belong to the set H 2 .

Remark 1. Compared with the case of tensor-product B-splines [START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF], the Algorithm Stiff-nessMatrixAssembly for HB-splines becomes more sophisticated due to the non-trivial task of generating the index sets I d (i d ) and N (i) in each loop. However, as analyzed in the next section, the complexity of the new algorithm remains unchanged under a mild assumption.

Complexity analysis

In order to discuss the computational complexity of our algorithm and compare it with Gauss quadrature, we denote the dimension and degree of the hierarchical spline space H by N and p, respectively. The number of non-zero entries in the system matrix, which equals O(N p 2 ) by Lemma 2 and under the assumption of admissible meshes, is a lower bound for the computational costs of matrix assembly.

A certain portion of the computation time is spent for memory operations (e.g., sparse matrix allocation and memory write), and for finding all the basis functions that possess a support intersection with a given B-spline β i in H, i.e., for computing the index sets N (i). This is needed by any approach towards matrix assembly for HB-splines. We do not include these costs in the analysis presented below, since they depend heavily on the data structures used for the sparse matrix representation and on the HB-spline implementation. We will assume that all the index sets N (i), i ∈ H, which have size |N (i)| = O(p 2 ) due to the admissibility assumption, have been precomputed and can be enumerated in O(p 2 ) time. We do, however, include these costs in the experimental part, which will demonstrate that they do not compromise the good behavior of the new method.

We now discuss the complexity of stiffness matrix assembly based on Gauss quadrature and the new method, respectively.

For the standard element-wise Gauss quadrature, one needs to iterate over O(N ) elements and evaluate O(p 4 ) values and derivatives of basis functions at each of the O(q 2 ) Gauss nodes per element (q = p + 1 is a reasonable choice for maintaining the optimal convergence rate [START_REF] Mantzaflaris | Integration by interpolation and look-up for Galerkin-based isogeometric analysis[END_REF]), resulting in a total complexity of O(N p 6 ).

The computational effort of the new proposed method is described by the following theorem.

Theorem 1. We consider an admissible HB-spline discretization of degree p and with N degrees of freedom. As N increases, the number of flops per dof required for assembling the stiffness matrix S (15) via the proposed method tends to O(p 3 ). Proof. Our algorithm, which is outlined in Section 3, consists of three stages: Approximating the integrals via spline projection, building three compact look-up tables and assembling the matrix via sum-factorization.

In the first step, we choose the local (T)HB-spline projector [START_REF] Giust | Local (T)HB-spline projectors via restricted hierarchical spline fitting[END_REF] as the projection operator. As stated in Section 4 (cf. Table 1), this type of quasi-interpolant requires asymptotically O(p 2 ) flops per dof.

We already noticed that the sizes of the built look-up tables for uniform B-splines and non-uniform B-splines are O(p 2 ) and no more than O(Lp 3 + 2 L+1 p 2 ), respectively. The number of Gauss quadrature nodes needed for calculating each of the entries is at most

3p + 1 2 (p + 1) = O(p 2 ),
hence the total numbers of quadrature nodes per dof (i.e., the numbers of flops per dof) required for generating the look-up tables for uniform and non-uniform B-splines are asymptotically equal to lim

N →∞ 1 N • O(p 2 ) • O(p 2 ) = 0 and lim N →∞ 1 N • O(Lp 3 + 2 L+1 p 2 ) • O(p 2 ) = 0,
respectively. Therefore, the computational cost for building the look-up tables is negligible compared to the overall effort. The task of the matrix assembly stage is to compute two auxiliary tensors in a recursive way. To estimate the cost of Algorithm StiffnessMatrixAssembly, it suffices to consider only the non-zero elements and terms in [START_REF] Hughes | Efficient quadrature for NURBS-based isogeometric analysis[END_REF] Remark 2. Although we restricted the presentation to the assembly of stiffness matrices, the advection matrix A (16) and the mass matrix M (17) can be dealt with analogously. The computational complexities of assembling these three kinds of matrices possess the same order, i.e., O(N p 6 ) for Gauss quadrature and O(N p 3 ) for our method.

Numerical experiments

In this section, we demonstrate the performance of the proposed algorithm by solving a 2D Poisson problem with Dirichlet boundary conditions on four different domains, which are illustrated in Figure 5. First, the theoretical complexity of our method in Theorem 1 is verified by exploring the dependence of the computational time on the number of degrees of freedom and on the polynomial degree. We show the superiority of our method by comparing it with the standard Gauss quadrature in terms of the computational time. Second, the convergence behavior for adaptive refinement with the new method is investigated. The algorithms are implemented in C++ using the G+Smo software [START_REF] Jüttler | Geometry + Simulation modules: Implementing isogeometric analysis[END_REF], and all the tests are run on a desktop PC with 32GB of RAM and an octa-core Intel i7 CPU @3.4GHz. For the purpose of obtaining a more accurate result, in each example we take the average of the computational times of 10 repeated experiments. In addition, although we do not consider the costs of memory access and computing the index set N (i) for any basis β i in the complexity analysis, the reported computational time in this section does include them. However, we do not consider the costs of the spline projection (first step of the algorithm) since the available implementation is not yet sufficiently optimized. We feel that it is still justified to present these experimental results since the overall computational complexity is dominated by the last step of the assembly procedure.

N -dependence of the computational time

We firstly investigate the relation between the complexity and the number of degrees of freedom N . The tests are performed on four computational domains (see Figure 5(a)-5(d)) using different HB-splines and polynomial degrees (p = 3, 4, 5).

On each domain, we report the assembly time for six hierarchical meshes with monotonically increasing numbers of degrees of freedom. These meshes are generated by repeatedly performing a dyadic refinement to all the cells of the initial hierarchical meshes depicted in Figure 6. Note that the post-processing step is indispensable so that all the hierarchical meshes are admissible of class 2.

Figure 7(a)-7(d) report the time needed for the system matrix formation on different domains with various N and degrees p. All these log-log plots reveal that the computational time of our method nearly grows linearly as the number of degrees of freedom increases (note that the slope is almost 1), which is consistent with the theoretical result presented in Theorem 1. It also shows that the costs of memory access and computing the index sets N (i) do not have an obvious impact on the N -dependence behavior.

Degree-dependence of the computational time

Next we explore the dependence of the complexity on the spline degree p both for Gauss quadrature and the new method. As suggested in [START_REF] Mantzaflaris | Integration by interpolation and look-up for Galerkin-based isogeometric analysis[END_REF][START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF], we choose p+1 nodes per element for Gauss quadrature, which suffices to obtain the optimal convergence rate.

The tests are conducted on four different computational domains. For each domain, we report the assembly time for a sequence of hierarchical meshes with four levels and spline degrees varying from 2 to 10. For each degree p, we choose the same hierarchical mesh for all the domains, and the corresponding hierarchical mesh is generated in the following way: a tensor-product mesh consisting of 4p × 4p cells is repeatedly refined by selecting the left-bottom 2p × 2p cells of each level to form the subdomain at the next level. Obviously, these meshes are admissible of class 2. As an example, we show the mesh for case p = 2 in Figure 8.

Figure 9(a)-9(d) report the required time per dof of these approaches, which is needed for assembling the system matrix. The slopes depicted in these log-log plots demonstrate that our method achieves the theoretical asymptotic behavior analyzed in Theorem 1. In contrast, the growth of the time of Gauss quadrature with respect to the degree p is lower than predicted by the theoretical analysis, i.e., O(N p 6 ), at least for p ≤ 10. This is probably related to the implementation of Gauss quadrature provided in the G+Smo library [START_REF] Jüttler | Geometry + Simulation modules: Implementing isogeometric analysis[END_REF], which has been particularly optimized for the low degree cases.

Convergence behavior

Besides the N -dependence and degree-dependence behavior, we also verify the convergence performance of the proposed method embedded into an adaptive refinement framework, which is fulfilled by solving a Poisson problem on the quarter annulus-shaped domain (Figure 5(a)) with Dirichlet boundary condition and known exact solution u(x, y) = e -100((x-1) 2 +(y-1) 2 ) .

With a uniform tensor-product B-spline as the initialization, an adaptive strategy for numerically solving the Poisson problem by our method repeats the following steps:

• Solve the discrete problem [START_REF] Calabrò | Quadrature Rules in the Isogeometric Galerkin Method: State of the Art and an Introduction to Weighted Quadrature[END_REF]. The boundary coefficients of the approximate solution u h are determined by least square fitting of the exact solution at boundary points with B-splines, the remaining coefficients are computed by solving a linear system using the preconditioned conjugate gradient method with incomplete Cholesky factorization.

• Estimate the error distribution. An element-based posteriori error estimator is employed and the error on each cell C ∈ Ω is defined in [START_REF] Kraft | Adaptive und linear unabhängige Multilevel B-Splines und ihre Anwendungen[END_REF]. • Mark the cells with large errors and refine them. The criterion for marking cells is introduced below.

• Apply the algorithms presented in [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF] to ensure the admissibility of the hierarchical mesh, which is a prerequisite for the proposed approach.

The procedure runs until the maximal number of levels L max is achieved.

To estimate the error, we consider at least C 1 -continuous discrete solutions u h (i.e., p ≥ 2). For a posteriori error estimator [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF], the estimate ε C of the local error on a cell C ∈ Ω is defined as

ε C (u h ) = h C f + ∆u h L 2 (C) , (38)
where h C denotes the diameter of cell C. A cell C is marked for refinement provided the criterion ε C ≥ τ is satisfied for a threshold τ . We use the second strategy that is discussed in [START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF] for choosing τ , i.e., a fixed percentage (denoted as η) of all cells should be marked in each refinement step. The parameter η is set as 25% in our experiment.

In the numerical test, we start with a uniform mesh consisting of 8×8 elements and repeat the refinement procedure five times, i.e., L max = 5. Figure 10 demonstrates the convergence plots both for quadratic and cubic splines. With respect to the degrees of freedom used to achieve a certain accuracy, the adaptive hierarchical approach is, as expected, superior to a uniform refinement. This is due to the fact that, the exact solution of the Poisson problem is a smooth function with a peak at the location (1, 1), once the error estimates are dominated by the ones around the peak, only few cells around the peak are marked for refinement in the adaptive approach (see the second row of Figure 10). In contrast, the uniform method refines all the cells, which results in an over-refinement of the domain. We also compare the adaptive refinement by the new method with the one by Gauss quadrature, which indicates that the accuracy of our method is maintained provided the spline projection is sufficiently accurate. Indeed, we cannot visually detect any difference between the two methods, which is due to the fact that the spline projections used by our methods are very accurate. For this reason, we provide Table 5 that shows the subtle differences.

Conclusions and future work

Based on the previous paper [START_REF] Pan | Fast formation of isogeometric Galerkin matrices via integration by interpolation and look-up[END_REF], this work introduces an efficient approach for assembling the matrices in IgA with bivariate HB-splines. The novel approach relies on three key ingredients: (1) projecting the common matrix-valued function into a suitable hierarchical spline space via quasi-interpolation; (2) building three compact look-up tables; (3) assembling the matrices by employing the sum-factorization technique. We present a detailed complexity analysis, which shows that the number of flops per dof required for assembling the system matrices via the proposed method tends to O(p 3 ) as N increases. In addition, the Step 0 1 2 3 4 Gauss 7.3226779e-2 2.1718522e-2 1.0762212e-3 1.0607612e-4 1.8194260e-05 Ours 7.3226785e-2 2.1718860e-2 1.0762424e-3 1.0607620e-4 1.8194257e-05 Table 5: The L 2 errors of the solution obtained using adaptive refinement combined with assembly via Gauss quadrature and our method, respectively. Here each number is displayed with 8 digits to highlight the slight differences between these two approaches. accuracy of our method is maintained provided the spline projection is sufficiently accurate. Finally, these theoretical results are verified by several numerical examples.

There are some major problems that are worthy of further investigation. Firstly, since the exposition in current work is restricted to bivariate HB-splines, we will consider the extension of the proposed approach to d-variate (d > 2) case. Already for d = 3, the problem gets more complicated, since the assembly procedure requires another auxiliary (intermediate) tensor, the construction and efficient evaluation of which would need further study, in particular for generalizing our computational complexity estimates.

Secondly, the present implementation has not been fully optimized, especially for the low-degree splines, see Figure 9. We will explore some options to accelerate the code. Thirdly, the extension to THB-splines [START_REF] Giannelli | THB-splines: An effective mathematical technology for adaptive refinement in geometric design and isogeometric analysis[END_REF][START_REF] Giannelli | THB-splines: The truncated basis for hierarchical splines[END_REF] is also a possible topic for future work, which is especially challenging since the factorization into univariate integrals (as used in this paper) is not possible. Finally, for large-scale problems, the solution of the linear system tends to dominate the total computational cost and matrix-free techniques [START_REF] Sangalli | Matrix-free weighted quadrature for a computationally efficient isogeometric k-method[END_REF] might be more efficient. Thus we will consider the extension of the present approach to matrix-free applications.

Figure 1 :

 1 Figure 1: The global indices (coordinates of the solid dots on the right) of the HB-splines (left) of degree 2 with three levels. Here the coordinates of all the dots represent the subspaces I × I ( = 2, 3, 4) of the global index system Z 2 , where I 2 = {8, . . . , 13}, I 3 = {14, . . . , 23} and I 4 = {24, . . . , 41}.

Figure 2 :

 2 Figure 2: An admissible mesh of class 2 for biquadratic HB-splines with N = 3. The subdomains Ω0 \ Ω 1 , Ω 1 \ Ω 2 , Ω 2 \ Ω3 and Ω 3 are highlighted in black, red, green and blue respectively.

Figure 3 :

 3 Figure 3: An admissible mesh (right) of class 2 for biquadratic HB-splines obtained from a given input mesh (left).

3 Figure 4 :

 34 Figure4: The B-splines in the tri-product integrals[START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF] belong to at most two successive levels. They all belong to the same level (a), two of them belong to the finer level (b), or two of them belong to the coarser level (c). Note that the knots of the B-splines in (22) may be different from the ones shown in these figures.

  This can be confirmed by exploiting the properties of B-splines plus index shifting and swapping operations. Here h d represents the length of the knot span of the basis functions β d,i d , β d,i d and β d,i d . • Case 2 -Figure 4(b): Two of the B-splines β d,i d , β d,i d and β d,i d belong to the finer level, the remaining one belongs to the coarser level. Suppose these B-splines are at level and + 1 respectively.

  37) respectively, where I -1 = ∅ and I L+1 = ∅. According to Lemma 1 and 2, the sizes of these index sets (I d (i d ) and N (i)) do not exceed (p + 1)/2 + 6(p + 1) = O(p) and 10(p + 1) 2 + ( (p + 1)/2 + (p + 1)) 2 = O(p 2 )

1 i 1 i 1 i 2 ,

 12 . The first loop, which evaluates the first tensor, computes 4 • N • O(p) • O(p) non-zero elements S µ,ν,requiring 2 flops per element. The computation of the second tensor with 4 • N • O(p 2 ) non-zero elements proceeds by summing over O(p) indices i 2 and needs 2 flops per term. Thus the cost of the matrix assembly step tends to O(p 3 ) per dof as N increases. Summing up, the computational effort per dof for the proposed algorithm asymptotically amounts to O(p 2 ) + 0 + O(p 3 ) = O(p 3 ).

Figure 5 :

 5 Figure 5: Four different computational domains used in the experiments. Each domain is parameterized by single-patch B-splines, and the corresponding iso-parametric curves are demonstrated.

Figure 6 :

 6 Figure 6: The initial hierarchical meshes for different domains. The left resp. right mesh is used for domain (a) and (c) resp. (b) and (d).

Figure 7 :

 7 Figure 7: N -dependence of the assembly time on four computational domains with different spline degrees (p = 3, 4, 5). (a)-(d) correspond to domain (a)-(d) (cf. Figure 5(a)-5(d)) respectively.

Figure 8 :

 8 Figure 8: An admissible mesh of class 2 for biquadratic HB-splines with four levels.

Figure 9 :

 9 Figure 9: p-dependence of the assembly time of Gauss quadrature and the proposed method on different computational domains. (a)-(d) correspond to domain (a)-(d) (cf. Figure 5(a)-5(d)) respectively.

  degree p = 2 (b) degree p = 3

Figure 10 :

 10 Figure 10: Convergence of error in L 2 norm with respect to degree p = 2 (left) and p = 3 (right) for the quarter annulus-shaped domain. The first row shows the error plots for uniform refinement by Gauss quadrature, adaptive refinements by Gauss quadrature and our method. Note that the latter two curves are visually indistinguishable. The second and last row present the hierarchical meshes on the computational domain obtained via adaptive refinement by Gauss quadrature and our method, respectively.

Table 1 :

 1 Comparison of various hierarchical spline QI schemes for the bivariate case with respect to their applicability to HB-or THB-splines, the spline projector property, and the computational costs per degree of freedom (dof).

	Method	applicable to	Spline	#smpl	#flops
		HB	THB	projector	per dof
	Kraft [38]	Yes	No	Yes	O(p 2 )	O(p 2 )
	Speleers and Manni [49]	No	Yes	Yes	O(p 2 )	O(p 2 )
	Speleers [48]	No	Yes	No	O(1)	O(p 2 )
	Giust et al. [29]	Yes	Yes	Yes	O(1)	O(p 2 )

Table 2 :

 2 The look-up table for tri-product integrals of quadratic uniform B-splines with respect to Case 1.

Table 3 :

 3

	1,1,0
	i3

A slice of the look-up table for tri-product integrals

[START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF]

of quadratic uniform B-splines with respect to Case 2, scaled by the factor Γ = 6, 720.

Table 4 :

 4 Three slices of the look-up table for tri-product integrals[START_REF] Evans | Discrete spectrum analyses for various mixed discretizations of the Stokes eigenproblem[END_REF] of quadratic uniform B-splines with respect to Case 3, scaled by the factor Γ = 26, 880.

	1,1,0
	ij

These B-splines β d,j d satisfy supp(β d,i d β d,j d ) = ∅.

The remaining integrals, which involve boundary B-splines, only account for a relatively small portion of

Acknowledgement

This work was supported by the Austrian Science Fund through the project "Geometry + Simulation" (NFN S11708), and by the European Research Council via the project "CHANGE" (694515). Maodong Pan was also supported by the Natural Science Foundation of China (No. 61972368) and the USTC Research Funds of the Double First-Class Initiative (No. YD0010002003). These supports are gratefully acknowledged.