Properties of Spline Spaces Over Structured Hierarchical Box Partitions

TL;DR

This paper explores how local refinements in hierarchical spline spaces can optimize the number of basis functions and improve numerical stability in computational applications.

Contribution

It demonstrates that local one-directional refinements in LRB spline spaces can reduce basis functions and enhance matrix condition numbers.

Findings

Refinements can minimize B-splines per element to polynomial degree limits.

Local modifications improve the condition numbers of mass and stiffness matrices.

LRB spaces can be constructed to contain and extend THB spline spaces.

Abstract

Given a spline space spanned by Truncated Hierarchical B-splines (THB), it is always possible to construct a spline space spanned by Locally Refined B-splines (LRB) that contains the THB-space. Starting from configurations where the two spline spaces are equal, we adress what happens to the properties of the LRB-space when it is modified by local one-directional refinement at convex corners of, and along edges between dyadic refinement regions. We show that such local modifications can reduce the number of B-splines over each element to the minimum prescribed by the polynomial bi-degree, and that such local refinements can be used for improving the condition numbers of mass and stiffness matrices.

Equations

$$\begin{array}{rl}B& :=B[\vec{x}]B[\vec{y}],\\ Q& :=B[\vec{s}]B[\vec{t}],\end{array}$$

$$\displaystyle\begin{split}\mathinner{B\bigr{\rvert}}_{\beta}(x,y)&=\frac{(\ell-x)^{3}(\ell-y)^{3}}{\ell^{6}},\\ \mathinner{Q\bigr{\rvert}}_{\beta}(x,y)&=\frac{(\ell-x)^{3}(\ell-y)^{3}}{36\ell^{6}}\,.\end{split}$$

$${\int }_{\beta }{B}^2=\frac{{\ell }^2}{49},{\int }_{\beta }{Q}^2=\frac{{\ell }^2}{49}\cdot \frac{1}{36^2}.$$

$$\beta ={\text{bigtimes}}_{i=1}^d{J}_i.$$

$$\mathcal{S}_{p}^{\mu}(\mathcal{M})\coloneqq\Big{\{}f\colon\Omega\to\mathbb{R}\mathrel{\mathop{\ordinarycolon}}\mathinner{f\rvert}_{\beta}\in\Pi_{p}\text{ for all }\beta\in\mathcal{E}\\ \text{ and }f\in C^{p_{k}-\mu(\gamma)}\text{ for all $k$-mesh-rectangles }\gamma\in\mathcal{M},\\ \text{ with }k=1,\ldots,d\Big{\}}.$$

$$\mathcal{M}=\left\{\gamma \cap {\Omega }_{\ell }:\gamma \in {\mathcal{M}}_{\ell }\text{ for }\ell =1,\dots ,M\right\},$$

$$B[\vec{t}]=\{\begin{array}{ll}1,& x\in [t_1,t_2);\\ 0,& \text{otherwise}.\end{array}$$

$$B[\vec{t}](x)=\frac{x-t_1}{t_{p+1}-t_1}B[{\vec{t}}^{-}](x)+\frac{t_{p+2}-x}{t_{p+2}-t_2}B[{\vec{t}}^{+}](x),$$

$${\vec{t}}^{+}=(t_2,\dots ,t_{p+2}),{\vec{t}}^{-}=(t_1,\dots ,t_{p+1}).$$

$$B[{\vec{t}}^1,\dots ,{\vec{t}}^d](\vec{x})=\prod _{i=1}^{d}B[{\vec{t}}^i](x_i),$$

$$\mathrm{supp}(B[{\vec{t}}^1,\dots ,{\vec{t}}^d])=\overline{\left\{\vec{x}\in {\mathbb{R}}^d:B[{\vec{t}}^1,\dots ,{\vec{t}}^d](\vec{x})\ne 0\right\}}.$$

$$\dim ({\Pi }_p(\beta ))=\prod _{i=1}^d(p_i+1).$$

$${\mathcal{M}}_{i+1}={\mathcal{M}}_i+{\gamma }_i$$

$$\mathcal{M}={\mathcal{M}}_N\supseteq {\mathcal{M}}_{N-1}\supseteq \cdots \supseteq {\mathcal{M}}_2\subseteq {\mathcal{M}}_1$$

$$f=\sum _{B_i\in {\mathcal{B}}_{\ell +1}}{c}_i^{\ell +1}(f)B_i,$$

$${\mathrm{trunc}}^{\ell +1}B:=\sum _{\begin{array}{c}{B}_i\in {\mathcal{B}}_{\ell +1}\\ \mathrm{supp}{B}_i⊈{\Omega }_{\ell +1}\end{array}}{c}_i^{\ell +1}(B)B_i.$$

$${\mathrm{trunc}}^{\ell +1}B=B-\sum _{\begin{array}{c}{B}_i\in {\mathcal{B}}_{\ell +1}\\ \mathrm{supp}{B}_i\subseteq {\Omega }_{\ell +1}\end{array}}{c}_i^{\ell +1}(B)B_i$$

$${\mathcal{H}}_{\ell +1}^{\mathrm{trunc}}$$

$${\mathcal{H}}_{\ell +1}^{\mathrm{new}}$$

$${\mathcal{H}}_{\ell +1}$$

$$\Vert e{\Vert }_{L^2}=\Vert u-f{\Vert }_{L^2}.$$

$${\int }_{\Omega }uvd\Omega ={\int }_{\Omega }fvd\Omega ,$$

$$\vec{M}\vec{c}=\vec{b},$$

$${\vec{M}}_{ij}={\int }_{\Omega }{\phi }_i{\phi }_{j}d\Omega ,{\vec{b}}_j={\int }_{\Omega }f{\phi }_{j}d\Omega .$$

$$\Delta u=f\text{ in }\Gamma ,$$

$$u=0\text{ on }{\Gamma }_D,\frac{\partial u}{\partial \vec{n}}=g\text{ on }{\Gamma }_N.$$

$${\int }_{\Omega }\nabla u\nabla vd\Omega ={\int }_{\Omega }fvd\Omega -{\int }_{{\Gamma }_N}gvdS,$$

$$\vec{A}\vec{c}=\vec{b},$$

$${\vec{A}}_{ij}={\int }_{\Omega }\nabla {\phi }_i\cdot \nabla {\phi }_{j}d\Omega ,{\vec{b}}_j={\int }_{\Omega }f{\phi }_{j}d\Omega -{\int }_{{\Gamma }_N}{\phi }_{j}gdS.$$

$$\mathrm{Cond}(\vec{\mathrm{B}}):=\Vert \vec{B}\Vert \Vert {\vec{B}}^{-1}\Vert ,$$

$$\mathrm{Cond}(\vec{\mathrm{B}})=\frac{|{\lambda }_1(\vec{B})|}{|{\lambda }_n(\vec{B})|}.$$

Full text

1 Introduction

The use of Hierarchical B-splines (HB) introduced in forseyHierarchicalBsplineRefinement1988 has gained much attention in Isogeometric Analysis (IgA) in recent years. Hierarchical B-splines are based on a dyadic sequence of grids determined by scaled lattices. On each hierarchical level a spline space is defined as the tensor product of univariate spline spaces spanned by uniform B-splines.

Hierarchical B-splines do not constitute a partition of unity, a much desired property in both Computer Aided Design (CAD) and IgA. As a remedy to this Truncated Hierarchical B-splines (THB) THBsplinesTruncatedBasis2012 ; VUONG20113554 were introduced, where B-splines on one hierarchical level are suitably truncated by B-splines from finer hierarchical levels when the support of a B-spline at a finer level is contained in the support of a B-spline at a coarser level.

An alternative to the THB-approach for forming a partition of unity came with the introduction of Locally Refined B-splines (LRB) dokkenPolynomialSplinesLocally2013 , where initial tensor product B-splines are split until only B-splines of minimal support remain. LRB permits dyadic refinement of hierarchical meshes while ensuring that all B-splines have minimal support. In the occasional case where meshlines at a dyadic level are too short to split an LR B-spline, the meshlines in question are extended. This fact ensures that the spline space spanned by THB-splines is either identical to or constitutes a subset of the LRB spline space.

In IgA open knot vectors are used to simplify the interpolation of boundary conditions, as reported for THB in GIANNELLI2016337 and for LRB in johannessenIsogeometricAnalysisUsing2014 . In open knot vectors the multiplicity at boundary knots is set to $m=d+1$. An alternative approach is to use B-splines with knot multiplicity of $m=1$ along the boundary. In order to force the partition of unity in this case, a ghost domain is added around the domain of interest, as seen in johannessenSimilaritiesDifferencesClassical2015a for both THB and LRB. This distinction is illustrated for univariate cubic splines in Figure 1.

In Section 2 we address the effects the choice of boundary knot multiplicity has on condition numbers. To distinguish between single multiplicity and open knots at the domain boundary we prefix any method using single knot multiplicity on the boundary with a ghosted domain by an “S”. Using this naming convention, the methods addressed in johannessenIsogeometricAnalysisUsing2014 are respectively S-THB and S-LRB. In this paper, we show that for the same tensor product spline space, THB and LRB are superior with respect to condition numbers of mass and stiffness matrices compared to respectively S-LRB and S-THB. We also explain the intriguing near constant behaviour of the condition numbers reported in johannessenSimilaritiesDifferencesClassical2015a , where S-LRB and S-THB were addressed, and condition numbers seemed to be nearly independent of the refinement level. We show that this is due to single knot multiplicity at domain boundaries for the examples presented in johannessenSimilaritiesDifferencesClassical2015a . Further it is showed that for more levels of refinement the condition numbers for the mass matrix for S-THB and S-LRB will meet and then follow the growing curves for respectively THB and LRB.

In HB and THB the refinement procedure (at an element level) consists of marking elements for splitting. Marked elements are subsequently split in both parameter directions. This contrasts with the refinement procedure LRB allows, namely that of splitting an element in a single parameter direction at the time, provided that at least one B-spline is split in the process. This can be used to modify the hierarchical refinement, and possibly improve the approximation properties of the resulting spline space. In the remaining sections we use open knot vectors at domain boundaries and address how such modifications influence the condition numbers for mass and stiffness matrices for bi-cubic spline spaces. The remaining sections are structured as follows:

Section 3

gives a lightweight introduction to box-partitions and spline spaces over such partitions. The starting point for THB and LRB refinement is a tensor product spline space. The key concept of element overloading is defined, the situation where more B-splines cover an element than are needed for spanning the polynomial space over the element. We briefly summarize some key properties. Subsequently, we recall the definitions of both LRB and THB splines. We also relate the refinement strategies for LRB to T-splines.

Section 4

takes a look at overloading. We look at how to reduce or completely remove overloaded regions in a mesh. We showcase some specific overloading patterns that occur for hierarchical refinement of THB and LRB. Furthermore, we show how local modifications to the LRB-mesh reduce overloading as well as condition numbers.

Section 5

provides a quantitative comparison between the methods. We conduct our numerical experiments using modified central and diagonal refinement examples from johannessenSimilaritiesDifferencesClassical2015a with a finer initial tensor-product mesh. This gives enough room on each hierarchical level for the local modifications to take place. The examples show that LRB with no overloading have smaller condition numbers for the mass matrix per degree of freedom than THB and LRB. However, the difference seems to be so small that in general all methods have a similar behaviour.

** Section 6**

summarizes the main results of this paper.

2 Condition numbers and knotline multiplicity at domain boundary

In johannessenSimilaritiesDifferencesClassical2015a , hierarchical refinement was performed for five levels of refinement, using S-THB and S-LRB. The results reported that the number of refinement levels had little to no influence on the evolution of condition numbers for stiffness and mass matrices. There were some minute differences between S-THB and S-LRB, but they followed the same trend. In Figure 2 we display the condition number of the mass matrix for up to eight refinement levels for S-THB, S-LRB, THB and LRB when run on a hierarchical mesh from johannessenSimilaritiesDifferencesClassical2015a . The relevant mesh at the fifth refinement level is displayed in Figure 4.

The results from johannessenSimilaritiesDifferencesClassical2015a is reproduced, and corresponds to the S-THB and S-LRB curves for the first five refinements. However, at the sixth refinement, the curve for the condition number of the mass matrix for S-LRB breaks off and grows exponentially following the curves of LRB that starts three orders of magnitude lower. In Figure 2 there are also two additional curves (S-LRB1 and LRB1). These are added to show that modifying the mesh by inserting additional knot lines in one parameter direction, with the effect of reducing overloading, significantly reduces the condition numbers of LRB-refinement. This modified mesh is shown in Figure 4b. We will discuss such modifications more closely in Section 4.

Multiplicity of domain boundary knot lines also influences the condition number of the stiffness matrix, as seen in Figure 3. Here we see that the condition numbers for single boundary knot multiplicity (S-THB, S-LRB and S-LRB1) are two orders of magnitude higher than the condition numbers for open knot vectors (THB, LRB, LRB-1).

2.1 Boundary knotline multiplicities

We now take a stab at explaining the drastic change in behaviour occuring at $n=6$ refinements for the S-LRB and S-THB methods as shown in Figure 2. Since the condition number of a matrix are computed in terms of its largest and smallest eigenvalues, we decided to take a look at the geometric localization of the eigenvectors corresponding to these eigenvalues. By coloring the hierarchical mesh based on the influence of each in terms of the corresponding coefficient in the eigenvector, we obtained a rudimentary geometric visualization of these eigenvectors. In Figure 5, we see the smallest eigenvector for the mass-matrix corresponding to LRB and S-LRB at the first, third and sixth refinement, and in Figure 6, the corresponding largest eigenvector.

These figures correspond to the behaviour observed in Figure 2 where the conditioning for LRB grows after only one refinement, whereas S-LRB needs six refinements before the behaviour in the refined region is registered.

The reason for this behaviour is due to the size of the B-splines defined along the boundary in comparison to the size of the B-splines defined in the interior of the domain. In order to illustrate this, we compute analytically the entries in the mass matrix corresponding to B-splines on various tensor product level and compare these values to the mass matrix entry corresponding to a B-spline defined in the corner of the domain.

Observation for the mass matrix

Over the domain $\Omega=[0,1]\times[0,1]$ we define a tensor product grid with element size $\ell$. In the case of bi-cubic spline spaces, the B-spline defined in the lower left corner of the domain can for LRB and S-LRB be written in terms of their knots as

[TABLE]

where $\vec{x}=\vec{y}=[0,0,0,0,\ell]$ and $\vec{s}=\vec{t}=[-3\ell,-2\ell,-\ell,0,\ell]$. In both cases, the two B-splines have only one element of support in the domain $\Omega$, namely $\beta\coloneqq[0,\ell]\times[0,\ell]$. To get a feel for the differences in influence on the mass matrix these B-splines have, we compute the corresponding diagonal elements in the mass matrix.

The polynomial restrictions of $B$ and $Q$ to the element $\beta$ is

[TABLE]

In other words, $\mathinner{B\rvert}_{\beta}=36\mathinner{Q\rvert}_{\beta}$. If we now compute the diagonal mass matrix entries corresponding to these two elements, we obtain the following:

[TABLE]

We here see that the matrix element corresponding to the corner B-spline $Q$ arising in S-LRB is three orders of magnitude smaller than the matrix element corresponding to B. Recall the disparity between the curves in Figure 2, where the differences in the condition numbers also were three orders of magnitude.

Remark 1

This relation between $Q$ and $B$ is both dimension and degree-dependent. The effect will be magnified for higher spatial dimension and higher polynomial degree.

3 Box partitions, meshes, and spline spaces

In this section we review the concept of box partitions and spline spaces over such partitions. While the construction generalizes to any dimension, we will gradually focus our attention to the two-dimensional case, as this is most relevant for our discussion. A fully general treatment can be found in dokkenPolynomialSplinesLocally2013 ; patriziLinearDependenceBivariate2018 . The fundamental building block of a box partition is the $d$-dimensional box.

Definition 1

A box $\beta$ in $\mathbb{R}^{d}$ (or $d$-box) is the Cartesian product of $d$ closed finite intervals $J_{1},\ldots,J_{d}$:

[TABLE]

The dimension of $\beta$ is defined to be the number of non-trivial intervals in its definition, and is denoted $\dim(\beta)$. We call a $d$-box of dimension $d$ an element, while a $d$-box of dimension $d-1$ is called a mesh-rectangle. To any mesh-rectangle, we associate an integer $k$ corresponding to which parametric dimension its trivial component resides, and we call the mesh-rectangle a $k$-mesh-rectangle if this has to be emphasized.

In the two-dimensional setting $(d=2)$, a meshline is a one-dimensional mesh-rectangle.

Remark 2

Note that these naming-conventions are independent of the dimension of the ambient space. Hence, a mesh-rectangle may very well be something different from a rectangle. As an example, a mesh-rectangle in $\mathbb{R}^{4}$ is an axis aligned box. Furthermore, the integer $k$ corresponding to any mesh-rectangle encodes the direction of the mesh-rectangle. In the two-dimensional case, where mesh-rectangles are lines, a $1$-mesh-rectangle is a vertical line, and a $2$-mesh-rectangle is a horizontal line.

As customary in discretizations of computational domains, a large domain is partitioned into a set of non-overlapping smaller geometrical entities. We call such a partition in this specific setting a box partition, and this is more precisely defined as follows:

Definition 2

Let $\Omega\subset\mathbb{R}^{d}$ be an element ($d$-box of dimension $d$). A finite collection $\mathcal{E}$ of elements is said to be a box partition of $\Omega$ if

${\kern 0.0pt\beta}^{\mathrm{o}}_{i}\cap{\kern 0.0pt\beta}^{\mathrm{o}}_{j}=\emptyset\text{ for all }\beta_{i},\beta_{j}\in\mathcal{E}\text{ where }\beta_{i}\neq\beta_{j}$, and 2. 2.

$\displaystyle\bigcup_{\beta\in\mathcal{E}}\beta=\Omega\,$.

In other words, a box partition is an interior-disjoint partition of $\Omega$ into a set of smaller elements.

Associated to any element $\beta$ is its boundary, which naturally consists of boxes of dimension one less, i.e., mesh-rectangles. Given a box partition of a larger element $\Omega$, it is therefore sensible to discuss the set of mesh-rectangles associated to this box partition.

Definition 3 (Informal)

Given a box partition $\mathcal{E}$ of a domain $\Omega$, we may naturally associate a set of mesh-rectangles $\mathcal{M}$ called a box mesh on $\Omega$ formed by taking unions of element boundaries.

Remark 3

The link between a box partition $\mathcal{E}$ and the associated box mesh $\mathcal{M}$ is such that by knowing one of them you may recover the other. The box mesh generated by a box partition is denoted $\mathcal{M}(\mathcal{E})$, and the box partition generated by a box mesh is denoted $\mathcal{E}(\mathcal{M})$.

As our ultimate goal is to define spline spaces based on tensor-product splines over box-partitions, we need to have a concept of knot multiplicity in this more general setting.

Definition 4

A box mesh with multiplicity is a pair $(\mathcal{M},\mu)$ where $\mu\colon\mathcal{M}\to\mathbb{N}$ associates to each mesh-rectangle $\gamma$ a positive integer $\mu(\gamma)$, called the multiplicity of the mesh-rectangle. Note that this is completely analogous to the notion of knot multiplicity for univariate B-splines.

Definition 5

Let a polynomial multi-degree $p=(p_{1},\ldots,p_{d})$ as well as a box mesh with multiplicity $(\mathcal{M},\mu)$ corresponding to the box-partition $\mathcal{E}$ of a $d$-dimensional domain $\Omega$ be given. The spline space of degree $p$ over $\mathcal{M}$ is defined as

[TABLE]

A dimension formula for general spline spaces over box partitions was presented in pettersenDimensionMultivariateSpline2013 . In general, the dimension depends on both the topological properties of the box partition and the parametrization of the box partition. In the two-dimensional case — with some requirements on the length of the constituent meshlines — the formula reduces to a formula depending only on the topological features of the mesh. We consider this outside the scope of this text, and refer the reader to pettersenDimensionMultivariateSpline2013 for details.

In order to compute with spline spaces over box partitions of the form above, we must be able to construct a set of basis functions that span this space. Several constructions has been studied. We will only be dealing with Truncated Hierarchical B-splines, and Locally Refined B-splines.

Before we move on, we define the notion of a hierarchical mesh, a type of box partition over which spline bases of the aforementioned type may be defined. This provides a common ground for comparison of the two methods. The construction is simple and relies on marking regions for which a tensor product mesh of various refinement levels is used.

Definition 6

Let $\Omega$ be a domain, and let $\mathcal{M}_{1}\subset\cdots\subset\mathcal{M}_{M}$ be a sequence of nested tensor product meshes on $\Omega$. Let $\Omega_{1}\supset\cdots\supset\Omega_{M}$ be a set of nested subsets of $\Omega$ whose boundaries $\partial\Omega_{\ell}$ align with the meshlines of the corresponding mesh on the coarser level $\mathcal{M}_{\ell-1}$ for $\ell=2,\ldots,M$. The hierarchical mesh $\mathcal{M}$ is defined as

[TABLE]

i.e., $\mathcal{M}$ consists of meshlines from each level intersected with the corresponding region, see Figure 7.

3.1 Tensor product splines

The foundation for all the locally refined spline spaces over box partitions addressed in this paper, is the tensor product B-spline. We start by glossing over some preliminary definitions.

Recall that a univariate B-spline of polynomial degree $p$ relies on exactly $p+2$ knots. This observation enables us to define B-splines locally without referring to some global knot vector.

Definition 7

Given a polynomial degree $p$ and a non-decreasing knot-vector $\vec{t}=(t_{1},\ldots,t_{p+2})$, we recursively define the univariate B-spline $B[\vec{t}]\colon\mathbb{R}\to\mathbb{R}$ as follows:

If $p=0$, then

[TABLE]

If $p>0$, then

[TABLE]

where $\vec{t}^{+}$ and $\vec{t}^{-}$ are obtained by dropping the first and last elements of $\vec{t}$ respectively:

[TABLE]

In the cases of a vanishing denominator, the whole term is taken to be zero.

Such univariate splines can be easily extended to higher dimensions through a tensor product construction.

Definition 8

Let the polynomial multi-degree $p=(p_{1},\ldots,p_{d})$ and the $d$ local knot vectors $\vec{t}^{1},\ldots,\vec{t}^{d}$ be given. The $d$-variate tensor product B-spline $B[\vec{t}^{1},\ldots,\vec{t}^{d}]\colon\mathbb{R}^{d}\to\mathbb{R}$ is then defined as

[TABLE]

where $\vec{x}=(x_{1},\ldots,x_{d})$.

The support of $B[\vec{t}^{1},\ldots,\vec{t}^{d}]$ is the closure of the area where the B-spline takes non-zero values, which we denote by:

[TABLE]

Since our B-spline construction is inherently local, we need to know when a tensor product B-spline has minimal support with respect to some box mesh.

Definition 9 (Informal)

A B-spline $B=B[\vec{t}^{1},\ldots,\vec{t}^{d}]$ has support on $(\mathcal{M},\mu)$ if all the knot lines of $B$ occurs as meshlines in $\mathcal{M}$. We say that $B$ has minimal support on $(\mathcal{M},\mu)$ if in addition, all the knot lines of $B$ occur consecutively in $(\mathcal{M},\mu)$.

One of the central concepts we will be addressing in this paper is the overloading of elements. We make this precise in the following definition.

Definition 10

Let a box partition $\mathcal{E}$ of a domain $\Omega$ and a polynomial multi-degree $p=(p_{1},\ldots,p_{d})$ be given. Assume that we construct a set $\mathcal{B}$ of B-splines degree $p$ over the mesh $\mathcal{M}$ corresponding to $\mathcal{E}$. We say that an element $\beta$ is overloaded with respect to $\mathcal{B}$ if the number of B-splines with support on $\beta$ is larger than the dimension of the corresponding space of polynomials over this element, namely

[TABLE]

We now proceed to review the definitions of LR B-splines and THB-splines.

3.2 Locally refined spline spaces

In preparation for the following discussion, we will adopt the notational convention as in johannessenSimilaritiesDifferencesClassical2015a in order to differentiate between the distinct types of basis functions. Depending on the underlying box partition, some of these types may coincide.

[TABLE]

Furthermore, in this and the following sections we will be dealing with box partitions and spline spaces in $\mathbb{R}^{2}$, unless otherwise stated.

LR-splines

Locally Refinable Splines (LRB or LR-splines) was introduced by dokkenPolynomialSplinesLocally2013 . The LR-spline framework permits the insertion of local splits in a tensor product mesh, and subsequently enables local refinement of the mesh. Being scaled tensor product B-splines, LR-splines admit a set of nice properties. The set of LR B-splines form a partition of unity. Their scaling weights are positive, meaning that they satisfy the convex hull property, and are therefore inherently stable in computations. Moreover, with some restrictions on the refinement process, linear independence of the resulting set of functions can be guaranteed.

LR-splines are defined over so-called LR-meshes, being special box partitions. Starting from an initial tensor product mesh, meshlines are inserted sequentially, yielding a sequence of box-meshes, where no meshline is allowed to terminate in the middle of an element. This is formalized in the following definition, and Figure 8 gives an example.

Definition 11

An LR mesh is a box mesh $\mathcal{M}=\mathcal{M}_{N}$ resulting from a sequence of meshline insertions in an initial tensor product mesh $\mathcal{M}_{1}$. That is

[TABLE]

for $i=1,\ldots,N-1$ where each intermediate mesh is a box mesh.

Remark 4

We often think of an LR-mesh as a sequence of intermediate meshes

[TABLE]

as each intermediate step is needed for the LR B-spline construction.

Over such an LR-mesh we may define the associated set of LR B-splines algorithmically. Starting from an initial space of tensor product B-splines, meshlines are inserted sequentially. Whenever a meshline completely traverses the support of a B-spline, the B-spline is split according to the knot insertion procedure, and two new B-splines are added. The B-spline that was split is removed.

Definition 12

Let $\mathcal{M}$ be an LR-mesh over a domain $\Omega$ and $p=(p_{1},p_{2})$ a polynomial bi-degree. We define the set $\mathcal{L}(\mathcal{M})$ of LR B-splines of degree $p$ over $\mathcal{M}$ algorithmically as in 1.

Remark 5

Note that all LR B-splines have minimal support on the resulting mesh. This is by construction. However, there is an important distinction to be made, namely that the set of LR B-splines differ from the set of minimal support B-splines on the resulting mesh. This is due to the LR refinement procedure putting some restrictions on the resulting mesh. A survey on the properties of LR-splines and minimal support B-splines are given in patriziLinearDependenceBivariate2018 .

Truncated hierarchical B-splines

Hierarchical B-splines, first introduced in forseyHierarchicalBsplineRefinement1988 , is a method for specifying locally refined spline spaces on hierarchical meshes. Recall that a hierarchical mesh consists of regions corresponding to various levels of tensor product grids. The hierarchical B-spline construction involves replacing any B-spline with support completely contained in a region of a finer level by B-splines at this finer level. This procedure will, however, lead to coarse B-splines partially overlapping the finer regions, and does not constitute a partition of unity.

A remedy to this problem came with the introduction of Truncated Hierarchical B-splines THBsplinesTruncatedBasis2012 , where B-splines on a coarse level is truncated by B-splines on a finer level. This leads to the resulting set of B-splines forming a partition of unity. The construction relies on the truncation operator. Recall that a spline $f\in\mathrm{span}(\mathcal{B}_{\ell})$ can be represented in terms of the finer basis $\mathcal{B}_{\ell+1}$:

[TABLE]

where $c^{\ell+1}_{i}(f)$ is the coefficient multiplying $B_{i}$ in the representation of $f$ in terms of $\mathcal{B}_{\ell+1}$. For uniform B-splines, this relation is often called the two-scale relation. The truncation operator is defined as follows:

Definition 13

Let $B\in\mathcal{B}_{\ell}$ be a coarse B-spline. The truncation with respect to the set of fine B-splines $\mathcal{B}_{\ell+1}$ and the corresponding region $\Omega_{\ell+1}$ is

[TABLE]

Remark 6

The definition above represents the truncation operator in an additive sense, where the contributions from the finer level are summed up. It is also possible to represent the truncation operator subtractively, by instead removing the bits of the representation that have been replaced by finer B-splines:

[TABLE]

Definition 14

Let $\mathcal{M}$ be a hierarchical mesh over a domain $\Omega$ (see Definition 6) and $p=(p_{1},p_{2})$ a polynomial bi-degree. On each level $\ell=1,\ldots,N$ we have a tensor product spline space $V_{\ell}$ spanned by a collection of B-splines $\mathcal{B}_{\ell}=\mathcal{B}(\mathcal{M}_{\ell})$. We define the set of THB-splines of degree $p$ over $\mathcal{M}$ algorithmically as in Algorithm 2.

Remark 7

Note that in the cases where a B-spline $B\in\mathcal{B}_{\ell}$ is truncated with respect to $\mathcal{B}_{\ell+1}$ and $\Omega_{\ell+1}$ and the support of $B$ happen to be entirely contained in $\Omega_{\ell+1}$, the truncation operator completely removes the coarse B-spline. In the THB-spline construction, this has the effect of replacing the coarse B-splines with fine B-splines defined in its support.

Remark 8

A simple framework for the implementation of truncated hierarchical B-splines is given in garauAlgorithmsImplementationAdaptive2018a , and this serves as a good introduction to the many ways such splines have been implemented in the literature. Efficient algorithms for the assembly of finite element matrices are also presented.

T-splines

While not directly addressed in this paper, we briefly mention T-splines as LRB with local modifications to the LR-meshes used in this paper is able to reproduce the spline space generated by semi-standard T-splines Sederberg-2004 and Analysis Suitable T-splines Scott-2012 . An example of an Analysis Suitable T-mesh in the index domain is displayed in Figure 9 to the left, with the corresponding LR-mesh to the right. This is a close up of the structure of a mesh similar to the one used in Figure 11b.

4 Local modification of meshes and the reduction of overloading

In this section we take a deeper look at the overloading of elements, and how local modifications to the mesh may be used to remedy this. Recall from the previous definition that an element $\beta$ in a box-partition $\mathcal{E}$ is said to be overloaded if the number of supported B-splines on the element exceeds the number needed to span the full polynomial space over the element.

We are interested in such overloaded regions, because by reducing or removing completely the overloading on elements we may

reduce the bandwidth of the resulting finite element matrices; and 2. 2.

improve conditioning of finite element matrices.

Such overloaded regions occur for LRB in convex corners of a fine hierarchical level, where a large B-spline from one hierarchical level overlaps several elements of a finer hierarchical level. For THB, overloading occurs along any border between two hierarchical levels. By coloring in elements with too many supported B-splines we obtain a visualization of this phenomenon, as seen in Figure 10 on a hierarchical mesh with three levels of refinement.

In order to reduce, or completely remove such overloaded regions, we may for LRB extend meshlines from the fine hierarchical level to the coarse level, in order to split the culprit B-splines. The length needed for this extended meshline depends on the polynomial degree of the B-spline to be split. In Figure 11 we see the effects of three types of meshline extension to the LRB-mesh from Figure 10 for a space of bi-cubic splines.

In order to capture what is happening, we take a closer look at overloading in a convex corner in Figure 12 where we show how B-splines from the coarse level of a hierarchical mesh may overlap with B-splines from the fine level in such a way that too many B-splines are active over a given element.

5 Numerical experiments

In order to compare the methods addressed in this paper, we assemble the mass and stiffness matrices associated to discretizations of partial differential equations using IgA or FEM. By computing the condition number of these matrices, we obtain a metric useful for comparison. These matrices arise amongst others in the discretizations of the Poisson equation, and in the computation of the $L^{2}$-projection of a function.

5.1 $L^{2}$-projection

Given a domain $\Omega$, a function $f\colon\Omega\to\mathbb{R}$ in some space of functions $V$, and a finite-dimensional subspace $V_{h}$ of $V$, we are interested in finding the function $u\in V_{h}$ that minimizes the $L^{2}(\Omega)$-error

[TABLE]

This can be reformulated as a variational equation by requiring $s$ to satisfy

[TABLE]

for all $v\in V_{h}$. By introducing a basis $\mathinner{\left\{\varphi_{1},\ldots,\varphi_{N}\right\}}$ for $V_{h}$, which in our case will be one of the THB or LRB-bases, we may write this as a linear equation

[TABLE]

where $\vec{M}$ is the mass matrix and $\vec{c}$ is the vector of coefficients representing $u$ in our chosen basis. The entries for $\vec{M}$ and the right-hand side $\vec{b}$ are given as

[TABLE]

5.2 The Poisson equation

A commonly encountered differential equation is the Poisson equation. Given a function $f\colon\Omega\to\mathbb{R}$, we wish to find a function $u$ in a space of admissible functions $V$ such that

[TABLE]

subject to the boundary conditions

[TABLE]

Here $\Gamma_{D}$ denotes the Dirichlet-boundary and $\Gamma_{N}$ the Neumann-boundary. We assume $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$ and $\Gamma_{D}\cap\Gamma_{N}=\emptyset$. Furthermore, $\vec{n}$ is the outward facing boundary normal to $\Omega$ and $g$ is the prescribed flux along the boundary.

As for the $L^{2}$-projection, the Poisson equation can be reformulated as a variational equation by multiplying the equation by a test-function, and integrating over the domain. The requirements on the smoothness of the sought solution $u$ can be relaxed, by moving some derivatives onto the test-functions. Again, we seek the solution $u$ in a subspace $V_{h}$ of $V$ spanned by a set of basis functions $\mathinner{\left\{\varphi_{1},\ldots,\varphi_{N}\right\}}$. The variational form of the Poisson-equation reads

[TABLE]

for all $v\in V$. Rewriting this in terms of the basis functions, we obtain the linear equation

[TABLE]

where $\vec{A}$ is the stiffness matrix of the problem. The entries of $\vec{A}$ and $\vec{b}$ are given as

[TABLE]

5.3 Condition Numbers

The condition number of a matrix $\vec{B}\in\mathbb{R}^{n\times n}$ quantifies how sensitive the solution $x$ to the linear system $\vec{B}\vec{x}=\vec{y}$ is to small perturbations both in $\vec{B}$ and the right-hand side $\vec{y}$ and is formally defined as

[TABLE]

where $\|\cdot\|$ is some matrix norm. Note that the condition number is norm-dependent, but all matrix norms are equivalent on $\mathbb{R}^{n\times n}$. We will be computing the condition numbers in the 2-norm, and in this specific setting for normal matrices the condition number can be computed as the ratio between the largest and smallest eigenvalue

[TABLE]

Here $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$, i.e., ordered in a decreasing fashion.

As in johannessenSimilaritiesDifferencesClassical2015a , we chose to estimate the condition numbers of the matrices before imposing any boundary conditions, as imposing boundary conditions can have a large impact on the conditioning of the matrix. The mass matrix $\vec{M}$ is non-singular, even with no imposed boundary condition. The stiffness matrix $\vec{A}$ however will be singular, and have a zero-eigenvalue of multiplicity one.

In addition to this, the computation of the smallest eigenvalue of a matrix is a numerically unstable procedure. We will therefore estimate the condition numbers as follows:

[TABLE]

using the second-smallest eigenvalue for the stiffness matrix.

5.4 Numerical results

Below we present results of the numerical simulations. As LRBNO generates higher dimensional spline spaces than THB and LRB we plot the condition numbers as a function of the degrees of freedom. Just plotting the condition numbers as a function of the levels provides less information. By including the dimension of the spline space we obtain a clearer distinction between the methods.

Central refinement

We assemble the stiffness- and mass-matrices on a sequence of meshes corresponding to central refinement, shown at the third refinement in Figures 10 and 11. The results are shown in Figures 13 and 14. Start by noting that for the mass matrix, THB performs better than LRB with no modifications, while for the stiffness matrices, the two methods are comparable with LRB having a slight advantage. The number of degrees of freedom are the same. By locally modifying the mesh, as is the case for LRBNO and T-LRBNO, we see that the number of degrees of freedom goes up, as expected. The condition number per degree of freedom is smallest for T-LRBNO.

Diagonal refinements

We now consider the case of diagonal refinement. Again, we use the same hierarchical mesh for LRB and THB. We will only consider one sequence of meshes with local modifications. In the diagonal refinement setting, the corners of the refined region are sufficiently close to each other so that we need to make a decision on which direction to refine in. The diagonally refined mesh is not compatible with a T-spline type mesh, and will therefore not be taken into consideration here.

We assembly stiffness- and mass-matrices on the meshes displayed in Figure 15. The results are shown in Figures 16 and 17. Note that for the diagonal refinement, the number of degrees of freedom generated when removing overloading, shown in the mesh in Figure 15c, is a fair bit larger than the unmodified counterparts. Despite this, LRBNO outperforms THB and LRB by a significant amount when it comes to the mass matrix. The conditioning of the stiffness matrix on the other hand grows approximately linearly with the number of degrees of freedom, and no significant effect of the overload-reduction can be seen.

6 Conclusion

We have addressed differences and similarities of Truncated Hierarchical B-splines (THB) and Locally Refined B-splines (LRB) on similar hierarchical meshes. The overall conclusion is that there are no big differences between the methods with respect to condition numbers of mass and stiffness matrices for the example meshes addressed.

• •

When THB and LRB are run on identical meshes THB has better conditions numbers for the mass matrix except for most complex example run, the diagonal example in Figures 16 and 15. The behaviour of the stiffness matrix is very similar for both methods.

• •

When making a mesh for LRB that has no overloading the condition numbers for the mass matrix of LRB are smaller than those of THB, with condition numbers of stiffness being similar. It should be noted that using meshes for LRB that has no overloading guarantees that the B-splines generated are linearly independent, and that the number of B-splines covering an element is the minimal needed for spanning the polynomial space over the element. For hierarchical meshes of bi-degree less than $(4,4)$ there is always linear independence in the set of LR B-splines generated,. For bi-degree $(4,4)$ and higher linear dependence can occur in very special configurations when the elements outside two opposing concave corners of a refinement region is covered by the same B-spline from a cruder level. This happens for bi-degree (4,4) when a refinement region is split if just one element from the cruder level is not refined, e.g., the refinement region is locally very narrow.

When trying to represent hierarchical refinements using T-splines as in Figure 9 there is a region of one directional refinement of length two just outside the boundary of the refinement region. This gives a smoother transition between refinement levels that can also be replicated by LRB. The results in Figures 13 and 14 show a better behaviour than going directly from one refinement level to the next. Having such an intermediate level of refinement if possible is advantageous. However, in situations such as the diagonal refinement in Figure 15 this is not possible.

Most often THB is described as based on dyadic sequences of grids determined by scaled lattices over which uniform B-spline spaces are defined. This implies that there is single knot multiplicity along domain boundaries. However, variants of THB are published GIANNELLI2016337 where open knots are used along the domain boundary. In Section 2 we have shown that open knot vectors are preferable, not only with respect to simplified interpolation of boundary conditions, but also to avoid that the condition number of the mass matrix is biased by the boundary B-splines. As we see the same effect for LR B-splines we have a strong recommendation that open knot vectors are used for locally refined splines, rather than single knot multiplicity at domain boundaries.

7 Acknowledgement

This project has received funding from the The Research Council of Norway under grant agreement No 270922.