Adaptive refinement with locally linearly independent LR B-splines: Theory and applications

In this paper we describe an adaptive refinement strategy for LR B-splines. The presented strategy ensures, at each iteration, local linear independence of the obtained set of LR B-splines. This property is then exploited in two applications: the construction of efficient quasi-interpolation schemes and the numerical solution of elliptic problems using the isogeometric Galerkin method.


Introduction
Since the '70s, curves and surfaces in engineering are usually expressed by means of computer aided design (CAD) technologies, such as B-splines and non-uniform rational Bsplines (NURBS).Thanks to properties like nonnegativity, local support and partition of unity, they allow for an easy control and modification of the geometries they describe, and this motivates their undisputed success as main modeling tools for objects with complex shapes in engineering; see, e.g., [21,7,18] and references therein.On the other hand, Bsplines also provide a very efficient representation of smooth piecewise polynomial spaces, and so are a popular ingredient in the construction of approximation schemes; see, e.g., [2,22,17] and references therein.
More recently, the advent of isogeometric analysis (IgA) has integrated spline and CAD technologies into finite element analysis (FEA); see, e.g., [8,1].IgA aims to unify the geometric description of the domain of the differential problem with its numerical resolution, in order to expedite the simulation process and gaining in accuracy.In addition to the properties listed above, B-splines and NURBS feature other qualities appreciated in this context, such as (local) linear independence and high global smoothness.
Tensor structures admit a simple but powerful multivariate extension of univariate splines and B-splines.On the other hand, they lack adequate local refinement.The constantly increasing demand for higher precision in simulations and reverse engineering processes requires the possibility to apply adaptive local refinement strategies, in order to reduce the approximation error while keeping the computational costs low.This request for adaptivity, triggered the interest in new formulations of B-splines and NURBS, still based on local tensor structures [12,23,9,5,13,10,11].All these new classes of functions are defined on locally refined meshes, in which T-vertices in the interior of the domain are allowed, the so-called T-meshes.
Locally refined B-splines, or in short LR B-splines [10], are one of these new formulations, and their definition is inspired by the knot insertion refinement process of tensor B-splines.These latter are defined on global knot vectors, one per direction.The insertion of a new knot in a knot vector corresponds to a line segment in the mesh crossing the entire domain.This refines all the B-splines whose supports are crossed.Instead, LR B-splines are defined on local knot vectors and the insertion of a new knot is always performed with respect to a particular LR B-spline.As a consequence, the LR B-spline definition is consistent with the tensor B-spline definition when the underlying mesh at the end of the process is a tensor mesh, and the formulation of LR B-splines remains broadly similar to classical tensor B-splines even though they allow for local refinements.This makes them one of the most elegant extensions of univariate B-splines on local tensor structures.
LR B-splines possess almost all the properties of classical tensor B-splines.Unfortunately, they are not always linearly independent.To this day, it is not yet known what are the precise conditions on the locally refined mesh to ensure a linearly independent set of LR B-splines, but some progress has been made in this direction.In [10] an efficient algorithm to seek and destroy linear dependence relations has been introduced, but it does not handle every possible locally refined mesh.In [20] a first analysis on the necessary conditions for encountering a linear dependence relation has been presented.There, it has also been proved that, for any bidegree, a linear dependence relation in the LR B-spline set involves at least 8 functions.In [3,4] a characterization of the local linear independence of LR B-splines has been provided.Such a strong property is guaranteed only on locally refined meshes with strong constraints on the lengths and positions of the line segments that yield particular arrangements of the LR B-spline supports.On the other hand, a practical adaptive refinement strategy to produce meshes with the local linear independence property is still missing in the literature.To the best of our knowledge, the only mesh construction that leads to a locally linearly independent set of LR B-splines can be found in [4].Such a construction, however, cannot be considered as a practical strategy because the regions to be refined and the maximal resolution, i.e., the sizes of the smallest cells in the domain induced by the mesh, must be chosen a priori.
In this paper, we describe a practical refinement strategy ensuring the local linear independence of the corresponding LR B-splines.Such a property is appealing as it admits, e.g., the construction of efficient quasi-interpolation schemes and the unisolvency of linear systems obtained by isogeometric discretization of differential problems based on such LR-splines.The remainder of the paper is divided into 5 sections.Section 2 contains the definition of LR B-splines and a summary of their main properties, whereas Section 3 describes the mesh refinement strategy and is the core of the paper.Sections 4 and 5 present applications of the refinement strategy in the context of quasi-interpolation and isogeometric Galerkin discretizations of elliptic problems.We end in Section 6 with some concluding remarks.
Throughout the paper, we assume the reader to be familiar with the definition and main properties of (univariate) B-splines, in particular with the knot insertion procedure.An introduction to this topic can be found, e.g., in the review papers [18,17] or in the classical books [2] and [22].

Locally refined B-splines
In this section, we introduce locally refined B-splines, or in short LR B-splines, and discuss several of their properties, following the terminology from [20].We denote by Π p the space of univariate polynomials of degree less than or equal to p, and by Π p p p the space of bivariate polynomials of degrees less than or equal to p p p = (p 1 , p 2 ) component-wise.Furthermore, we denote by B[x x x, y y y] the bivariate B-spline defined on the (local) knot vectors x x x = (x 1 , . . ., x p 1 +2 ) and y y y = (y 1 , . . ., y p 2 +2 ), where x i ≤ x i+1 and y i ≤ y i+1 for all i.The bidegree of B[x x x, y y y] is p p p = (p 1 , p 2 ) and is implicitly specified by the length of x x x and y y y.
In order to define LR B-splines, we first introduce the concept of box-partition.
Definition 2.1.Given an axis-aligned rectangle Given a box-partition E, we define the vertices of E as the vertices of its elements.A meshline is an axis-aligned segment contained in an edge of an element of E, connecting two and only two vertices of E located at its end-points.The collection of all the meshlines of the box-partition is called mesh, and denoted by M. A meshline can be expressed as the Cartesian product of a point in R and a finite interval.Let α ∈ R be the value of such a point and let k ∈ {1, 2} be its position in the Cartesian product.If k = 1 the meshline is vertical and if k = 2 the meshline is horizontal.We sometimes write k-meshline to specify the direction of the meshline, and (k, α)-meshline to specify exactly on which axis-parallel line in R 2 the meshline lies.A vertex of E is called T-vertex if it is the intersection of two collinear meshlines and another meshline, say γ, orthogonal to them.We call the T-vertex vertical if γ is vertical, and horizontal otherwise.
For defining splines of a certain bidegree p p p = (p 1 , p 2 ) and smoothness across the meshlines, we also need the notion of multiplicity of a meshline.This is a positive integer associated with every meshline in M. For a k-meshline this number is assumed to be maximally p k + 1.A meshline in M has full multiplicity if its multiplicity is maximal, and we say that M is open if every boundary meshline has full multiplicity.If all the meshlines of the box-partition have the same multiplicity m we say that M has multiplicity m.When the T-vertices of E occur only on ∂Ω and all collinear meshlines have the same multiplicity, the corresponding mesh is called tensor mesh.Figure 1 shows an example of a box-partition and its associated mesh.
Furthermore, we say that B[x x x, y y y] has minimal support on M if • it has support on M, • the multiplicities of the interior meshlines in M(x x x, y y y) are equal to the multiplicities of the corresponding meshlines in M, and • there is no collection γ of collinear meshlines in M\M(x x x, y y y) such that supp B[x x x, y y y]\γ is not connected.
Figure 2 shows examples of B-splines of bidegree (2, 2) with support on a mesh of multiplicity 1.In particular, the B-splines in (b)-(c) have minimal support, whereas the support of the B-spline in (d) can be disconnected by the collection of meshlines γ, visualized by dashed lines in the figure .Given a mesh M and a B-spline B[x x x, y y y] with support in M, assume that it does not have minimal support on M.Then, there exists a collection of (k, α)-meshlines γ such that supp B[x x x, y y y]\γ is not connected and either γ is in M\M(x x x, y y y) or γ ⊆ M(x x x, y y y), i.e., α is an internal knot of x x x for k = 1 or y y y for k = 2, but its meshlines have lower multiplicities in M(x x x, y y y) than in M. Assume that the meshlines in γ have all the same multiplicity m in M. Denoting by µ(α) ≥ 0 the number of times α appears in the knot vector of B[x x x, y y y] in the k-th direction, then m − µ(α) is strictly positive as B[x x x, y y y] has support, but not minimal support, on M. One could consider such α as an extra knot, of multiplicity m − µ(α), with respect to the knot vector of B[x x x, y y y] in the k-th direction (in x x x if k = 1 and in y y y if k = 2), and perform knot insertion on B[x x x, y y y].If α was already a knot of B[x x x, y y y], so µ(α) ≥ 1, this means rising its multiplicity by m − µ(α).The resulting generated B-splines will still have support on M and eventually they will also have minimal support on M. As an example, the collection γ highlighted with dashed lines in Figure 2(d) is made of (2, α)-meshlines, for some α, of multiplicity 1.Such α can be inserted as new knot of multiplicity 1 in the knot vector in the y-direction of the considered B-spline to refine it in two B-splines via knot insertion.
The LR B-splines are generated by means of the above procedure.We start by considering a coarse tensor mesh and we refine it by inserting collections of collinear meshlines, one at a time, of the same multiplicity.On the initial mesh we consider the standard tensor B-splines and whenever a B-spline in our collection has no longer minimal support during the mesh refinement process, we refine it by using the knot insertion procedure.The LR B-splines will be the final set of B-splines produced by this algorithm.In the following definitions we formalize this by describing the mesh refinement process in our framework.Definition 2.3.Given a box-partition E and an axis-aligned segment γ, we say that γ traverses β ∈ E if γ ⊆ β and the interior of β is divided into two parts by γ, i.e., β\γ is not connected.A split is a finite union of contiguous and collinear axis-aligned segments γ = ∪ i γ i such that every γ i either is a meshline of the box-partition or traverses some β ∈ E. A mesh M has constant splits if each split in it is made of meshlines of the same multiplicity.
Like for meshlines, we sometimes write k-split with k ∈ {1, 2} to specify the direction of the split or (k, α)-split to specify on what axis-parallel line in R 2 the split lies.
When a split γ is inserted in a box-partition E, any traversed β ∈ E is replaced by the two subrectangles β 1 , β 2 given by the closures of the connected components of β\γ.The resulting new box-partition will be denoted by E + γ and its corresponding mesh by M + γ.We also assume that a positive integer µ γ has been assigned to any split γ.The multiplicities of the meshlines in M ∩ (M + γ) and not contained in γ are unchanged.Contrarily, the multiplicities of the meshlines in γ that were already in M are rised by µ γ , and the new meshlines in γ have multiplicity equal to µ γ on M + γ.
The LR B-splines are defined on a class of meshes with constant splits, called LR-meshes.Thus, from now on, we restrict our attention to meshes that have constant splits.In particular, we note that when refining a mesh M by inserting a split γ, either γ is made solely of (a) (b) Figure 3: Two meshes.Assume that the boundary has a multiplicity large enough so that it is possible to define a B-spline of bidegree p p p on it.Then, the mesh in (a) is not an LR-mesh because it cannot be built as a sequence of split insertions.The mesh in (b) is an LR-mesh similar to the one in (a).
new meshlines or it is made solely of meshlines already on M, in order for M + γ to have constant splits.Definition 2.4.Given a mesh M with constant splits, a B-spline B[x x x, y y y] with support on M and a split γ, we say that γ traverses B[x x x, y y y] if the interior of supp B[x x x, y y y] is divided into two parts by γ, i.e., supp B[x x x, y y y]\γ is not connected and either γ is in M\M(x x x, y y y) or γ ⊆ M(x x x, y y y) but the multiplicity of its meshlines is lower in M(x x x, y y y) than in M.
We are now ready to define the mesh refinement process and the LR B-splines.The meshes generated by this procedure will be called LR-meshes.Definition 2.5.Given a bidegree p p p = (p 1 , p 2 ), let M 1 be a tensor mesh such that the set of standard tensor B-splines of bidegree p p p on M 1 is nonempty, and denote it by B 1 .We then define a sequence of meshes M 2 , M 3 , . . .and corresponding function sets B 2 , B 3 , . . .as follows.For i = 1, 2, . .., let γ i be a split such that M i+1 := M i + γ i has constant splits and such that the support of at least one B-spline in B i is traversed by a split in M i+1 .On this refined mesh M i+1 , the new set of B-splines B i+1 is constructed by the following algorithm.

Initialize the set by
As long as there exists B[x x x j , y y y j ] ∈ B i+1 with no minimal support on M i+1 : (a) Apply knot insertion: The mesh generated at each step is called LR-mesh and the corresponding function set is called LR B-spline set.
Obviously not every mesh is an LR-mesh.For instance, one could consider meshes that do not have constant splits or meshes that cannot be built through a sequence of split insertions as the mesh depicted in Figure 3(a).In general, the mesh refinement process producing a given LR-mesh M = M N is not unique.Indeed, the split insertion ordering can often be changed.Nevertheless, the LR B-spline set on M is well defined because it is independent of such insertion ordering, as proved in [10,Theorem 3.4].
Given an LR-mesh, the corresponding LR B-splines have several desirable properties for applications.By their definition, it is clear that • they are nonnegative, • they have minimal support, and • they can be expressed by the LR B-splines on finer LR-meshes using nonnegative coefficients (provided by the knot insertion procedure).
Furthermore, it is possible to scale them by means of positive weights so that they also form a partition of unity; see [10,Section 7].Unfortunately, they are not always linearly independent.Figure 4 shows an example of linear dependence among the LR B-splines of bidegree (2, 2) defined on an LR-mesh of multiplicity 1.To this day, it is not yet known what are the precise conditions on the LR-mesh to ensure a linearly independent set of LR B-splines.
In [4] a characterization of the local linear independence of LR B-splines has been provided.Such a strong property is guaranteed only on LR-meshes with strong constraints on the split lengths and positions that yield particular arrangements of the LR B-spline supports.This last statement is formalized in the following.Definition 2.6.Given a mesh M, let B[x x x 1 , y y y 1 ] and B[x x x 2 , y y y 2 ] be two different LR B-splines defined on M. We say that B[x x x 2 , y y y 2 ] is nested in B[x x x 1 , y y y 1 ], and we write y y y 2 ] has a higher (or equal) multiplicity when considered in M(x x x 1 , y y y 1 ) than in M(x x x 2 , y y y 2 ).
An open mesh where no LR B-spline is nested is said to have the non-nested support property, or in short the N 2 S-property.
The definition of nested LR B-splines was formulated for the first time in [3].Definition 2.6 is different but equivalent to it (see Appendix A). Figure 5 shows an example of an LR Bspline nested into another.The following result, presented in [4], relates the local linear independence of the LR B-splines to the N 2 S-property of the mesh.4. The LR B-splines form a partition of unity, without the use of scaling weights.
An element of E for which item 3 of Theorem 2.7 holds is said to be non-overloaded.Note that (p 1 + 1)(p 2 + 1) is the dimension of the polynomial space over the element.
In [4] one can also find an algorithm to construct LR-meshes so that the N 2 S-property is fulfilled.This approach, however, has a relevant drawback for practical purposes: the regions to be refined and the maximal resolution have to be chosen a priori.Moreover, the algorithm cannot be stopped prematurely, before having inserted all the splits determined initially.In practice, one rarely knows in advance where the error will be large and how fine the mesh has to be chosen to reduce it under a certain tolerance.
In the next section, we present an alternative way to generate LR-meshes so that the N 2 S-property is guaranteed.

N 2 S-structured mesh refinement strategy
In this section, we define a local refinement strategy that ensures the N 2 S-property for the obtained meshes.It consists of two steps.First, we apply the so-called structured mesh refinement, defined in [15], to the LR B-splines whose contribution to the approximation error is larger than a given tolerance.Then, we slightly modify the obtained mesh by prolonging As opposed to the classical finite element method, in which the refinement is applied to the box-partition elements, the structured mesh refinement is a refinement applied to the function space, i.e., we select for refinement the LR B-splines contributing more to the approximation error rather than the box-partition elements where a larger error occurs.This approach is justified by the fact that on an LR-mesh, any new split inserted must traverse at least the support of one LR B-spline.If we choose to select the elements where the error is larger, then the refinement has to be extended anyway to traverse the support of at least one LR B-spline containing the elements, resulting in a refinement of the LR B-spline basis.Moreover, since several LR B-splines contain such elements, those chosen for the refinement extension could be not those contributing more to the error, resulting in a suboptimal refinement, or we could refine more LR B-splines than necessary, wasting degrees of freedom.
Once the LR B-splines are selected, we refine them by halving the interval steps in their knot vectors.This corresponds to the insertion of a net of meshlines in the B-spline supports on the mesh.We therefore perform the LR B-spline generation algorithm described in Definition 2.5.Every selected LR B-spline is fragmented into LR B-splines with smaller support and replaced by them.The LR-mesh obtained in this way will be called a structured LR-mesh.
In summary, the structured mesh refinement consists of two steps: 1. LR B-splines are selected to be refined and not box-partition elements; The LR B-splines with support in the region highlighted in (c) are linearly dependent.In particular, the region corresponds to the support of an LR B-spline that has many nested LR B-splines in it.One can prove the existence of the linear dependence relation by computing the spline space dimension and the number of LR B-splines defined on the mesh as explained in the examples of [20].This configuration can be reproduced for any bidegree (p 1 , p 2 ) with p k ≥ 4 for k = 1, 2.
2. the interval steps of their knot vectors are halved to obtain the new LR-mesh.
Figure 6 shows two iterations of such refinement.In general, the structured mesh refinement does not generate LR-meshes with the N 2 S-property.The LR-mesh in Figure 6(f) is an example as explained in Figure 7. Furthermore, the structured mesh refinement may produce linearly dependent sets of LR B-splines.Figure 8 shows an example for bidegree (4,4).On the other hand, the standard B-splines defined on a plain tensor mesh are locally linearly independent, and the meshes generated by the structured mesh refinement are locally tensor meshes far from the boundary of the region where the structured mesh refinement is applied.The LR B-splines defined in these zones of the mesh behave like the standard Bsplines, and therefore are locally linearly independent.On the boundary of the region where the refinement has been applied, LR B-splines with smaller support can be nested in LR B-splines with larger support.Hence, in such case the resulting LR-mesh does not have the N 2 S-property.
The idea for our refinement strategy, which will be called N 2 S-structured mesh refinement, is therefore to recover the N 2 S-property in the mesh by slightly modifying it in these transition regions.When an LR B-spline B[x x x 2 , y y y 2 ] is nested into another LR B-spline B[x x x 1 , y y y 1 ], one could prolong the splits in M(x x x 2 , y y y 2 ) in some direction to traverse entirely supp B[x x x 1 , y y y 1 ].This, by Definition 2.5, would refine B[x x x 1 , y y y 1 ] in LR B-splines that turn out not to have nested LR B-splines in their supports anymore.This last statement is formalized in Corollary 3.3.To this end, we first prove the N 2 S-property for LR-meshes with a particular structure.
Definition 3.1.An LR-mesh M on the domain Ω is said to be tensorized in the k-th direction, for k ∈ {1, 2}, if all the internal k-meshlines in M are contained in k-splits crossing Ω entirely, i.e., there are no vertical, if k = 1, or horizontal, if k = 2, T-vertices in the interior of Ω. Proposition 3.2.Let M be an LR-mesh tensorized in the k-th direction for some k ∈ {1, 2}.Then, the LR B-splines defined on M are all non-nested.
Proof.Without loss of generality, we can assume that M is tensorized in the first direction, i.e., the vertical meshlines are all contained in vertical splits crossing the domain entirely.This means that in M no vertical meshline ends in the interior of the domain and therefore in the interior of the support of any LR B-splines defined on M. We now proceed by contradiction and assume that there exists an LR B-spline in M, say B 2 = B[x x x 2 , y y y 2 ], nested in another, say B 1 = B[x x x 1 , y y y 1 ]; see Definition 2.6.Because of the tensorization in the first direction, this can only happen if they share the same knot vector in the x-direction, x x x 1 = x x x 2 .In particular, their supports have the same extreme values in the x-direction.This implies that all the horizontal splits, counting the multiplicities, of M traversing supp B 2 must traverse supp B 1 as well.Since B 2 is nested in B 1 and B 1 has minimal support (it is an LR B-spline), it follows that y y y 1 = y y y 2 , and as a consequence we have B 1 = B 2 .This is a contradiction and concludes the proof.Corollary 3.3.Given an LR-mesh M, let B = B[x x x, y y y] and B 1 = B[x x x 1 , y y y 1 ], . . ., B n = B[x x x n , y y y n ] be LR B-splines defined on M such that B 1 , . . ., B n B. Let N be the mesh defined by the restriction of M to the meshlines of M(x x x, y y y), M(x x x 1 , y y y 1 ), . . ., M(x x x n , y y y n ).Then,

there are at least one horizontal T-vertex and one vertical T-vertex of N in the interior
of supp B; 2. by extending all the splits of N in some direction to cross supp B entirely, B is refined, by Definition 2.5, in LR B-splines that do not have any nested LR B-splines anymore. Proof.
1. Assume that there are no vertical T-vertices of N in the interior of supp B. Then, N would be tensorized in the first direction.By Proposition 3.2, it would imply that all the LR B-splines defined on N are non-nested, which is a contradiction.Analogously, one can prove that at least one horizontal T-vertex of N must be in the interior of supp B.  2. Since B contains the support of other B-splines, N = M(x x x, y y y) and in particular there exist at least one horizontal T-vertex and one vertical T-vertex by the previous item.We now focus on the vertical T-vertices, but of course the same argument can also be carried out for the horizontal T-vertices.We extend all the vertical splits in N to cross supp B entirely, and denote this new mesh as N .By Definition 2.5, the extensions trigger a refinement of B via knot insertions.N is tensorized in the first direction and, by Proposition 3.2, no LR B-spline defined on N is nested into another.
The extension of the splits considered in item 2 of Corollary 3.3 will be called a onedirectional tensor expansion of B 1 , . . ., B n in B. An example is illustrated in Figure 9.
The N 2 S-structured mesh refinement is defined algorithmically as follows.We start from a structured mesh refinement to obtain a new set of LR B-splines.We then collect in a set B all those LR B-splines that have nested LR B-splines in their supports.If B is non-empty, we select an LR B-spline B in B and we apply a one-directional tensor expansion to it.This triggers a refinement of the LR B-spline set, and therefore it changes also the set B. We repeat this procedure till B becomes empty.In Theorem 3.4 we shall prove that this always happens in a finite number of steps.This procedure is sketched in Algorithm 1.The one-directional tensor expansions are performed by alternating the direction for i even and odd, respectively, in order to bound the thinning of the box-partition elements in a specific direction and preserve the uniformity of the mesh as much as possible.The LR-mesh obtained in this way will be called an N 2 S-structured LR-mesh, or in short N 2 S 2 LR-mesh.Theorem 3.4.Given an axis-aligned rectangular domain Ω ⊆ R 2 , let B 1 be the set of standard bivariate B-splines defined on the open tensor mesh whose meshlines are the edges of ∂Ω.Then, 1. the LR B-spline sets B i provided by Algorithm 1 are well defined, i.e., the set B of the algorithm becomes empty in a finite number of iterations, for every index i ≥ 2; 2. all the LR B-splines in B i are non-nested, for every i ≥ 1.
Proof.Without loss of generality, we can assume that Ω = [0, 1] × [0, 1].We proceed by induction on the index of the B-spline set.For i = 1, B 1 is the set of standard B-splines on the open tensor mesh equal to the domain's boundary and we know they are locally linearly independent.By Theorem 2.7 this is equivalent to be all non-nested.Assume now that B i for some i ≥ 1 is well defined and that the functions in it are all non-nested.Let us then prove that also B i+1 is well defined and there is no LR B-spline nested into another LR B-spline of it.At every loop iteration in the algorithm, the LR B-splines that have a nested LR B-spline in their support are collected in the set B. Therefore, whenever we can show that B becomes empty after a certain iteration of the loop, we can immediately conclude both statements in the theorem.
By Corollary 3.3, all the one-directional tensor expansions performed to define the set B i+1 can be done in the same direction k ∈ {1, 2}, which is therefore fixed once and for all by the index i + 1.The length of the LR B-spline supports in the (3 − k)-th direction at any iteration of the loop cannot become shorter than 2 −(i+1) regardless of the number of one-directional tensor expansions applied until then.This is because the (3 − k)-splits on the LR-meshes defined in the loop are fixed by the structured mesh refinement performed on B i at the beginning of the process and the minimal length of the box-partition elements in the (3 − k)-th direction is 2 −(i+1) .Therefore, the split extensions applied when performing a one-directional tensor expansion in the k-th direction have lengths bounded from below by 2 −(i+1) in all the steps of the loop.This means that in a finite number of one-directional tensor expansions a k-split could be extended up to the domain's boundary, if needed, to remove nestedness issues, as these extensions cannot become arbitrarily small.In the worst case scenario, we must extend all the k-splits to cross entirely the domain.However, in this case, the resulting LR-mesh would be tensorized in the k-th direction.By Proposition 3.2, there are only non-nested LR B-splines on this LR-mesh and thus B becomes empty in a finite number of loop iterations.
In practice, the loop related to B stops quickly and the N 2 S 2 LR-meshes are far from being entirely tensorized in one direction.In Figure 10 we depict (a) the structured LRmesh, (b) the corresponding N 2 S 2 LR-mesh, and (c) the LR-mesh proposed in [4], obtained by performing 7 iterations of diagonal refinement in [0, 1] 2 , using bidegree (2, 2).For ease of comparison, we also indicate the number of LR B-splines defined on each of these meshes.We recall that the LR B-splines are not locally linearly independent on the structured LR-mesh, Algorithm 1: N 2 S-structured mesh refinement.whereas they are on the N 2 S 2 LR-mesh and the LR-mesh proposed in [4].
In Figure 10(b) and Figure 11 one can see how the refinement in the N 2 S 2 LR-meshes propagates from the region where the structured mesh refinement has been applied.In all the considered cases, the refinement does not heavily spread out.It is important to highlight, however, that the prolongation of the splits needed to recover the N 2 S-property is not unique.Indeed, when refining an LR B-spline to remove nestedness issues, the inserted split prolongations refine not only the considered LR B-spline but in general also other LR B-splines in the neighborhood.Then, some of the newly introduced neighboring LR Bsplines might not need a one-directional tensor expansion anymore and the ordering used for removing nestedness has thus an effect on the resulting mesh.
One might consider to treat all the LR B-splines "in parallel", i.e., first collect all the split extensions needed to remove nestedness in all the LR B-splines requiring a treatment and then insert all of them at the same time to refine the function basis.This could result in a more uniform propagation of the refinement out of the region where the structured mesh has been applied.On the other hand, by doing this, some split extensions could be unnecessary for recovering the N 2 S-property.Therefore, in general, also the number of LR B-splines on these meshes would be higher than the number obtained when treating one LR B-spline at a time.In the examples presented in this paper we do not remove nestedness "in parallel".Hence, the resulting N 2 S 2 LR-meshes depend on the order used when the onedirectional tensor expansions are applied.On the other hand, the number of LR B-splines will be closer to the number of LR B-splines obtained when performing only the structured mesh refinement, i.e., closer to the "optimal" number of LR B-splines needed to reduce the error while preserving the local linear independence.
We finally remark that one might also opt for full tensor expansions in the supports, instead of one-directional tensor expansions, to solve nestedness issues.The proof of Theorem 3.4 could be easily rephrased for the case of full tensor expansions.The key is that we only prolong splits provided by the structured mesh refinement performed at the beginning of the process.Therefore, if we do full tensor expansions, in the worst case scenario we would end up with a standard tensor mesh of size h = 2 −(i+1) to define the set B i+1 , instead of an LR-mesh tensorized in one direction.In such case, B would still become empty in a finite number of loop iterations.However, we decided to do the expansion of the splits only in one direction at a time because it results in the least propagation.

Application I: Quasi-interpolation
A quasi-interpolation method is a procedure to compute the coefficients assigned to the basis elements of a prescribed function space, with the aim of approximating a given arbitrary function or data set.The resulting approximant is called a quasi-interpolant (QI).The computation of any of such coefficients may depend only on the data/function restricted to the corresponding basis element's support (local method), and perhaps some neighboring other basis elements' supports, or it can depend on the data/function in the entire domain (global method), as in the least-squares method.Given a function f and an approximation space, whose basis is denoted by B, we write a related QI in the form where λ B (f ) is the coefficient of the basis element B ∈ B computed by the selected method.Definition 4.1.A quasi-interpolation method such that Qf = f for all f in a space V is said to reproduce the space V .
When using spline spaces of bidegree p p p as approximation spaces, a common requirement is that the polynomial space Π p p p is reproduced by the quasi-interpolation method, in order to ensure good approximation properties.A general recipe for constructing local quasiinterpolation methods for tensor spline spaces, with the polynomial reproduction property, can be found in [16].
Recipe 4.2.Let f be a given function defined on the rectangle Ω.Given a bidegree p p p, let M be an open tensor mesh on Ω, and let B( M) be the set of tensor B-splines of bidegree p p p on M. We compute the coefficient λ B (f ), for every B = B[x x x, y y y] ∈ B( M), as follows.
1. Let U ⊆ R 2 be an open set that intersects the interior of supp B (for instance, U can be a box-partition element of M(x x x, y y y)), and let B(U ) be the subset of B( M) consisting of all the tensor B-splines not identically zero on U . 2. Choose a local polynomial approximation method P U such that P U g = g for all g ∈ Π p p p defined on U (typical choices are least-squares or interpolation methods).Letting g |U be the restriction of g to U , we can write for some coefficients b B (f ) provided by the chosen local approximation method.We have tested the quasi-interpolation strategy described in Recipe 4.3 on N 2 S 2 LRmeshes to approximate polynomials and transcendent functions.Given an N 2 S 2 LR-mesh M, this recipe requires the construction of a QI based on B( M B ) for each of the LR B-splines B ∈ B LR (M) of bidegree p p p = (p 1 , p 2 ).Following Recipe 4.2, we have used interpolation as local approximation method in the computation of these QIs.More precisely, we have selected a unisolvent set of (p 1 + 1)(p 2 + 1) interpolation points, organized in a tensor grid over a single box-partition element of M B , and then we have set a linear system by evaluating f and the tensor B-splines in B( M B ) at these points.This guarantees polynomial reproduction (actually spline reproduction) of the quasi-interpolation method in the spaces B( M B ), for B varying in B LR (M).Therefore, also the resulting quasi-interpolation method on B LR (M) has the polynomial reproduction property thanks to Proposition 4.4.Indeed, in all the tests with polynomial functions of bidegree at most p p p, the maximum error was in the order of the machine precision, regardless of the number of iterations performed to construct the N 2 S 2 LR-mesh.The maximum error was computed on a uniform 150 × 150 grid.
As test with a transcendent function, we have considered which is characterized by three steep peaks on the square [−1, 1] 2 located at (−0.3, −0.3), (0, 0) and (0.3, 0.3); see Figure 12.This function has also been used in [25] to investigate the approximation power of a similar quasi-interpolation method developed for THB-splines.In Table 1, we compare the number of basis functions of bidegree (2, 2) when considering global tensor meshes and local N 2 S 2 LR-meshes for different levels of maximal resolution (for level , the smallest box-partition elements on the mesh have length 2 − ).The N 2 S 2 LR-mesh is obtained by refining the LR B-splines whose supports contain one of the three points where a peak occurs via structured mesh refinement and then by recovering the N 2 S-property via onedirectional tensor expansions.For a given maximal resolution level, the optimal maximum error, i.e., the maximum error when using the global tensor mesh, is preserved by the N 2 S 2 LR-mesh.However, the number of B-splines is significantly different and the discrepancy exponentially grows with the maximal resolution level.

Application II: Isogeometric analysis
Isogeometric analysis (IgA), introduced in [14], is a technique to perform numerical simulations on complex geometries.The numerical solution is represented by means of the same functions used for the domain modeling.Nowadays, complex geometries are expressed in terms of computed aided design (CAD) technologies, such as B-splines, non-uniform rational B-splines (NURBS) and their generalizations to address adaptive refinements.
In this section, we adopt the IgA approach, using our LR refinement strategy, to approximate the solution of the Poisson problem on Ω = [0, 1] 2 , whose exact solution is u(x, y) = arctan 100 (x − 1.25) 2 + (y + 0.25) 2 − π 3 ; see Figure 13.This example is a good benchmark for numerical schemes, as the sharp interior layer of the exact solution highlights the approximation quality, and it has been used extensively in the literature; see, e.g., [19,15,6].
In the context of Galerkin discretizations, two properties are desirable: • (local) linear independence of the space generators, • refinement adaptivity.
The linear independence of the functions used as building blocks of the numerical solution avoids the numerical complexity posed by the singularity of the matrices associated to the problem discretization.The refinement adaptivity is desired for balancing accuracy and computational cost as it allows for a higher precision, only there where it is needed to reproduce fast variations of the exact solution.LR B-splines on N 2 S 2 LR-meshes are suitable candidates as both the (local) linear independence of the space generators and the adaptivity of the refinement are guaranteed.In Figure 14 we compare the L ∞ -norm and the L 2 -norm of the error (Figures 14(a) and 14(b) respectively), using bidegree (2, 2) with global tensor meshes and local N 2 S 2 LRmeshes for different levels of maximal resolution to approximate the solution of the Poisson problem (2).The N 2 S 2 LR-mesh is computed by first applying the structured mesh to the LR B-splines whose supports intersect the curve where the sharp interior layer in the exact solution occurs, and then by performing one-directional tensor expansions to recover the N 2 S-property.The error norms, which are computed discretely on a uniform grid of 1000 × 1000 points, are plotted in log-log scale with respect to the number of LR B-splines on the mesh.The solid line with circular markers shows the decay when using global tensor meshes, whereas the dashed line with star markers shows the decay for the N 2 S 2 LR-meshes.In the figures, the first marker corresponds to the 4 × 4 tensor mesh, for maximal resolution level = 2, and it is the maximal level for which the LR B-spline and standard tensor B-spline sets coincide.When considering a comparable number of functions, the N 2 S 2 LR-mesh leads to a significant reduction of both the L ∞ -norm and the L 2 -norm of the error with respect to the tensor mesh, thanks to the adaptivity of the refinement.

Conclusion
LR B-splines are one of the most elegant extensions of univariate B-splines on local tensor structures that allow local refinement.They possess almost all the properties of classical B-splines, but they are not always linearly independent.Recently, a characterization of LR-meshes ensuring local linear independence of the corresponding LR B-splines has been presented in the literature.However, a practical adaptive refinement strategy for LR-meshes that maintain such a property was missing.In this paper, we have filled this gap by describing an adaptive refinement strategy that generates LR-meshes where the corresponding LR B-splines are locally linearly independent.Subsequently, we have exploited the local linear independence of the LR B-splines to construct efficient quasi-interpolation schemes and to solve elliptic problems using the isogeometric Galerkin method.supp B 2 ⊆ supp B 1 , we have ]x 2  1 , x 2 p 1 +2 [ ⊆ ]x 1 1 , x 1 p 1 +2 [.If z / ∈ x x x 1 , then µ x x x 1 (z) = 0 and therefore µ x x x 2 (z) ≥ µ x x x 1 (z) for any value of µ x x x 2 (z).If z ∈ x x x 1 , then it must be also in x x x 2 , otherwise the (1, z)-split {z} × [y 1  1 , y 1 p 2 +2 ] would traverse B 2 , which would not have minimal support.For the same reason, it must also hold that µ x x x 2 (z) = µ x x x 1 (z).This proves item 1 of Definition A.1.Assume now z / ∈ [x 1 1 , x 1 p 1 +2 ].Since supp B 2 ⊆ supp B 1 , we have x 1 1 ≤ x 2 1 and x 2 p 1 +2 ≤ x 1 p 1 +2 .Therefore, µ x x x 1 (z) = µ x x x 2 (z) = 0.If z ∈ {x 1  1 , x 1 p 1 +2 } but z / ∈ x x x 2 , then trivially µ x x x 2 (z) ≤ µ x x x 1 (z).If z ∈ x x x 2 , then z corresponds to (1, z)-meshlines in ∂supp B 1 ∩ ∂supp B 2 .By assumption, these meshlines have a higher (or equal) multiplicity in M(x x x 1 , y y y 1 ) than in M(x x x 2 , y y y 2 ), which means µ x x x 2 (z) ≤ µ x x x 1 (z).The same line of arguments also applies to the knots in the second direction.This proves item 2 of Definition A.1.

Figure 1 :
Figure 1: Example of a box-partition E of a rectangle Ω in (a), and the mesh corresponding to E in (b).The meshlines are identified by squares showing the associated multiplicities.

Figure 2 :
Figure 2: Support of B-splines of bidegree (2, 2) on a mesh M of multiplicity 1.The mesh is shown in (a).The B-splines whose supports are depicted in (b) and (c) have minimal support on M. The tensor meshes defined by the B-spline's knots are highlighted with thicker lines.On the other hand, the B-spline in (d) does not have minimal support on M: the collection of meshlines contained in the dashed line disconnects its support.

Figure 4 :
Figure 4: Example of linear dependence in the LR B-spline set.The parameterization of an LR-mesh M of multiplicity 1 is considered in (a), and the linear dependence relation among some of the LR B-splines of bidegree (2, 2) defined on M is illustrated in (b).The LR B-splines are represented by means of their supports on the mesh and the tensor meshes induced by their knots are highlighted with thicker meshlines.

Figure 5 : 1 .
Figure 5: Example of nested LR B-splines on the mesh M shown in (a).All the meshlines have multiplicity 1 except those in the left edge of M, highlighted with a double line, which have multiplicity 2. In (b)-(d) three LR B-splines, B[x x x 1 , y y y 1 ], B[x x x 2 , y y y 2 ], B[x x x 3 , y y y 3 ] respectively, of bidegree (2, 2) with minimal support on M are represented by means of their supports and the tensor meshes induced by their knots.All the knots of these LR B-splines have multiplicity 1 except x 31 which has multiplicity 2. Therefore, B[x x x 2 , y y y 2 ] B[x x x 1 , y y y 1 ] but B[x x x 3 , y y y 3 ] B[x x x 2 , y y y 2 ] and B[x x x 3 , y y y 3 ] B[x x x 1 , y y y 1 ], despite that supp B[x x x 3 , y y y 3 ] ⊆ supp B[x x x 2 , y y y 2 ] and supp B[x x x 3 , y y y 3 ] ⊆ supp B[x x x 1 , y y y 1 ], because the shared meshlines in the left edge of supp B[x x x 3 , y y y 3 ], supp B[x x x 2 , y y y 2 ] and supp B[x x x 1 , y y y 1 ] have multiplicity 2 in M(x x x 3 , y y y 3 ) and multiplicity 1 in M(x x x 2 , y y y 2 ) and M(x x x 1 , y y y 1 ).

Figure 6 :
Figure6: Two iterations of the structured mesh refinement of bidegree(2,2).We consider the initial open tensor mesh with internal meshlines of multiplicity 1 in (a). Figure (b) shows the support of an LR B-spline selected for refinement.We refine it by halving the interval steps in its knot vectors.This results in the insertion of a net of meshlines in the LR-mesh as shown in (c).In (d) we select another LR B-spline in the new set of LR B-splines and we refine it as illustrated in (e). Figure (f) depicts the final mesh obtained.

Figure 7 :Figure 8 :
Figure 7: An LR B-spline nested in another LR B-spline, which in turn is nested in another LR B-spline on a structured LR-mesh for bidegree (2, 2).Consider again the mesh in Figure 6(f).In (a)-(c) we depict the supports of three LR B-splines on this mesh.The support in (a) is contained in the interior of the support in (b) and (c), and the support in (b) is contained in the interior of the support in (c).Therefore, the LR B-spline considered in (a) is nested both in the LR B-splines in (b) and (c), and the LR B-spline in (b) is nested in the LR B-spline in (c).Hence, the considered mesh does not have the N 2 S-property.

Figure 9 :
Figure 9: Example of a vertical tensor expansion.We consider five LR B-splines of bidegree (2, 2), namely B and B 1 , . . ., B 4 , with supp B 1 , . . ., supp B 4 contained in the upper left corner of supp B. The mesh N , of multiplicity 1, generated by the meshlines of B, B 1 , . . ., B 4 is depicted in (a) and the supports of the LR B-splines are shown in (b).In (c) we perform a vertical tensor expansion of B 1 , . . ., B 4 in B. In (d) the supports of the new set of LR B-splines are shown: none of them has a nested LR B-spline anymore.

Figure 11 :
Figure 11: Meshes of bidegree (2, 2) obtained by performing 8 iterations of mesh refinement on 3 different regions resulting in a structured LR-mesh (left column) and an N 2 S 2 LR-mesh (right column).

3 .
Since B ∈ B(U ), set λ B (f ) := b B (f ).B-splines containing B. Since, according to Recipe 4.3, any Q B reproduces all polynomials in Π p p p we have g B = λ B (g), ∀B ∈ B LR (M), g ∈ Π p p p , which completes the proof.

Figure 14 :
Figure 14: Decay of the error in (a) the L ∞ -norm and (b) the L 2 -norm, when approximating the solution of problem (2) with B-splines of bidegree (2, 2) on tensor meshes (solid line) and N 2 S 2 LR-meshes (dashed line) for different levels of maximal resolution.
For instance, the (1, x i )-meshlines {x i }×[y jn , y j n+1 ] for n = 1, . . ., s − 1 have all the same multiplicity equal to the multiplicity of x i in x x x.