Multigrid methods for isogeometric discretization.

We present (geometric) multigrid methods for isogeometric discretization of scalar second order elliptic problems. The smoothing property of the relaxation method, and the approximation property of the intergrid transfer operators are analyzed. These properties, when used in the framework of classical multigrid theory, imply uniform convergence of two-grid and multigrid methods. Supporting numerical results are provided for the smoothing property, the approximation property, convergence factor and iterations count for V-, W- and F-cycles, and the linear dependence of V-cycle convergence on the smoothing steps. For two dimensions, numerical results include the problems with variable coefficients, simple multi-patch geometry, a quarter annulus, and the dependence of convergence behavior on refinement levels [Formula: see text], whereas for three dimensions, only the constant coefficient problem in a unit cube is considered. The numerical results are complete up to polynomial order [Formula: see text], and for [Formula: see text] and [Formula: see text] smoothness.

Introduction

Isogeometric method (IGM), introduced in 2005 [28], aims to bridge the gap between finite element method (FEM) and computer aided design (CAD). The main idea of IGM is to directly use the geometry provided by the CAD system, and following the isoparametric approach, to approximate the unknown variables of differential equation by the same functions which are used in the CAD system. IGM offers several advantages when compared to classical FEM. For example, some common geometries arising in engineering and applied sciences, such as circles or ellipses, are represented exactly, and complicated geometries are represented more accurately than traditional polynomial based approaches. Another noteworthy advantage of IGM over classical FEM is the higher continuity. It is a difficult and cumbersome (if not impossible) task to achieve even inter-element continuity in FEM, whereas IGM offers up to continuity, where p denotes the polynomial order and m denotes the knot-multiplicity.

A primary goal of IGM is to be geometrically precise at the coarsest discretization level. In particular, the description of the geometry, taken directly from the CAD system, is incorporated exactly at the coarsest mesh level. This eliminates the necessity of further communication with the CAD system when mesh refinement is carried out. Thereby, the mesh refinement does not modify the geometry. There are several computational geometry technologies that could serve as a basis for IGM. However, non-uniform rational B-splines (NURBS) are the most widely used and well established computational technology in CAD, which we shall also pursue in this work. In last several years IGM has been applied to a variety of problems, e.g., fluid dynamics, electromagnetics, structural mechanics, etc. with promising results. For a detailed discussion see early papers on IGM [2,8-11,18,19] and the book [17]. Since the introduction, most of the IGM progress has been focused on the applications and discretization properties. Nevertheless, when dealing with large problems, the cost of solving the linear system of equations arising from the isogeometric discretization becomes an important issue. Clearly, the discretization matrix A gets denser with increasing p. Therefore, the cost of a direct solver, particularly for large problems, becomes prohibitively expensive. This necessitates the development and use of fast and efficient iterative solvers. It is known that the performance of iterative solvers depends on the condition number of the matrix A. Let (i.e. ratio of largest to smallest eigenvalues) denote the spectral condition number of A. In Table 1, we present of the Laplace operator. We consider a unit square domain and a uniform mesh of elements (open knot-spans for IGM) with mesh-size h. This also serves as a comparison between FEM with Lagrange basis1 and IGM. For a fair comparison, we take continuity in IGM as this results in the same problem size for both the methods. Though the condition number for both the methods reaches asymptotically, however, for IGM the polynomial order p clearly affects the range of the mesh when asymptotic behavior is reached. For example, for IGM with , the asymptotic behavior is not reached up to a reasonably refined mesh. On one hand, this is an advantage as the condition numbers are moderate towards the finer spectrum of the mesh, but on the other hand, this is a serious disadvantage towards the coarser spectrum of the mesh. Note that the condition number rapidly increases with p, and it can reach ∼109 for even for . This is also reflected by the bound of which behaves like , see [24], where d denotes the dimension of the problem domain.

Table 1

Comparison of .

n_0	p=2		p=5
	FEM	IGM	FEM	IGM
2	14	7	581	11094
4	55	12	2317	12951
8	216	36	9263	13680
16	859	140	37050	13886
32	3434	554	148198	13939
64	13734	2215	592789	13952

To the best of authors’ knowledge, so far there are only very few papers [12,13,16,29] which address the performance of linear algebra solvers. In Ref. [16], the authors study the performance of direct solvers which are clearly not suitable for large problems, specially in three-dimensions. In Ref. [29], the tearing and interconnecting approach of finite element methods is used in the context of isogeometric analysis, and the numerical tests (in absence of any theoretical study) suggest almost optimal (with a logarithmic factor) convergence rates of the proposed isogeometric tearing and interconnecting method. The only paper which provides rigorous theoretical study, supported by extensive numerical examples, is by Beirao et al. [12] where the authors discuss the overlapping Schwarz methods. The same authors have also proposed BDDC preconditioners for isogeometric analysis in [13].

In this paper we address another class of linear algebra solvers with optimal complexity, namely multigrid methods. During the last five decades (first paper by Fedorenko in 1961), these methods have been established as a powerful and efficient tool for solving linear system of equations arising in a variety of problems [5,26,38]. The key idea of multigrid goes back to R.P. Fedorenko in the early 60s [22,23], who developed the first multigrid method for solving the Poisson equation on a unit square. The first rigorous convergence proof was provided by Bakhwalov [4]. In early 70s, the multigrid idea was generalized to variational finite difference equations and general finite element equations by Astrachancev [1] and Korneev [30]. However, the huge potential of multigrid methods was realized due to the works of Brandt [6] and Hackbusch [25,26]. A few years later, in the early eighties, algebraic multigrid methods were introduced by Brandt et al. [7], which rebuild the multigrid algorithm based on the information that is accessible via the system of (linear) algebraic equations only. For a more recent exposition of multigrid methods in a multilevel block factorization framework, see also [39].

Our focus in this paper is on multigrid methods for solving the linear system of equations arising from the isogeometric discretization of scalar second order elliptic problems in a single patch. We first prove the condition number estimates of the discrete system for the h-refinement, and provide the supporting numerical results for all levels of smoothness (from to ). These results suggest the expected behavior from the two-(multi-) grid solver. We then prove both the components of the two-grid solver, namely the approximation property of the inter-grid transfer operators, and the smoothing property of the classical Gauss–Seidel (symmetric as well as non-symmetric) method. Together, these two components establish the h-independence of the two-grid solver. For the multi-grid solver, which uses the two-grid solver recursively, we recall the h-independent convergence estimates from [26].

Following the terminology of traditional FEM, we will call the open knot-span as element wherever appropriate. Moreover, as most of the NURBS based designs in engineering use polynomial order p = 2 and 3, throughout this article we will confine ourselves up to p = 4. Furthermore, throughout this article we use the notation (respectively ) to denote (respectively ) where the constant c is independent of the mesh parameter h and the inequality arguments, but it may depend on the polynomial order p.

The contents of this article are organized as follows. In Section 2 we briefly recall the notations for B-splines and NURBS. The geometry mapping and the function spaces are also introduced there. In Section 3 we describe the model problem and recall error estimates. Furthermore, the properties of the discrete system and the norm equivalences are also studied there. In Section 4 we discuss the two-grid method. The multigrid method is then discussed in Section 5. Numerical results on four model problems are presented in Section 6. Finally, some conclusions are drawn in Section 7.

Notations

To keep the article self-contained, we briefly recall the definitions of B-splines and NURBS. For the properties of B-splines and NURBS, which are related to our problem, the reader is referred to [17,28]. For a detailed exposition see, e.g., [32,34,36]. Let p be a non-negative integer denoting the polynomial order, and n be the number of basis functions (B-splines or NURBS). With , denoting the knot index, we assume that the knot vector is a sequence of non-decreasing knots . The knot vector is uniform if the knots are equally spaced, and it is non-uniform when the knots are unequally placed. It is also possible for more than one knot to have the same value, wherein they are called multiple knots. A knot vector is said to be open if its first and last knot values appear times.

The B-spline basis functions, denoted by , are defined recursively as follows:Note that for non-repeated internal knots the support of a B-spline basis function of order p is always knot spans, and every knot span is shared by B-spline basis functions, see Fig. 1 where we plot B-spline basis functions for open, uniform knot vector {} with order 2 and 8. The basis functions formed from open knot vectors are interpolatory at the ends of the parameter space interval . In general, basis functions of order p have continuous derivatives across knot , where is the multiplicity of the value in the knot vector. When the multiplicity of an internal knot value is exactly p, the basis is interpolatory at that knot. This is an important property of B-spline basis functions, in particular, from analysis point of view. Moreover, in IGM the geometry is fixed at the coarsest level of discretization, and any subsequent refinement (whether or ) does not change it. For example, if a partition of is given with the knot vector , then the refined partition can be obtained from via a regular subdivision of knot vector into , where . Further refinements are similarly carried out.

Fig. 1

B-spline functions for open uniform knot vector.

Given n basis functions , and corresponding control points , a piecewise polynomial B-spline curve is given byUsing a tensor product of one-dimensional B-spline functions, a B-spline surface is defined as follows:where , denote the control points, is the tensor product of B-spline basis functions and , and and are the corresponding knot vectors.

Furthermore, let be a set of control points for a projective B-spline curve in with knot vector . For the desired NURBS curve in , the weights and the control points are derived by the relationswhere is called the ith weight and is the dth-dimension component of the vector . Now let the weight function w be defined asThen, the NURBS basis functions and curve are defined byThe NURBS surfaces are analogously defined as follows, see e.g. [17,28] for details,where is the tensor product of NURBS basis functions and . For details related to the B-spline (NURBS) solids (), see e.g. [17,28].

To deal with the tensor-product structure in d-dimensions, we introduce the dimension index set , and the index set for knot vectors . Also, let be the index set of number of basis functions in each dimension, and and be the index set of polynomial order and number of basis functions, respectively, in all dimensions. Now let be an open parametric domain which we will refer as a patch. Assume that d open knot vectors , are given such that and for all . Associated with , we partition the patch into a meshwhere Q is a d-dimensional open knot-span whose diameter is denoted by . We consider a family of quasi-uniform meshes on , where denotes the family index, see [8].

Furthermore, let denote the B-spline space associated with the mesh . Since we do not consider p-refinements, we will use to denote the mesh family for all polynomial orders. The functions in are piecewise polynomials of order in the dth coordinate. Given two adjacent elements and , by we denote the number of continuous derivatives across their common -dimensional face . In the analysis, we will use the following Sobolev space of order where has the usual meaning of ith-order partial derivative, and is the usual Sobolev space of order m. The space is equipped with the following semi-norms and normClearly, for all nested meshes we have for all , where refers to the initial mesh. To a non-empty element we associate the support extensionwhich is the union of supports of those basis functions whose support intersects Q. The restriction of to the support extension is denoted by , and is equipped with the following semi-norms and normThe NURBS space on the patch , associated with the mesh , will be denoted by . When no ambiguity should arise, we will use the notation to represent the polynomial space of either B-splines or NURBS.

Moreover, let the NURBS geometrical map , which is a parametrization of the physical domain , be given by (7) with suitable control points. We assume that is invertible, with smooth inverse, on each element . Therefore, each element is mapped into an element , and the support extension is mapped into . Thereby, in the physical domain we introduce the mesh , where h denotes the maximum element size (hereinafter called the mesh-size) in the domain . Note that the notation h is used for parametric domain as well as physical domain, however, it is a different quantity in both the contexts. Wherever needed, by we will denote the element size in the physical domain. On the physical domain , we denote the space of B-splines by and the space of NURBS by , which are defined asWhen no ambiguity should arise, we will collectively denote and by , and and by , respectively. We will denote the number of elements (open knot-spans with non-zero measure) for a one-dimensional uniform knot vector by . Furthermore, let denote the cardinality of the space . Note that for with order , for all , and continuity, we have .

Finally, we associate a reference support extension to through a piecewise affine map such that each element is the image of a unit hypercube , where . For brevity reasons, we omit further details (including the related spaces) related to the map and refer the reader to [8].

Model problem

Let , be an open, bounded and connected Lipschitz domain with Dirichlet boundary . In this article we consider the scalar second order elliptic equation as our model problem:where is a uniformly bounded function for . Let denote the space of functions which vanish on . By we denote the finite-dimensional spaces of the B-spline (NURBS) basis functions.

Introducing the bilinear form and the linear form asthe Galerkin formulation of this problem reads:

Find such thatIt is well known that (16) is a well-posed problem and has a unique solution.

Error estimates

To keep the article self-contained, we recall some results from [8,36]. By l and m we shall denote integer indices with .

Approximation property of the spline space : The following result is analogous to the classical result by Bramble and Hilbert.

Lemma 1 [8, Lemma 3.1]. Given , the support extension as defined in (10), and , there exists an such that

Projection operators (quasi-interpolants): Let be a projection operator on the spline space , which is defined as follows, see [36, Chapter 12]:where are dual basis functionals defined asThe projection operator on the NURBS space is defined as, see [8],where the weight function w is defined in (5). Collectively, the projection operators and will be denoted by . Finally, the projection operator on the physical space is defined as, see [8],Lemma 2 [36, Theorem 12.6]. The projection operator has the following properties:

Interpolation error estimates: The following lemmas concern the interpolation error estimates.

Lemma 3. Let the projection operator , defined by (18), satisfy (21). Then the following estimate holds for all , see and .For the projection operator the following result is valid for all , see :For the physical domain we have the following result:

Lemma 4 [8, Theorem 3.2]. For the projection operator , the following estimate holds for all .Note that the constants in (23) and (24) depend on the weight function w (and hence on the shape of the parametric domain).

Now assuming sufficient regularity (for the dual problem), a classical convergence analysis and the duality argument (Aubin-Nitsche’s trick) easily give the following result.

Theorem 5. The solution of the problem (16) satisfies the following error estimates

The discrete system

By approximating and using B-spline (NURBS) basis functions , where , the weak formulation (16) is transformed into a set of linear algebraic equationswhere denotes the stiffness matrix obtained from the bilinear form denotes the vector of unknown degrees of freedom (DOF), and denotes the right hand side (RHS) vector from the known data of the problem. In the following Lemma we show the equivalence of the Euclidean norm and the maximum norm for the B-spline (NURBS) space. In this section, for ease of notations we assume uniform polynomial order in each dimension, i.e. for all , although the results are easily generalizable for non-uniform order case.

Let be the space of B-spline (NURBS) basis functions. Let , where are arbitrary. Then the following relation holds for all

We only consider the non-trivial case, i.e. there exists some i for which . For any , there are at most basis functions with non-zero support. Let denote the index set for the basis functions which have non-zero support in K. Also, let . Invoking the non-negativity and the partition of unity properties of basis functions, we haveFurthermore, since , using the stability of B-Spline basis functions [35,31], we obtain the right hand side inequality with a constant . □

Using Sobolev inequalities, see [8], and following the standard FEM approach, see e.g. [3], we obtain the following bounds on the condition number of the matrix .

Let the basis satisfy (28). Then the following relation holdswhere denotes the spectral condition number of .

From (29), we also note thatwhere denotes the spectral norm.

In Table 2, we present the extremal eigenvalues and the spectral condition number of for . We consider all levels of smoothness, i.e. minimum to maximum for the polynomial orders . For h-refinement (knot insertion), we see that the extremal eigenvalues satisfy the theoretical estimates (29) for the discrete system of second order elliptic problems, i.e. maximum eigenvalues are constant, and the minimum eigenvalues are of asymptotically, see e.g. [15]. As mentioned earlier, for reducing the smoothness we insert multiple knots. Note that, due to a high condition number of the B-Spline basis (see proof of Lemma 1), for a given mesh size a higher polynomial order adversely affects the condition number of the matrix , see [24].

Table 2

, and for . Smoothness from to .

		n_0
		2	4	8	16	32	64
p=2
C^0	λ_max	2.1726	2.5607	2.6436	2.6612	2.6653	2.6663
	λ_min	0.2929	0.2008	0.0726	0.0190	0.0048	0.0012
	κ(A¯_h)	7.4169	12.755	36.405	140.01	555.00	2215.0
C^1	λ_max	1.4222	1.4238	1.4896	1.4951	1.4991	1.4997
	λ_min	0.3556	0.3556	0.2855	0.0756	0.0192	0.0048
	κ(A¯_h)	4.0000	4.0044	5.2173	19.768	78.142	311.58
p=3
C^0	λ_max	2.1297	2.2415	2.2844	2.2961	2.2992	2.2999
	λ_min	0.0284	0.0210	0.0190	0.0085	0.0021	0.0005
	κ(A¯_h)	75.111	106.56	120.34	269.99	1075.4	4297.2
C^1	λ_max	0.8962	1.1705	1.1910	1.2078	1.2129	1.2142
	λ_min	0.0386	0.0386	0.0386	0.0191	0.0048	0.0012
	κ(A¯_h)	23.234	30.346	30.878	63.200	252.68	1008.4
C^2	λ_max	1.0384	1.3698	1.5247	1.5627	1.5720	1.5743
	λ_min	0.0336	0.0464	0.0522	0.0547	0.0191	0.0048
	κ(A¯_h)	30.927	29.509	29.192	28.561	82.102	327.22
p=4
C^0	λ_max	2.1002	2.1105	2.1174	2.1195	2.1200	2.1202
	λ_min	0.0024	0.0019	0.0018	0.0017	0.0012	0.0003
	κ(A¯_h)	881.41	1099.7	1189.1	1214.8	1761.9	7041.3
C^1	λ_max	0.8752	1.0840	1.1452	1.1606	1.1644	1.1654
	λ_min	0.0030	0.0030	0.0030	0.0030	0.0021	0.0005
	κ(A¯_h)	293.90	364.01	384.55	389.72	545.08	2177.9
C^2	λ_max	0.6780	0.9178	0.9847	1.0059	1.0118	1.0133
	λ_min	0.0040	0.0048	0.0051	0.0052	0.0047	0.0012
	κ(A¯_h)	167.95	191.78	193.17	192.21	213.24	842.40
C^3	λ_max	0.9369	1.3334	1.7182	1.8111	1.8311	1.8357
	λ_min	0.0028	0.0050	0.0072	0.0081	0.0085	0.0048
	κ(A¯_h)	339.92	269.23	240.26	222.54	215.00	381.73

Before proceeding further, we need to introduce some more notations which are needed for two-(multi-)grid analysis. Let , denote the level of mesh , and be the associated mesh size. The discrete space of B-spline (NURBS) basis functions at level k is denoted by . We assume that the meshes are nested and that . The mesh-dependent inner product on is defined bywhere and denote the approximation coefficients of functions v and w, respectively, with respect to the basis of . The operator is defined byNote from (31) and (32) that . In terms of the operator , the discrete system (27) can be equivalently written aswhere satisfiesSince is symmetric positive definite (SPD) with respect to , we define the following mesh-dependent normwhere denotes the sth-power of the SPD operator for any . Note that the norm coincides with the energy norm . Moreover, . For the equivalence of the norm with the -norm we have the following result.

For we have

Let , where are arbitrary. For any , there are at most basis functions with non-zero support. Let and the index set be as defined in Lemma 1. Also, let . Using the positivity and the partition of unity properties of the basis functions, we know that . Therefore, using , we haveFor the right hand side inequality we havewhere the equivalence constant, say , is , see [31,35]. The result then follows by using . □

Note that the equivalence constant is of the same order as in Lemma 1, and can be improved up to , see [35] for details.

To bound the spectral norm of the matrix we proceed as follows. For SPD matrices we know that the eigenvalues can be estimated in terms of the Rayleigh quotients. Therefore, using the norm (35), the norm-equivalence relation (36), and the inverse inequality , we obtain

Two-grid analysis

In this section we present a two-grid analysis for solving the linear system (27). The purpose of this analysis is to show that the rate of convergence of the two-grid method for IGM is independent of the mesh-size h.

In a two-grid method, the solution of the system (27) is first approximated on the fine grid using a simple stationary iterative method (e.g., Jacobi or Gauss–Seidel), which is often referred to as relaxation process (or smoother because it smooths the error). Then, since on a coarser grid the smooth error can be well represented, and computations are cheaper, the resulting residual equation is transferred to the coarse grid and an error correction (by solving the residual equation) is computed. This error correction is then transferred back to the fine grid where it is added to the approximate solution obtained by the relaxation process. This is called the coarse-grid correction step. Finally, post-relaxation helps to further improve the fine-grid approximation by smoothing error components that may have been contaminated during the inter-grid transfer (from the coarse to the fine grid). The convergence rate of any two-grid method like this depends on the efficiency of the relaxation method (smoother) and on the approximation properties of the inter-grid transfer operators, and on how well smoothing and coarse-grid correction complement each other.

For the two-grid analysis, we shall use the conventional notations h and H to denote the mesh size at the fine level and the coarse level, respectively. Together with the space of basis functions V, the SPD operator A, and the linear functional f, these notations shall be used to reflect the mesh level.

Let be the identity matrix and be the smoothing iteration matrix. Furthermore, let be the orthogonal projection operator (called restriction operator) with respect to , i.e.Another projection operator , called prolongation operator, is analogously defined. We know that the convergence of the two-grid method depends on the iteration matrix [26]where , and and denote the number of pre- and post-smoothing steps, respectively. For simplicity sake (only in analysis), we take . Then, for , the Eq. (39) can be written as [26]This break-up into two separate parts, and , greatly helps the convergence analysis, see [26]. The factor is related to the approximation property and the factor is related to the smoothing property.

In the following two sections we study the approximation property of the inter-grid transfer operators, and the smoothing property of the relaxation method. The h-independent convergence of the two-grid method, i.e.where is defined in (53), is then an immediate consequence of (44), (54) and (57).

Approximation property

To establish the approximation property we first prove the following Lemma, see e.g. [14].

Let . Then the following estimates hold for all .

Using the triangle inequality, we haveThe inequality (42a) is then easily obtained by the equivalence of discrete norms and their continuous counter-parts, using (26), and noting that for quasi-uniform nested meshes. For (42b) we proceed as follows.which gives the desired result. □

Combining (42a) and (42b) we getHence, the quality of approximation of by , where , can also be measured in terms ofEquivalently, albeit in a different terminology, see [14,26] for details, the estimate (43) readsIn Table 3, we present the spectral norm of , which confirms the estimate (44).

Table 3

Illustration of the approximation property, i.e. .

p	n_0
	8	16	32	64
2	2.8125	2.8125	2.8125	2.8125
3	19.1435	18.2758	17.9280	17.8227
4	139.6540	122.8700	117.4090	116.4410

Smoothing property

In this section we recall the smoothing property of the symmetric Gauss–Seidel method.

Let be the decomposition of the matrix , where denotes the diagonal matrix formed from the diagonal of , and and denote strictly lower and strictly upper triangular matrices, respectively. From it follows that .

Now consider the symmetric Gauss–Seidel iterationwhere the preconditioner is given byand the iteration matrix is given byIt is easy to see thatNote that if is SPD (denoted by since for all ) then the matrix and the preconditioner are SPD, and we have the estimateMoreover,because by using a Poincare type inequality on the domain of characteristic size , it can be shown that . Note that the inequality constant also depends on the stability constant and . Therefore, using we getWe also note that , where and denote the entries of the matrices and , respectively, and c is the maximum number of non-zero entries per row (which depends on the polynomial order p). Similarly, it can be shown that . Therefore, using , we getFrom (50) and (51) we getWe are now in a position to prove the following lemma.

LetThe symmetric Gauss–Seidel method (45) satisfies the smoothing propertywhere the function as .

Let . From (48) it follows that . Also, from [26, Lemma 6.2.1] we have for . Hence, using (52) we obtainwhich completes the proof.  □

For the non-symmetric (forward) Gauss–Seidel method, with , we proceed as follows.

Let be a matrix norm corresponding to a vector norm. Let be the iteration matrix of the smoother, and be some matrix. AssumeThen for the following smoothing property holds

We have , and . Therefore,Therefore, . Now using for some matrix with (from Reusken’s Lemma [33], see also [27, Theorem 10.6.8,.6.9]), we get the desire result. □

In Table 4, we list the spectral norm of for , symmetric Gauss–Seidel iterations, which confirms the estimate (55). To compare the smoothing property of symmetric Gauss–Seidel iterations with forward Gauss–Seidel iterations, since the latter is practically advantageous, in Table 5, we list the spectral norm of for , forward Gauss–Seidel iterations. As one might expect, the effect of one symmetric (one forward followed by one backward) Gauss–Seidel iteration is almost the same as for two forward Gauss–Seidel iterations. In fact, we see that for higher p and smaller , the forward Gauss–Seidel iterations perform better than the symmetric version. Due to this reason, we will use forward Gauss–Seidel iterations in our numerical tests for multigrid convergence.

Table 4

Illustration of the smoothing property, i.e. , for symmetric Gauss–Seidel method, .

ν	n_0
	p=1				p=2
	8	16	32	64	8	16	32	64
1	0.3523	0.3789	0.3889	0.3915	0.1317	0.1534	0.1597	0.1620
2	0.1312	0.1468	0.1516	0.1535	0.0462	0.0596	0.0622	0.0632
3	0.0856	0.0894	0.0929	0.0941	0.0181	0.0346	0.0376	0.0388
4	0.0563	0.0662	0.0669	0.0678	0.0071	0.0266	0.0276	0.0280
	p=3				p=4
1	0.1948	0.1947	0.1947	0.1947	0.3775	0.3878	0.3904	0.3911
2	0.0530	0.0521	0.0520	0.0520	0.0917	0.0918	0.0918	0.0918
3	0.0253	0.0251	0.0251	0.0251	0.0364	0.0360	0.0360	0.0360
4	0.0158	0.0157	0.0163	0.0164	0.0222	0.0222	0.0222	0.0222

Table 5

Illustration of the smoothing property, i.e. , for forward Gauss–Seidel method, .

ν	n_0
	p=1				p=2
	8	16	32	64	8	16	32	64
1	0.8917	0.9508	0.9674	0.9716	0.3508	0.3817	0.3946	0.3982
2	0.3496	0.3783	0.3888	0.3915	0.1262	0.1525	0.1595	0.1619
3	0.2007	0.2134	0.2206	0.2229	0.0738	0.0861	0.0905	0.0919
4	0.1314	0.1466	0.1516	0.1535	0.0447	0.0599	0.0622	0.0632
5	0.1065	0.1138	0.1153	0.1168	0.0260	0.0447	0.0477	0.0481
6	0.0862	0.0895	0.0930	0.0941	0.0147	0.0348	0.0377	0.0389
7	0.0697	0.0760	0.0783	0.0788	0.0082	0.0305	0.0323	0.0324
8	0.0561	0.0666	0.0671	0.0678	0.0045	0.0267	0.0277	0.0281
	p=3				p=4
1	0.4897	0.4918	0.4945	0.4951	0.6895	0.7160	0.7218	0.7230
2	0.1758	0.1731	0.1729	0.1729	0.2766	0.2833	0.2843	0.2845
3	0.0868	0.0856	0.0854	0.0854	0.1240	0.1247	0.1257	0.1260
4	0.0510	0.0502	0.0501	0.0501	0.0743	0.0730	0.0730	0.0730
5	0.0342	0.0333	0.0332	0.0332	0.0486	0.0483	0.0483	0.0483
6	0.0249	0.0243	0.0242	0.0242	0.0349	0.0345	0.0345	0.0345
7	0.0193	0.0190	0.0190	0.0194	0.0269	0.0263	0.0263	0.0263
8	0.0159	0.0156	0.0163	0.0164	0.0214	0.0212	0.0211	0.0211

Multigrid convergence

In this section we summarize some important consequences of the smoothing and approximation properties on the convergence of the classical multigrid algorithm in the setting of the isogeometric discretization. Since the proofs of the quoted convergence results can be found in [26], we confine ourselves to a short discussion without repeating any proofs.

For convenience we first consider the symmetric case which is the simplest to analyze. Letdenote the stiffness matrix and the interpolation matrix at level k, respectively, where . Further, let the coarse grid matrix satisfy the Galerkin relationMoreover, assume that the preconditioner is SPD, i.e.and that the smoothing iteration at level k is defined via the iteration matrixThen the iteration matrix of the classical multigrid algorithm with pre- and post-smoothing steps at level k can be recursively defined viawhere , cf. [26, Lemma 7.1.4]. Note that the choices and in (62) correspond to the classical V-cycle and W-cycle multigrid methods, respectively.

W-cycle convergence

Consider the iteration matrix (62) of the W-cycle method, i.e. the case . Further, for convenience, let on all levels k where . Then the following convergence result holds true, cf. [26, Theorem 7.2.3].

Convergence of W-cycle

Let (58)–(61) hold, and the approximation property (44) be satisfied on all levels with a constant , i.e.If and , the contraction number of the W-cycle method () with pre- and post-smoothing steps can be estimated byOtherwise, the smallest root of satisfiesfor all .

V-cycle convergence

Next, consider the iteration matrix (62) of the V-cycle method, i.e. the case . For the case of equal number of pre- and post-smoothing steps, i.e. , we have the following convergence estimate for the V-cycle, cf. [26, Theorem 7.2.2].

Convergence of V-cycle

Under the assumptions of Theorem 7 the V-cycle method () is convergent. In the case its contraction number can be estimated by

For the more general case of pre- and post-smoothing steps, see [26, Theorem 7.2.5]. The numerical results in the next section indicate, however, that these estimates are somewhat pessimistic, and that one obtains better convergence rates in practice.

Numerical results for multigrid convergence

To test the multigrid solvers’ performance, we consider the following test problems, whose discretizations are performed using the Matlab toolbox GeoPDEs [20,21].

Let . Together with , and homogeneous Dirichlet boundary conditions, the right hand side function f is chosen such that the analytical solution of the problem is given by .

Let . Together with , and homogeneous Dirichlet boundary conditions, the right hand side function f is chosen such that the analytical solution of the problem is given by .

The domain is chosen as a quarter annulus in the first Cartesian quadrant with inner radius 1 and outer radius 2, see [20]. Together with , and homogeneous Dirichlet boundary conditions, the right hand side function f is chosen such that the analytic solution is given by .

Let . Together with , and homogeneous Dirichlet boundary conditions, the right hand side function f is chosen such that the analytical solution of the problem is given by .

Furthermore, the operator is chosen such that the coarse basis functions are exactly represented in the space of fine basis functions. At the finest level (largest problem size), the parametric domain is divided into equal elements in each direction. The initial guess for (iteratively) solving the linear system of equations is chosen as a random vector. Let denote the initial residual vector and denote the residual vector at a given multigrid iteration . The following stopping criteria is usedThe average convergence factor reported in the following tables is defined as . In the tables, by M we collectively denote the multigrid method which is specified by the choice of the cycle, i.e. V- or W- or F-cycle. Moreover, , and denote the polynomial order, number of pre- and post-smoothing steps, and the number of mesh refinement levels, respectively. For all the test cases we take the polynomial order .

For Example 1, since the geometry mapping is identity, it suffices to choose the basis functions as B-splines. To evaluate the integrals computationally, we use the Gauss-quadrature formulae with number of quadrature points in each direction. This number is sufficient since the Jacobian from the mapping is constant.

We first present the -dependence of two-grid V-cycle method in Table 6. We see that for a fixed polynomial order p and a fixed mesh size h, the number of iterations of two-grid V-cycle inversely depends on the number of smoothing steps .

Table 6

Poisson problem in a unit square: -dependence of two-grid V-cycle.

ν	n_0
	8		16		32		64
	P	n_it	ρ	n_it	ρ	n_it	ρ	n_it
p=2
1	0.1639	11	0.1869	11	0.1819	11	0.1833	11
2	0.0286	6	0.0320	6	0.0338	6	0.0350	6
4	0.0010	3	0.0009	3	0.0010	3	0.0011	3
8	1.0e−06	2	3.0e−06	2	3.0e−06	2	3.0e−06	2
p=3
1	0.6052	37	0.5864	35	0.5987	36	0.6039	37
2	0.3659	19	0.3494	18	0.3716	19	0.3584	18
4	0.1197	9	0.1195	9	0.1385	10	0.1278	9
8	0.0212	5	0.0172	5	0.0179	5	0.0180	5
p=4
1	0.8790	143	0.8645	127	0.8586	121	0.8598	122
2	0.7763	73	0.7611	68	0.7418	62	0.7392	61
4	0.5487	31	0.5614	32	0.5611	32	0.5502	31
8	0.3293	17	0.3281	17	0.3069	16	0.3043	16

To study the effect of refinement levels on the convergence of V-cycle method, in Table 7, we present the average convergence factor and number of iterations against the number of refinement levels for a fixed , and with and smoothness. As predicted by the theoretical estimates on the optimality of the V-cycle method, it is not surprising to see that and are practically same for all refinement levels. We do not repeat this study for W- and F-cycles, which are also of optimal order and their results for are presented in Table 8.

Table 7

Poisson problem in a unit square: V-cycle convergence, (and ) versus .

ℓ	C^0						C^p-1
	p=2		p=3		p=4		p=2		p=3		p=4
	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it
2	0.0349	6	0.4051	21	0.8143	90	0.0358	6	0.3569	18	0.7420	62
3	0.0349	6	0.4050	21	0.8144	90	0.0358	6	0.3569	18	0.7420	62
4	0.0349	6	0.4050	21	0.8144	90	0.0358	6	0.3569	18	0.7420	62
5	0.0349	6	0.4050	21	0.8144	90	0.0358	6	0.3569	18	0.7420	62
6	0.0349	6	0.4050	21	0.8144	90	0.0358	6	0.3569	18	0.7420	62
7	0.0349	6	0.4050	21	0.8144	90	0.0358	6	0.3569	18	0.7420	62
8	0.0349	6	0.4050	21	0.8144	90	0.0358	6	0.3569	18	0.7420	62

Table 8

Poisson problem in a unit square: multigrid convergence, .

M(ℓ)	n_0
	8		16		32		64
	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it
p=2,C^0
V(2)	0.0236	5	0.0337	6	0.0340	6	0.0341	6
V(4)	0.0236	5	0.0338	6	0.0340	6	0.0341	6
W(4)	0.0236	5	0.0337	6	0.0340	6	0.0341	6
F(4)	0.0039	4	0.0062	4	0.0062	4	0.0063	4
p=2,C^p-1
V(2)	0.0290	6	0.0351	6	0.0347	6	0.0356	6
V(4)	0.0290	6	0.0351	6	0.0347	6	0.0356	6
W(4)	0.0290	6	0.0351	6	0.0347	6	0.0356	6
F(4)	0.0049	4	0.0066	4	0.0065	4	0.0067	4
p=3,C^0
V(2)	0.3762	19	0.3922	20	0.4068	21	0.4043	21
V(4)	0.3761	19	0.3922	20	0.4067	21	0.4043	21
W(4)	0.3762	19	0.3922	20	0.4068	21	0.4043	21
F(4)	0.2335	13	0.2506	14	0.2595	14	0.2571	14
p=3,C^p-1
V(2)	0.3589	18	0.3468	18	0.3465	18	0.3546	18
V(4)	0.3589	18	0.3468	18	0.3465	18	0.3546	18
W(4)	0.3589	18	0.3468	18	0.3465	18	0.3546	18
F(4)	0.2150	12	0.2042	12	0.2040	12	0.2111	12
p=4,C^0
V(2)	0.8101	88	0.8145	90	0.8122	89	0.8139	90
V(4)	0.8103	88	0.8147	90	0.8122	89	0.8140	90
W(4)	0.8103	88	0.8147	90	0.8122	89	0.8139	90
F(4)	0.7299	59	0.7353	60	0.7315	59	0.7343	60
p=4,C^p-1
V(2)	0.7494	64	0.7679	70	0.7278	58	0.7387	61
V(4)	0.7493	64	0.7679	70	0.7278	58	0.7387	61
W(4)	0.7493	64	0.7679	70	0.7278	58	0.7387	61
F(4)	0.6496	43	0.6736	47	0.6220	39	0.6358	41

In Table 8, we present the average convergence factor and number of iterations for V-, W-, and F-cycle multigrid methods. The mesh size varies from to in each direction. We consider both the extreme cases of smoothness, namely, and . As all the cycles are of optimal order, to present a comparative study of all the cases in a concise manner, we consider here only . We make the following observations.

For all polynomial orders, all the approaches exhibit optimal convergence with respect to the mesh refinement, which confirm the theoretical estimates (41) for two-grid method, and (64)–(66) for multigrid methods.

For a fixed mesh size, since the condition number rapidly increases with increasing polynomial order, this affects the two-(multi-) grid convergence.

For smoothness, for any given polynomial order, the convergence factor is slower (for requires more number of iterations) as compared to the problem with smoothness. This phenomenon, which is more prominent for higher polynomial orders, may be attributed to an increased problem size.

Since the V-cycle method is optimal, we see that the performance of the W-cycle for four-grids, i.e. , is only as good as the V-cycle method. Moreover, there is a consistent improvement of a factor about in the number of iterations in the F-cycle as compared to the number of V-cycle iterations. This compensates the additional computational cost in F-cycle to a good extent.

We now study the performance of V-cycle multigrid solver on a multi-patch geometry. This simple model case is produced by p times repetition of the knot at (in both directions). Thereby, we get four patch fully-conforming geometry which has smoothness at interfaces and smoothness elsewhere. The coarsest mesh is fixed with elements in each direction and for both the refinements (2-level and 4-level). The results presented in Table 9 show that the convergence behavior fits nicely between the convergence behavior for global and smoothness, with a bias towards smoothness.

Table 9

Poisson problem in a unit square: V-cycle convergence on a multi-patch geometry, .

p	n_0
	8		16		32		64
	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it
ℓ=2
2	0.0217	5	0.0284	6	0.0353	6	0.0343	6
3	0.3925	20	0.3790	19	0.3727	19	0.3651	19
4	0.8082	87	0.7756	73	0.7558	66	0.7485	64
	32		64		128		256
ℓ=4
2	0.0353	6	0.0343	6	0.0352	6	0.0357	6
3	0.3727	19	0.3651	19	0.3578	18	0.3577	18
4	0.7558	66	0.7485	64	0.7423	62	0.7448	63

We now consider Example 2 with variable coefficients. In Table 10, we present the results for V-cycle multigrid convergence for and . We take the number of quadrature points in both the directions so that the integrals with respect to x-variable are evaluated exactly. However, due to the exponential function, exact integration is not possible with respect to y-variable. We note that the results are qualitatively same as those with constant coefficients case (see Table 8).

Table 10

Variable coefficients elliptic problem in a unit square: V-cycle convergence, .

p	n_0
	8		16		32		64
	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it
C^0
2	0.0177	5	0.0241	5	0.0290	6	0.0322	6
3	0.3162	16	0.3872	20	0.3887	20	0.3910	20
4	0.8005	83	0.7977	82	0.8104	88	0.8121	89
C^p-1
2	0.0342	6	0.0199	5	0.0306	6	0.0357	6
3	0.3067	16	0.3737	19	0.3556	18	0.3516	18
4	0.8146	90	0.7870	77	0.7257	58	0.7260	58

We now consider Example 3 with curved boundary. The geometry for this example is represented by NURBS basis functions of order 1 in the radial direction and of order 2 in the angular direction, see [20]. Since the Jacobian of the geometry mapping is no more a constant, for exact integral evaluations it does not suffice to take the number of Gauss quadrature points in each direction (which is clear from simple heuristic arguments). Therefore, we choose . From numerical experiments, it is found that this is sufficient (for up to ) to keep the approximation error (44) smaller than the -norm of the discretization error (which otherwise would contaminate the accuracy of two-(multi-) grid solver). Note however that this is not detrimental to the optimality of any of the methods, which can be seen from the variable coefficients case presented in Table 10. In Table 11, we present the -dependence of two-grid V-cycle method. In Table 12, we present the convergence factor and the number of iterations for V-, W-, and F-cycle multigrid methods. The mesh size again varies from to in each direction, and both the extreme cases of smoothness, namely, and are considered. All the results are qualitatively similar to that of Example 1 with square domain.

Table 11

Poisson problem in a quarter annulus: -dependence of two-grid V-cycle.

ν	n_0
	8		16		32		64
	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it
p=2
1	0.1926	12	0.2823	15	0.3052	16	0.3319	17
2	0.0371	6	0.0810	8	0.0931	8	0.1126	9
4	0.0014	3	0.0066	4	0.0087	4	0.0136	5
8	3.0e−06	2	4.3e−05	2	7.5e−05	2	2.3e−04	3
p=3
1	0.5858	35	0.6118	38	0.5977	36	0.6036	37
2	0.3477	18	0.3741	19	0.3575	18	0.3670	19
4	0.1196	9	0.1437	10	0.1277	9	0.1383	10
8	0.0159	5	0.0206	5	0.0181	5	0.0191	5
p=4
1	0.8703	133	0.8594	122	0.8604	123	0.8617	124
2	0.7564	66	0.7384	61	0.7408	62	0.7425	62
4	0.5767	34	0.5475	31	0.5488	31	0.5513	31
8	0.3331	17	0.3046	16	0.3054	16	0.3083	16

Table 12

Poisson problem in a quarter annulus: multigrid convergence, .

M(ℓ)	n_0
	8		16		32		64
	ρ	n_it	ρ	n_it	ρ	n_it	ρ	n_it
p=2,C^0
V(2)	0.0716	7	0.0977	8	0.0985	8	0.1071	9
V(4)	0.0716	7	0.0976	8	0.0985	8	0.1071	9
W(4)	0.0716	7	0.0977	8	0.0985	8	0.1071	9
F(4)	0.0189	5	0.0314	6	0.0325	6	0.0346	6
p=2,C^p-1
V(2)	0.0371	6	0.0810	8	0.0931	8	0.1126	9
V(4)	0.0371	6	0.0810	8	0.0931	8	0.1126	9
W(4)	0.0371	6	0.0810	8	0.0931	8	0.1126	9
F(4)	0.0071	4	0.0225	5	0.0302	6	0.0378	6
p=3,C^0
V(2)	0.3904	20	0.4046	21	0.3977	20	0.4046	21
V(4)	0.3903	20	0.4045	21	0.3975	20	0.4045	21
W(4)	0.3903	20	0.4046	21	0.3976	20	0.4046	21
F(4)	0.2415	13	0.2573	14	0.2556	14	0.2574	14
p=3,C^p-1
V(2)	0.3477	18	0.3741	19	0.3575	18	0.3670	19
V(4)	0.3469	18	0.3741	19	0.3575	18	0.3670	19
W(4)	0.3472	18	0.3741	19	0.3575	18	0.3670	19
F(4)	0.2046	12	0.2311	13	0.2138	12	0.2246	13
p=4,C^0
V(2)	0.8158	91	0.8184	92	0.8232	95	0.8230	95
V(4)	0.8145	90	0.8183	92	0.8229	95	0.8228	95
W(4)	0.8145	90	0.8183	92	0.8229	95	0.8228	95
F(4)	0.7351	60	0.7413	62	0.7461	63	0.7459	63
p=4,C^p-1
V(2)	0.7564	66	0.7384	61	0.7408	62	0.7421	62
V(4)	0.7582	67	0.7384	61	0.7408	62	0.7422	62
W(4)	0.7581	67	0.7384	61	0.7408	62	0.7422	62
F(4)	0.6609	45	0.6355	41	0.6368	41	0.6411	42

Finally, we consider the three-dimensional problem described in Example 4. The results for V-cycle multigrid method are presented in Table 13, which confirm the h-independence and optimality of the solver. The entries marked by represent the cases where the results could not be obtained due to limitation on computational resources. As shown by the results of two-dimensional examples, the W- and F-cycle methods will not offer any improvement in convergence results, and are thus not repeated here.

Table 13

Poisson problem in a unit cube: V-cycle multigrid convergence, .

p	n_0
	8		16		32
	ρ	n_it	ρ	n_it	ρ	n_it
C^0,ℓ=2
2	0.3578	18	0.4073	21	0.4066	21
3	0.8221	147	0.8929	163	0.8947	166
4	0.9879	1514	0.9881	1540	†	†
C^p-1,ℓ=2
2	0.2874	15	0.3383	17	0.3692	19
3	0.8582	121	0.8403	106	0.8431	108
4	0.9728	669	0.9751	731	0.9745	713
C^0,ℓ=4
2	0.3685	19	0.3977	20	0.4076	21
3	0.8891	157	0.8923	162	0.8942	165
4	0.9877	1493	0.9881	1543	†	†
C^p-1,ℓ=4
2	0.3339	17	0.2356	18	0.3700	19
3	0.8572	120	0.8556	110	0.8422	112
4	0.9772	797	0.9738	695	0.9740	698

For all the examples, we also tested the multigrid convergence for intermediate continuities , i.e. , and found that the results lie nicely between the results of and continuities. However, they are not reported here for brevity reasons. We also remark the following on the numerical results of high polynomial orders and where the exact solution has reduced regularity.

It is known from finite elements literature that standard h-multigrid, which is the focus of this article, is not suited for high polynomial order. Most of the literature is for first and second order polynomials only. This fact is related to the smoothing properties of the classical smoothers like Jacobi, Gauss–Seidel or Richardson methods. These methods work effectively only when the error function is oscillatory, whereas the error function gets smoother with increasing polynomial order. For high polynomial order, either p-multigrid should be used or different smoothers should be devised, both of which are beyond the scope of this article. Nevertheless, since IGM in engineering applications mostly utilize second or third order polynomials, in this study we considered polynomial order up to .

In the presence of discontinuities in the coefficients, or due to the irregular geometry (e.g., L-shaped domain), the exact solution of elliptic problems has reduced regularity and lies only in , where depends on the strength of the singularity. Firstly, in such cases the single-patch isogeometric approach with global continuity (for ) is not so attractive. Secondly, the standard (geometric) multigrid methods are not tailored for such general problems and need special treatment. The reduced regularity negatively affects the approximation property of Lemma 4, and thus the overall convergence behavior of solver. Though specific problems can be treated to obtain optimal order convergence (which involves more technical results), however, this is beyond the scope of this article. For such problems, the multi-patch techniques, such as the tearing and interconnecting approach of Kleiss et al. [29] or BDDC approach of Beirao et al. [13], are more suitable where the multigrid solver can be used within each sub-patch.

Conclusions

We have presented multigrid methods, with V-, W- and F-cycles, for the linear system arising from the isogeometric discretization of the scalar second order elliptic problems. For a given polynomial order p, all multigrid cycles are of optimal complexity with respect to the mesh refinement. Despite that the condition number of the stiffness matrix grows very rapidly with the polynomial order, these excellent results exhibit the power of multigrid methods. Nevertheless, this study can only be regarded as a first step towards utilizing the power of multigrid methods in IGM. In our forthcoming work, we will study the multigrid techniques as preconditioners in conjugate gradient method, and also address the Fourier analysis of multigrid methods. Another solver approach with optimal complexity, but with more generality, namely, algebraic multilevel iteration method, is the subject of our current focus for isogeometric discretization of elliptic problems.