The analysis and development of adaptive schemes is one of the most active areas of research in the context of isogeometric analysis (IGA), a recent methodology for the solution of partial differential equations with high continuity splines. The main idea of adaptive methods is to obtain a good accuracy of the solution with less computational effort by applying local mesh refinement, hence adaptive IGA schemes require the use of different spline spaces that break the tensor-product structure of B-splines. The most popular alternatives in the IGA research community are T-splines[Sederberg_Zheng], LR-splines [LR-splines, Johannessen2013] and hierarchical splines [Vuong_giannelli_juttler_simeon].
In particular, in this paper we focus on T-splines, introduced by Sederberg et al. in [Sederberg_Zheng] and applied in IGA for the first time in [Bazilervs_Calo_Cottrell_Evans, Dorfel_Juttler_Simeon]
. T-splines are constructed from a T-mesh, a rectangular tiling with hanging nodes, and T-spline blending functions are defined from their local knot vectors, which are computed from the tiling. The mathematical research on T-splines has been very active in recent years, and it has led to the introduction of analysis-suitable (or dual compatible) T-splines[LZSHS12, BBCS12, BBSV2013, Li-Scott], a subset of T-splines with good approximation properties that provide local linear independence.
A standard adaptive scheme based on mesh refinement can be written in a loop of the form [CNX]
and suitable strategies for all the steps of the adaptive scheme are needed in order to guarantee its efficiency. In this paper we focus on the solution of the linear system arising in the SOLVE module for T-splines, and study the optimality of a suitable BPX preconditioner. Several domain decomposition preconditioners have been recently studied in the IGA context: overlapping Schwarz methods [BePa13], balancing domain decomposition by constraint (BDDC) methods [BePa_BDDC, BPSWZ] and dual primal isogeometric tearing and interconnecting methods [KPJ_FETIDP, PS_FETIDP]. Multilevel preconditioners for IGA have been extensively analyzed in the tensor-product setting: the BPX preconditioner in [BHKS13], and multigrid preconditioners [Gahalaut_MG, Hofreither2015_1, Hofreither2015_2, Hofreither2016, Hofreither-Takacs, Manni_MG]. Moreover, preconditioners based on fast solvers for Sylvester-like equations have been proposed in the recent paper [Sangalli_Tani].
The T-splines obtained with the refinement strategy in [Morgenstern_Peterseim] are analysis-suitable by construction, and present a multilevel structure, which makes them very appealing to apply multilevel preconditioners. In the present paper, we first present theoretically optimal multilevel preconditioners for IGA on locally refined T-meshes, extending the results in [BHKS13] to T-splines.
For the study of the optimality of the BPX preconditioner we follow [CNX], writing the preconditioner in the framework of the parallel subspace correction (PSC) method. In this framework, the optimality follows from two basic properties: a stable space decomposition, and the strengthened Cauchy-Schwarz inequality. The proof of these two properties is the core of this work, and as a consequence we obtain that the BPX preconditioner gives a uniformly bounded condition number, which is independent of the mesh size , but depends on the degree as in [BHKS13].
The construction of the BPX preconditioner as in [CNX] is performed by adding a new edge to the T-mesh by bisection, and the new level is defined by the functions appearing or modified by the insertion of this edge. An alternative construction is also proposed, adding all the edges of the same generation at once, and defining the functions of the new level as those appearing or modified after the insertion of all edges. The theoretical optimality for this alternative construction, that we name macro decomposition, follows from the previous one, but the numerical results show an improved performance.
The paper is organized as follows. In Section 2, we introduce the framework of parallel subspace correction (PSC) method and present its convergence theory based on the two properties mentioned above. We briefly review the basics of univariate/multivariate B-splines in Section 3. In Section 4, we give a new definition of T-meshes by bisections and discuss -admissible T-meshes and their fundamental properties as in [Morgenstern_Peterseim]. In Section 5, we construct a space decomposition on -admissible T-meshes and then prove that the two aforementioned properties are satisfied. In Section 6 we obtain the optimality result for the BPX preconditioner, and also deal with the macro decomposition, showing that is also an efficient space decomposition for the purpose of implementation. Some numerical results that validate our theory are presented in Section 7.
2.1 Problem setting
We are interested in the second order Laplacian with homogeneous Dirichlet boundary conditions,
where denotes the boundary of and . The isogeometric approximation to the solution of (1) is the function such that
and is the isogeometric discrete space. Defining a linear operator by
and also by , we have to solve the linear operator equation
for some .
In the rest of the paper, we will adopt the following compact notation. Given two real numbers we write , when for a generic constant independent of the knot vectors and the mesh size but depending on the spline degree and the geometric map (defined below), and we write when and .
2.2 Preconditioned conjugate gradient method
Let be a symmetric positive definite (SPD) operator. Applying it to both sides of , we get an equivalent equation
The preconditioned conjugate gradient (PCG) method can be viewed as a conjugate gradient method applied to solving where is called a preconditioner (see, e.g., [CNX, Xu1992] for an extensive description).
Let , be the solution sequence of the PCG algorithm. Then the following error estimate is well-known:
which implies that the PCG method converges faster with a smaller condition number .
2.3 The method of parallel subspace corrections
The method of parallel subspace corrections (PSC) provides a particular construction of the iteration operator . The starting point is a suitable decomposition of :
where are subspaces of . The model problem (1) can be split into sub-problems in each with smaller size. Throughout this paper, we use the following operators, for :
the projection in the inner product ;
the natural inclusion;
the projection in the inner product ;
the restriction of to the subspace ;
an approximation of .
This method performs the correction on each subspace in parallel, with the operator defined by
where denote the adjoint of with respect to . It is readily to check that and with .
The convergence analysis of (PSC) is based on the following two important properties [CNX]:
(A1) Stable Decomposition For any , there exists a decomposition such that
(A2) Strengthened Cauchy-Schwarz (SCS) inequality For any
For a space decomposition satisfying both properties, one can prove the following result on the preconditioned linear system [CNX]:
Let be a space decomposition satisfying (A1) and (A2), and let be SPDs for such that
Then defined by (2) is SPD and
The goal of the paper is the construction of a preconditioner like (2) for T-splines, and the proof that (A1) and (A2) are satisfied.
In this section we recall the definition and main properties of B-splines, mainly to fix the notation. For a more extensive description on the use of splines in isogeometric analysis, the reader is referred to [Hughes_Cottrell_Bazilevs, IGA-book], and to [IGA-acta] for a mathematical perspective.
3.1 Univariate B-splines
Given two positive integers and , we say that is a -open knot vector if
where repeated knots are allowed. From the knot vector , univariate B-spline basis functions of degree are defined recursively using the Cox-De Boor formula (see [DeBoor]). The definition of each B-spline , is determined by a local knot vector . Whenever necessary, we will stress it by adopting the following notation:
Thus, the basis function has support
Throughout the paper, we assume that the maximum multiplicity of the internal knots is less than or equal to the degree , that is, the B-spline functions are at least continuous.
Let us select from a subset of non-repeated knots, or breakpoints, with , . We point out that the local mesh size of the element is called for . Moreover, to the element , that can be written as for a certain , we associate the support extension defined by
We denote by
the univariate B-splines space spanned by those B-splines of degree . The functions are a partition of unity, as shown in [Schumi].
Following [Schumi, Theorem 4.41], we define suitable functionals , for , which are dual to the B-splines basis functions, that is
where is the Kronecker delta. The following estimate of the functionals will be useful in the sequel.
If , with then
Proof. The proof can be found in [Schumi, Theorem 4.41].
We note that these dual functionals are locally defined and depend only on the corresponding local knot vector, namely,
Let be the projection that is built from the dual basis as detailed in [Schumi, Theorem 12.6], that is,
is locally quasi-uniform, that is, there exists a constant such that the mesh sizes satisfy the relation , for .
For any non-empty knot span , we have
where the constant depends only on the degree . Moreover, if Assumption 3.1 holds, we have
where the constant depends only on and and denotes the classical Sobolev norm.
Proof. The proof can be found in [IGA-acta, Proposition 2.2].
3.2 Multivariate splines
Multivariate B-splines can be constructed by means of tensor products. We discuss here the bivariate case, the higher-dimensional case being analogous.
For , assume that , the degree and the -open knot vector are given. We set the polynomial degree vector and . We introduce a set of multi-indices and for each multi-index , we define the local knot vector
Then we can introduce the set of multivariate B-splines
The spline space in the parametric domain is then
We also introduce the set of non-repeated interface knots for each direction , , which determine the intervals , for . These intervals lead to the Cartesian grid (or simply ) in the unit domain , also called the Bézier mesh:
For a generic element , we also define its support extension
with the univariate support extension by (5).
3.2.1 Multivariate quasi-interpolants
The main drawback of B-splines is their tensor-product structure, which prevents local refinement as required by adaptive methods. One of the alternatives is the use of T-splines [Sederberg_Zheng], a superset of B-splines that allows for local refinement. The use of T-splines in IGA was first explored in [Bazilervs_Calo_Cottrell_Evans, Dorfel_Juttler_Simeon], and it has led to a growing interest for the analysis of their mathematical properties. In this section we are collecting mathematical results from [LZSHS12, BBCS12, BBSV2013] (linear-independence, dual basis and projectors), [Li-Scott, Bressan_Buffa_Sangalli] (nestedness and space characterization) and [IGA-acta, Section 7] (local linear independence), following mainly the notation in [IGA-acta].
We restrict ourselves to T-splines where refinement is always performed by bisection, and the multiplicity is never reduced. A more general setting can be considered, but it would only add technical difficulties without adding more insight.
4.1 T-mesh defined by bisection
As in the previous section, let us assume that we are given the degrees , the integers and the open knot vectors , and let us denote for . For simplicity we assume that the internal knots in are not repeated and equally spaced, so the element size in each parametric direction can be denoted by and . Our starting point is the index mesh in the index domain , defined as the Cartesian grid of unit squares
which is associated to the tensor-product B-spline space .
To define T-splines by bisection, we start introducing, for any integer and for , the set of rational indices
and we notice that if . We also define the ordered knot vectors
in a recursive way: starting from , for and for any new index , we define the knot
which is well defined because . Clearly, for , and the interval size is . Notice that in this procedure we do not introduce new knots between the repeated knots of the open knot vector.
We also define, for an arbitrary rectangular element in the index domain , with indices , the following bisection operators (see [Morgenstern_Peterseim, Definition 2.5]):
Notice that the bisection operators split the element in two adding a vertical and a horizontal edge, respectively, only when the corresponding length in the parametric domain is greater than zero, that is, when .
Starting from the Cartesian grid , we define a T-mesh by successive applying bisection of elements, in the form
where we use the formal addition (see [Morgenstern_Peterseim, Definition 2.6] and also [CNX, Section 3])
and the bisection operator can be either or . Moreover, we define the finest level as the minimum integer such that the th coordinate of all vertices in the T-mesh belongs to , for .
We notice that a T-mesh defined with this procedure automatically satisfies the first condition in the definition of admissible mesh in [IGA-acta, Definition 7.10], that is, the first lines closer to each boundary are completely contained in the mesh. Instead, the second condition of not having T-junctions in the so-called frame region is not satisfied, because T-junctions may appear on the interface between the frame and the active region. In any case, these T-junctions do not affect the definition of the T-spline functions, and the mesh can still be considered “admissible”.
4.2 Analysis suitable T-splines
To construct the blending functions associated to a T-mesh , we have to define first the set of anchors, that we denote by , and that depends on the degree, see [IGA-acta, Definition 7.13]. These are either the set of vertices ( and odd), elements ( and even), horizontal edges ( even, odd) or vertical edges ( odd, even) in the active region, which is the rectangle
see the examples in Figure 1. We also define the horizontal (resp. vertical) skeleton of the mesh, and denote it by (resp. ), as the union of all horizontal (resp. vertical) edges and vertices. The union of and will be called skeleton. Then, for each anchor we construct an ordered horizontal index vector of indices, tracing a horizontal line from the anchor and collecting the closest indices leftwards and rightwards where the line intersects the vertical skeleton of the mesh, plus the index of the anchor if the degree is odd, see [IGA-acta, Definition 7.14] and Figure 1. A vertical index vector of indices is constructed in an analogous way tracing a vertical line passing through the anchor.
Then, given an anchor with horizontal and vertical index vectors
we define the local knot vectors
and from these local knot vectors we define the associated bivariate function
We will denote by the space generated by the T-splines.
A sufficient condition to guarantee linear independence of the T-spline functions is the dual compatibility condition, introduced in [BBCS12, BBSV2013] and equivalent to the analysis-suitability condition in [Li-Scott]. The former was generalized to arbitrary dimension in [IGA-acta]. We say that two local knot vectors overlap if both of them can be written as subvectors, with consecutive indices, of the same global knot vector. Then we say that the T-mesh is dual compatible if for each pair of anchors , with , there exists a direction such that the local knot vectors and are different and overlap. See [IGA-acta, Section 7.1] for details.
Assuming that the T-mesh is dual compatible, then the functionals
form a dual basis for (see [IGA-acta, Proposition 7.3]). Moreover, we can build the projection operator by
For the analysis of the projector properties, we make use of the Bézier mesh, that we define as in [IGA-acta, Section 7.3]. The Bézier mesh is different from the T-mesh, and plays a similar role to the mesh in finite elements. We start recalling that T-junctions are internal vertices of the T-mesh with valence equal to three, that can be grouped in horizontal () and vertical () T-junctions. For a T-junction of type with index coordinates , we define the extension as the minimal horizontal line that, passing through the T-junction, intersects (closed) vertical edges to its left, and to its right. The extensions in the other cases are defined in a similar way, using for vertical T-junctions. Then, we define the extended T-mesh in the index domain adding to all the T-junction extensions. The Bézier mesh (or simply ) associated to is then defined as the collection of non-empty elements in the domain of the form
as depicted in Figure 2.
For a Bézier element , and in general for any subdomain , we define the support extension as the union of the supports of the functions whose support intersects , that is
where int() denotes the interior of a set . Moreover, we define as the smallest rectangle in containing . The following result holds (see [IGA-acta, Proposition 7.7]):
Let be a dual compatible T-mesh. Then there exists a constant , depending only on , such that for any Bézier element the projector (11) satisfies
Finally, we notice that for each anchor the index vectors (10) define a local Cartesian grid of cells, called tiled floor in [BBCS12, BBSV2013]. Moreover, we also define the parametric tiled floor, as the set of non-empty cells
We remark that in general, the cells of the tiled floor do not coincide with the elements of the T-mesh, and the cells of the parametric tiled floor do not coincide with the elements of the Bézier mesh.
4.3 Analysis suitable T-splines by bisection
In order to apply BPX preconditioners to analysis-suitable T-splines, it is necessary to define a suitable refinement procedure that provides a multilevel structure. In this section we adopt the refinement strategy introduced and analysed in [Morgenstern_Peterseim], and present a new local quasi-uniformity result necessary for the analysis of multilevel preconditioners. The idea in that work is to refine by bisection alternating the refinement direction, and whenever a new edge is added, a recursive algorithm is called to refine the elements in the neighborhood, ensuring that the condition of analysis-suitability is preserved. One of the advantages of this refinement algorithm is that it allows to associate a level (or generation) to each element and function, as required by multilevel methods.
and setting the generation for (see [Morgenstern_Peterseim, Definition 2.6]). Moreover, we say that the bisection has generation . Without loss of generality (see also the proof of Theorem 3.6 in [Morgenstern_Peterseim]), in the following we will assume that the bisections in (8) are ordered by their generation, that is, if then , and the same relation holds for the elements generated by the two bisections. In the following, we will denote by the finest generation (or level), that is, the generation of the last bisection.
We recall that elements with zero length in one parametric direction are not bisected in that direction, see Section 4.1. However, when applying the bisection operator their generation is increased by one, in such a way that the next time the bisection operator is applied, they will be refined in the other direction.
To define the generation of the functions, we denote by the functions of the tensor product spline space in the coarsest mesh, and then define for the collection of T-spline functions newly appeared or modified after the bisection as
An example of the definition of is shown in Figure 3. Functions in have generation 0, whereas functions in have the same generation of the bisection . To alleviate notation, in the following we will denote by the generation (or level) of functions in , that is
with the convention that . Notice that the subscript varies from to , while takes values from 0 to .
To each , we associate a subspace
We also define the support of functions in and its support extension as
Notice that, since we alternate the directions of refinement, and recalling the notation of Section 4.1, the local knot vectors of a function of generation are contained in
For convenience, we also introduce the corresponding set of rational indices, , and the Bézier mesh in the parametric domain . For each generation we have that the mesh size is . This important relation between generation and mesh size can be also represented using additional notation as
Apart from the generations, that are necessary to provide the multilevel structure, we also need some definitions to adapt the refinement algorithm from [Morgenstern_Peterseim] to the case of having an open knot vector. Given two points we define the distance between them componentwise as the vector
Moreover, let us define, for a point in the index domain , its translated version given by
Given a bisection T-mesh and an element , we denote its middle point as . We can define the distance of a point to the element , and the distance between two elements and as
respectively. Then, we define the -neighborhood of as (see [Morgenstern_Peterseim, Definition 2.4] and Figure 4)
We also need to define the set (see [Morgenstern_Peterseim, Corollary 2.15])
Given a bisection T-mesh and an element , we say that the bisection of is -admissible or simply admissible, if all satisfy . Moreover, we say that a bisection T-mesh is -admissible, if it can be obtained as in (8) with a sequence of admissible bisections, see the example of Figure 4. It has been proved in [Morgenstern_Peterseim, Theorem 3.6] that admissible T-meshes are also analysis suitable, and therefore the dual basis and the projector of the previous section can be built.
The result in [Morgenstern_Peterseim] does not take into account the repeated knots of the open knot vector. However, the same ideas apply using the definitions above, because the bisection of zero measure elements only adds new lines in the “safe” direction, without causing intersection of T-junction extensions. The use of the translated points in practice forces that a line arriving at a repeated knot, which is on the boundary of the parametric domain, should continue until the boundary of the index domain.
Besides the analysis suitability condition, we need a local quasi-uniformity result, for which it is necessary to use the following auxiliary lemmas.
Let be an anchor of a -admissible T-mesh . Then, any cell in the parametric tiled floor of contains at most two Bézier elements of .
Proof. Given a cell of the tiled floor, from [BBSV2013, Lemma 3.2] it does not contain any vertex of in its interior, and any line of the extended mesh in its interior belongs to a T-junction extension. Since the mesh is analysis suitable, only vertical or horizontal lines can be found in its interior, but not both. Since a cell in the parametric tiled floor corresponds to a cell in the tiled floor, the result holds because we are refining by bisection alternating the refinement directions.
Given a -admissible T-mesh and an element , for any it holds .
Proof. See [Morgenstern_Peterseim, Lemma 2.14]
Let be an anchor of a -admissible T-mesh , associated to a function of generation . Then the length of the cells of the parametric tiled floor of in each parametric direction is equal to