Optimal multilevel preconditioners for isogeometric collocation methods

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation xxx (xxxx) xxx www.elsevier.com/locate/matcom

Original articles

Optimal multilevel preconditioners for isogeometric collocation methods Durkbin Cho Department of Mathematics, Dongguk University, Pil-dong 3-ga, Jung-gu, Seoul, 04620, South Korea Received 6 August 2018; received in revised form 12 March 2019; accepted 5 August 2019 Available online xxxx

Abstract We present optimal additive and multiplicative multilevel methods, such as BPX preconditioner and multigrid V-cycle, for the solution of linear systems arising from isogeometric collocation discretizations of second order elliptic problems. These resulting preconditioners, accelerated by GMRES, lead to optimal complexity for the number of levels, and illustrate their good performance with respect to the isogeometric discretization parameters such as the spline polynomial degree and regularity of the isogeometric basis functions, as well as with respect to domain deformations. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Isogeometric analysis; Collocation methods; Multigrid; Preconditioners; GMRES; Multilevel methods

1. Introduction Isogeometric Analysis (IGA), introduced by Hughes et al. in [25], is a novel numerical methodology for partial differential equations (PDEs), for which the same basis functions are used both for the geometric description of the domain and for the numerical approximation of the solution of PDEs, following an isoparametric paradigm. From this perspective, IGA was proposed to fill the gap between computer aided design (CAD) and finite element methods (FEM), by providing an easier and exact representation of the CAD domain and the approximation spaces. Since then, isogeometric techniques have been investigated in a broad scope of application fields: see [6,14] for a more general introduction to IGA. In its original formulation [25], IGA had employed the Galerkin approach, i.e., the discretization of the variational form of the problem. The highly continuous IGA Galerkin solutions exhibit increased accuracy and robustness on a per-degree-of-freedom in comparison to standard finite element Galerkin solutions, but are more costly per-degreeof-freedom, due to the denser structure of isogeometric matrices [2,12,13,26]. Within IGA framework, collocation methods have been recently proposed as an alternative to standard Galerkin approaches. IGA collocation methods are characterized by a significantly reduced computational cost but still guaranteeing higher order convergence rates [3,32]. IGA collocation methods have been extensively used in diverse fields: linear dynamics [4], beam analysis [5,43], Cahn–Hillard problem [24] and contact [16]. As in the Galerkin formulation, collocation stiffness matrix has a condition number that grows quickly with the inverse of mesh size h, see [22,42]. As a consequence, the cost of solving the linear system of equations E-mail address: [email protected]. https://doi.org/10.1016/j.matcom.2019.08.003 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

2

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

arising from the isogeometric collocation discretization becomes an important issue. Therefore, there is currently a growing interest in the design of efficient preconditioners/solvers for IGA Galerkin/collocation methods, in both the mathematical and the engineering communities. Preconditioners and domain decomposition methods [19,35,37] have been successful in IGA Galerkin schemes (for instance, overlapping domain decompositions [39,40], FETItype nonoverlapping preconditioners [28,30], BDDC [41], BDDC deluxe [44,45], BPX [9], multigrid methods [18,21,23]), while there are a handful of results in the literature for preconditioners/solvers or domain decomposition methods for IGA collocation methods (see [18,42]). In the present paper, after constructing our additive and multiplicative multilevel methods for IGA collocation methods, we will study their performance through a variety of numerical tests in two and three spatial dimensions. Since IGA collocation matrices are not generally symmetric, we solve the resulting collocation problems by generalized minimal residual (GMRES) method with BPX and multigrid V-cycle preconditioners. First we not only demonstrate the lower number of iterations to solve the BPX and V-cycle preconditioned systems as compared with unpreconditioned system but also show the dependence of iteration counts on the number J of levels. A particular attention is paid to the optimality of our multilevel preconditioners for IGA collocation problems. Also we investigate the good convergence behavior of the preconditioners with respect to the polynomial degree p and the regularity k. In addition, the good enough results in the presence of domain deformation are reported. The present paper is organized as follows. In Section 2, we introduce the framework of subspace correction methods on IGA Galerkin methods. In Section 3, we briefly review the basics of both B-splines and NURBS, and present a convergence theory of multilevel iteration methods on IGA Galerkin discretization. Section 3 also gives an overview of the isogeometric collocation methods. In Section 4, we describe our proposed multilevel methods for IGA collocation methods, i.e., BPX preconditioners and multigrid V-cycle. In Section 5, we exhibit numerical experiments that justify the optimality of our proposed methods in practice. 2. Problem setting based on Galerkin approach We are interested in the second order elliptic problem with homogeneous Dirichlet boundary conditions, in Ω ⊂ Rd ,

− ∆u = f

on ∂Ω ,

u=0

(1)

where ∂Ω denotes the boundary of Ω and f ∈ L 2 (Ω ). The isogeometric Galerkin approximation to the solution of (1) is the function u ∈ V such that a(u, v) = ( f, v) ∀v ∈ V

(2)

where a(u, v) :=

∫ ∇u · ∇v d x, Ω

( f, v) :=

∫

f v d x, Ω

and V is the isogeometric discrete space that will be defined later. Defining a linear operator A : V → V by (Au, v) = a(u, v), ∀u, v ∈ V and also b ∈ V by (b, v) = ( f, v), ∀v ∈ V, we have to solve the linear operator equation (3)

Au = b

for some u ∈ V. In the rest of the paper, we will adopt the following notation. Given two real numbers a, b, we write a ≾ b when a ≤ Cb for a generic constant C independent of the knot vectors and the mesh size h but depending on the spline degree p, the regularity k and the geometric map F, that will all be defined below. 2.1. Space decomposition and method of subspace corrections Following [46], the starting point is a suitable decomposition of the discrete space V: V=

J ∑

Vi ,

i=0

Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

3

where Vi are subspaces of V. The model problem (1) can be split into sub-problems in each Vi with smaller size. Throughout this paper, we use the following operators, for i = 0, 1, . . . , J : • • • • •

Q i : V → Vi the projection in the L 2 inner product (·, ·); Ii : Vi → V the natural inclusion; Pi : V → Vi the projection in the inner product (·, ·) A := (A·, ·); Ai : Vi → Vi the restriction of A to the subspace Vi ; Ri : Vi → Vi an approximation of Ai−1 ;

From the definitions, we can obtain a simple but important result Ai Pi = Q i A. Parallel subspace correction (PSC) This method performs the correction on each subspace in parallel. In operator form, it reads u k+1 = u k + B(b − Au k ), where B :=

J ∑

Ri Q i .

(4)

i=0

Denoting by Iit the adjoint of Ii with respect to (·, ·), the fact that Ri Q i = Ii Ri Iit implies that B=

J ∑

Ri Q i =

i=0

J ∑

Ii Ri Iit .

i=0

Remark 2.1. ∑ J It is well known (e.g., in [8,10,46]) that if each Ri is a symmetric and positive definite (SPD) operator, then B = i=0 Ri Q i becomes SPD and thus can be used as a preconditioner. Hence the conjugate gradient method can be applied to the preconditioned system (5)

B Au = B f. Among various PSC preconditioners, see Section 3.3 for a description of BPX preconditioners. Successive subspace correction (SSC) This method performs the correction in a successive way. In operator form, it reads v0 = u k , v v

(i+1)/(2J +2)

(J +i+2)/(2J +2)

u

k+1

= v i/(2J +2) + Ii Ri Iit (b − Av i/(2J +1) ), =v

(J +i+1)/(2J +2)

+

I J −i R J −i I Jt −i (b

− Av

i = 0, 1, . . . , J, (J +i+1)/(2J +2)

),

i = 0, 1, . . . , J,

=v . 1

Defining Ti = Ri Q i A : V → Vi , Ti can also be represented as Ri Ai Pi and the error operator for SSC method becomes E sJ = (I − T0 )(I − T1 ) · · · (I − TJ )(I − TJ )(I − TJ −1 ) · · · (I − T0 ). We remark that the SSC algorithm is equivalent to the standard V-cycle multigrid algorithm that is based on a nested sequence of subspaces: V0 ⊂ V1 ⊂ · · · ⊂ V j · · · ⊂ V J , (see Section 5.2.1 or [46, Algorithm 3.8] and the references therein for more details). 3. NURBS-based isogeometric collocation method For the sake of clarity and simplicity, throughout this section we only consider open uniform knot vectors where internal knots have multiplicity one, i.e., the associated B-splines (or NURBS) that will be defined below have global C p−1 regularity. Also, throughout the paper we follow the notation and basic concepts in [42]. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

4

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

3.1. Univariate B-splines and NURBS curves Assume the parametric interval ˆ I := (0, 1). Letting p ≥ 3 be the polynomial degree and j0 be a fixed positive integer, we introduce for any ℓ ∈ N0 := {0, 1, 2, . . . , J } a uniform subdivision by splitting ˆ I into n ℓ = j0 2ℓ 1 −ℓ 1 subintervals of length h ℓ := n = j 2 by uniform dyadic refinement, namely, ℓ

0

ℓ

I := {0, h ℓ , 2h ℓ , . . . , (n ℓ − 1)h ℓ , 1}. On these grids, we define for any integer ℓ ∈ N0 an open knot vector Ξ ℓ with the first and last knots repeated p + 1 times each: Ξ ℓ = Ξ ℓ (I ℓ ) := {ξ1 , . . . , ξ p+1 , ξ p+2 , . . . , ξn ℓ + p , ξn ℓ + p+1 , . . . , ξn ℓ +2 p+1 } = {0, . . . , 0, h ℓ , 2h ℓ , . . . , (n ℓ − 1)h ℓ , 1, . . . , 1}.       p+1 times

p+1 times

ℓ

Given an open knot vector Ξ for ℓ ∈ N0 , univariate B-spline basis functions with maximum smoothness are defined recursively starting with p = 0 (piecewise constants) { 1 if ξi ≤ ξ < ξi+1 0 Nℓ,i (ξ ) = (6) 0 otherwise, for p ≥ 1, ⎧ ξ −ξ ξi+ p+1 − ξ i p−1 p−1 ⎨ Nℓ,i (ξ ) + N (ξ ) p ξi+ p+1 − ξi+1 ℓ,i+1 Nℓ,i (ξ ) = ξi+ p − ξi ⎩ 0

if ξi ≤ ξ < ξi+ p+1

(7)

otherwise, p

where we adopt the convention 0/0 = 0. Thus, the general basis function Nℓ,i has support p

supp(Nℓ,i ) = (ξi , ξi+ p+1 ),

i = 1, 2, . . . , n ℓ + p.

Now we can define the spline space of degree p with maximum smoothness, p Sˆℓ,h ℓ = span{Nℓ,i (ξ ), i = 1, . . . , n ℓ + p}.

(8)

We remark that the spaces are nested for fixed spline degree p, that is, Sˆℓ−1,h ℓ−1 ⊂ Sˆℓ,h ℓ and the number of degrees of freedom roughly doubles in each refinement step. nℓ + p p To obtain a NURBS curve in R2 using basis functions {Nℓ,i (ξ )}i=1 with fixed integer ℓ, we start by introducing the NURBS basis functions of degree p p

p

Nℓ,i (ξ )ωi Nℓ,i (ξ )ωi p = , (9) Rℓ,i (ξ ) = ∑n ℓ + p p w(ξ ) ˆı=1 Nℓ,ˆı (ξ )ωˆı ∑n ℓ + p p where the denominator w(ξ ) = ˆı=1 Nℓ,ˆı (ξ )ωˆı ∈ Sˆℓ,h ℓ is called the weight function. The NURBS curve in Rd is then defined by nℓ + p

C:ˆ I → C(ˆ I ),

C(ξ ) =

∑

p

Rℓ,i (ξ )Ci ,

(10)

i=1

where Ci ∈ Rd are control points (see [20,31,33] for a complete description). 3.2. Multivariate B-splines and NURBS: tensor product By means of tensor products, a multi-dimensional B-spline region can be constructed. We discuss here the case ˆ := (0, 1) × (0, 1) be the twoof a two-dimensional region, the higher-dimensional case being analogous. Let Ω dimensional parametric space. Assume that for d = 1, 2, we are given the degrees pd , the fixed positive integers jd,0 and denote n d,ℓ := jd,0 2ℓ and h d,ℓ := n 1 = j 1 2−ℓ . Now we can define two uniform subdivisions by uniform d,ℓ d,0 dyadic refinements: I 1,ℓ = {0, h 1,ℓ , 2h 1,ℓ , . . . , (n 1,ℓ − 1)h 1,ℓ , 1}, and I 2,ℓ = {0, h 2,ℓ , 2h 2,ℓ , . . . , (n 2,ℓ − 1)h 2,ℓ , 1}. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx ξ

5

ξ

With these grids, we associate open knot vectors Ξ1,ℓ = Ξ1,ℓ (I 1,ℓ ) = {0 = ξ1 , . . . , ξn 1,ℓ +2 p1 +1 = 1} and η η Ξ2,ℓ = Ξ2,ℓ (I 2,ℓ ) := {0 = η1 , . . . , ηn 2,ℓ +2 p2 +1 = 1} and an (n 1,ℓ + p1 ) × (n 2,ℓ + p2 ) net of control points Ci, j . p p One dimensional basis functions Nℓ,i1 (ξ ) and Nℓ,2j (η) with i = 1, . . . , n 1,ℓ + p1 and j = 1, . . . , n 2,ℓ + p2 of degrees η ξ ˆ is then p1 and p2 , are defined from the knot vectors Ξ1,ℓ and Ξ2,ℓ , respectively. The bivariate spline basis on Ω defined by the tensor product construction p ,p

p

p

2 1 1 2 Bℓ,i, j (ξ, η) = Nℓ,i (ξ ) Nℓ, j (η).

η

ξ

Note that the two knot vectors Ξ1,ℓ and Ξ2,ℓ generate a mesh of rectangular elements in the parametric space in a natural way. Analogous to (8), we define p1 , p2 Sˆℓ,h ℓ = span{Bℓ,i, j (ξ, η), 1 ≤ i ≤ n 1,ℓ + p1 , 1 ≤ j ≤ n 2,ℓ + p2 }

(11)

We remark that the spaces are nested for fixed spline degree p, that is, Sˆℓ−1,h ℓ−1 ⊂ Sˆℓ,h ℓ , and the number of degrees of freedom roughly quadruples in each refinement step. ˆ = (0, 1) × (0, 1) Analogously to B-splines, NURBS basis functions on the two-dimensional parametric space Ω are defined as p1 , p2 Rℓ,i, j (ξ, η)

p ,p

1 2 Bℓ,i, j (ξ, η)ωi, j = ∑n 1,ℓ + p1 ∑n 2,ℓ + p2 p , p Bℓ,ˆ1ı,ˆȷ 2 (ξ, η)ωˆı,ˆȷ ˆı=1 ȷˆ=1

p ,p

=

1 2 Bℓ,i, j (ξ, η)ωi, j

w(ξ, η)

,

(12)

where ωi, j = (Ci,ω j )3 and the denominator is the weight function denoted with w(ξ, η). Observe that the continuity and support of NURBS basis functions are the same as for B-splines. NURBS spaces are the span of the basis functions (12). NURBS regions with fixed integer ℓ are defined in terms of the basis functions (12). In particular a single-patch domain Ω is a NURBS region associated with the (n 1,ℓ + p1 ) × (n 2,ℓ + p2 ) net of control points Ci, j , and we ˆ → Ω given by introduce the geometrical map F : Ω F(ξ, η) =

n 1,ℓ + p1 n 2,ℓ + p2 ∑ ∑ i=1

p ,p

1 2 Rℓ,i, j (ξ, η)Ci, j .

(13)

j=1

Following the isoparametric approach, the space of NURBS scalar fields on the domain Ω is defined, component by component as the span of the push-forward of the basis functions (12) p ,p

−1 1 2 Nℓ,h ℓ := span{Rℓ,i, j ◦ F , with i = 1, . . . , n 1,ℓ + p1 ; j = 1, . . . , n 2,ℓ + p2 }.

(14)

Remark 3.1. In one dimensional case, the space of NURBS scalar fields on the NURBS curve C(ˆ I ) is similarly defined as p

Nℓ,h ℓ := span{Rℓ,i ◦ C−1 , with i = 1, . . . , n ℓ + p}. Let us consider two sequences of nested isogeometric spaces: ˆ0 ⊂ V ˆ1 ⊂ · · · V ˆj ⊂ · · · ⊂ V ˆJ , V

and V0 ⊂ V1 ⊂ · · · ⊂ V j · · · ⊂ V J ,

ˆℓ := Sˆℓ,h and Vℓ := Nℓ,h , 0 ≤ ℓ ≤ J are living in parameter space and in physical space, respectively. where V ℓ ℓ Therefore the isogeometric space used in the discretization will be ˆ = SˆJ,h V J

and V = N J,h J .

ˆ on unit domain In other words, the isogeometric space that we consider in the following numerical tests would be V (square in 2D and cube in 3D), and V on the deformed domain. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

6

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

3.3. Convergence analysis of multilevel methods The convergence analysis of subspace correction methods on IGA Galerkin discretization has been carried out by several previous works. · The following result on BPX preconditioned linear system is shown in [9,11]: ∑J – Let V = i=0 Vi be a space decomposition, R0 = A−1 or symmetric Gauss–Seidel 0 , Ri , be either Jacobi∑ J smoothers for 1 ≤ i ≤ J . Then κ(B A) ≾ 1 with BPX preconditioners B = i=0 Ri Q i defined by (4). · The following convergence estimate for multigrid V-cycle is presented in [21]: ∑J – Let V = i=0 Vi be a space decomposition and let the error operators I − Ti be contractions for 0 ≤ i ≤ J . Then we have ∥E sJ ∥ A ≾ 1. The optimality of BPX preconditioners for IGA Galerkin method is shown in [9], while the optimality of BPX preconditioner for IGA collocation methods that will be discussed in the subsequent subsection section is theoretically still missing (see [3–5,32] for details). Moreover the uniform convergence of IGA collocation multigrid V-cycle is still an open issue. The goal of the paper is the design of multilevel iteration methods for isogeometric collocation methods. Taking a closer look at these properties established in IGA Galerkin methods, we observe that the optimal convergence of both preconditioners is independent of the mesh size h J , the number of levels J , but dependent on the polynomial degree p, the regularity k and the geometrical map F. We may expect that the convergence rate of our IGA collocation multilevel preconditioners is comparable to that of their Galerkin counterparts (see the numerical section for more details). Remark 3.2. When the multigrid V-cycle method is used as a preconditioner for conjugate gradient iteration, the proposed multigrid method not only maintains robustness but also leads to significant speedup (see numerical results shown in Table 1 for using the multigrid V -cycle as a preconditioner). Remark 3.3. The multigrid preconditioned conjugate method has the multigrid V-cycle or W-cycle as a preconditioner (see [27,36] for details). It is clear that the W-cycle is a direct consequence of the V-cycle. Thus our results hold true for the multigrid W-cycle preconditioned conjugate method. Although in this paper our attention is paid to the multigrid V-cycle preconditioned conjugate gradient method, we need to clarify that it is not a single method but a representative of a class of methods. 3.4. IGA collocation method In this section, we give a brief review of the isogeometric collocation method for the model problem (1), see [3,4,34] for a complete discussion. To simplify the discussion, we consider only the case of the traditional IGA collocation where collocation points are defined using the Greville abscissae of the knot vectors, see [7]. For fixed integer ℓ, let ξ ℓ,i , i = 1, . . . , n 1,ℓ + p1 , denote the Greville abscissae associated to the knot vector {ξ1 , . . . , ξn 1,ℓ +2 p1 +1 }: ξ ℓ,i := (ξi+1 + ξi+2 + · · · + ξi+ p1 )/ p1 .

(15)

Analogously, let ηℓ, j , j = 1, . . . , n 2,ℓ + p2 , be the Greville abscissae related to the knot vector {η1 , . . . , ηn 2,ℓ +2 p2 +1 }. It is easy to see that ξ ℓ,1 = ηℓ,1 = 0, ξ ℓ,n 1,ℓ + p1 = ηℓ,n 2,ℓ + p2 = 1, while all the remaining collocation points are in (0, 1). For i = 1, . . . , n 1,ℓ + p1 , j = 1, . . . , n 2,ℓ + p2 , we define the collocation points τℓ,i j ∈ Ω by the tensor product structure ( ) ˆ , τℓ,i j = F(ˆ τℓ,i j ) , ˆ τℓ,i j := (ξ ℓ,i , ηℓ, j ) ∈ Ω ( ) ( ) ˆ denotes the closure of the parametric space Ω ˆ , that is, Ω ˆ = [0, 1] × [0, 1]. where Ω Then, the isogeometric collocation problem with Greville abscissae reads: find u h ∈ V such that ⎧ ⎪ ⎨ −∆u h (τi j ) = f (τi j ) i = 2, . . . , n 1,ℓ + p1 − 1, j = 2, . . . , n 2,ℓ + p2 − 1 , u h (τi j ) = 0 (i, j) ∈ ({1, n 1,ℓ + p1 } × {1, . . . , n 2,ℓ + p2 }) (16) ⎪ ⎩ ∪ ({1, . . . , n 1,ℓ + p1 } × {1, n 2,ℓ + p2 }) . Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

7

Throughout the paper, we focus only on homogeneous Dirichlet boundary conditions for notational simplicity and different kinds of boundary conditions can also be treated in the similar fashion. Other possible collocation points can be selected by using the Demko abscissae [17] or Galerkin superconvergent points [29]. Numerical performance in Section 5 will show that our proposed preconditioners work equally well for various choices of collocation points. Remark 3.4. The approximation properties of IGA collocation problem (16) have been studied, only for the onedimensional case, in [1,3]. This theory cannot be extended to higher spatial dimensions. Nevertheless, extensive numerical experience from the references mentioned in the introduction illustrates that IGA collocation method is stable and convergent in various problems of practical interest. 4. Multilevel preconditioners for IGA collocation methods The original problem (16) is associated to the following discrete system on the isogeometric discrete space V: Au = b.

(17)

For each i ∈ {0, 1, . . . , J }, let Ii be the matrix form of the interpolation operator (the natural inclusion) Ii from the space Vi to the finest space V = V J and Iit the transpose of Ii . We define the stiffness matrix at each level by Ai = Iit A Ii .

(18)

Then the matrix form of the BPX preconditioner becomes T = BA, where B is the BPX preconditioner t B = I0 A−1 0 I0 +

J ∑

Ii Ri Ii ,

(19)

i=1

where Ri , 1 ≤ i ≤ J , is either Di−1 or (Di − Ui )−1 Di (Di − Li )−1 with Di , Li , Ui the diagonal matrix, the lower and upper parts of Ai . We remark that Ri = Di−1 or Ri = (Di − Ui )−1 Di (Di − Li )−1 are Jacobi smoothers and symmetric Gauss–Seidel smoothers, respectively. The use of the B preconditioner (19) for the iterative solution of the discrete system Au = b can also be considered as the preconditioned system Tu = g,

with g = Bb,

which is accelerated by GMRES iteration. Remark 4.1. Although the numerical tests presented in the next section are implemented on a single-core machine, our BPX preconditioners (as additive multilevel preconditioners) are well suited for a parallel implementation on multi-core processors. We can derive the algebraic expression of IGA collocation multigrid V-cycle for solving (17): v0 = uk i = 0, 1, . . . , J,

v(i+1)/(2J +2) = vi/(2J +2) + Ii Ri Iit (b − Avi/(2J +1) ), (J +i+2)/(2J +2)

v

(J +i+1)/(2J +2)

=v

+

I J −i R J −i ItJ −i (b

(J +i+1)/(2J +2)

− Av

),

i = 0, 1, . . . , J,

(20)

uk+1 = v1 . 5. Numerical results In the present section, we test the optimality of the isogeometric collocation BPX preconditioners (19) and the uniform convergence of IGA collocation multigrid V-cycle (20) for model elliptic problems (1) on both parametric domain (reference square or cube) and physical domain (quarter-ring) with Dirichlet boundary conditions in 2D or ˆ on the 3D. The model problem is discretized by isogeometric NURBS spaces V on the deformed domain (or V −ℓ unit domain) with associated levels J relative to the grid sizes h 0 = 1, h ℓ = 2 , 1 ≤ ℓ ≤ J , polynomial degree p and regularity k, using a starting basis in the Matlab isogeometric library GeoPDEs [15,38]. We will consider Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

8

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

Table 1 Unpreconditioned (NP) and BPX preconditioners with symmetric Gauss–Seidel (GS) smoothers and Jacobi (Ja) smoothers (first three parts); Multigrid V-cycle (MG V-cycle) and GMRES preconditioned with V-cycle (MG PGMRES) with GS smoothers and Ja smoothers at each fine level (last three parts); On both the unit square domain and the deformed domain: GMRES iteration counts as a function of the number of levels. Fixed p = 3, k = 2. BPX PGMRES, p = 3, k = 2 Levels

iter

2

5

6

GS

Ja NP

GS

Ja

NP

GS

Ja NP

GS

Ja

NP

GS

Ja

19

7

10 36

9

12

70

10

12 150

11

12 297

11

12

2

3

On the square domain, Demko abscissae 4

5

6

NP

GS

Ja NP

GS

Ja

NP

GS

Ja NP

GS

Ja

NP

GS

Ja

18

7

10 35

9

11

68

10

12 136

11

12 327

11

12

Levels

iter

On the square domain, Greville abscissae 4

NP

Levels

iter

3

On the quarter of ring domain, Greville abscissae 3 4

2 NP

GS

31

7

Ja NP 9

59

GS 8

Ja

NP

12 124

5

GS

Ja NP

GS

9

13 358

9

6 Ja

NP

12 895

GS

Ja

9

12

MG PGMRES / MG V-cycle, p = 3, k = 2 Levels

2 NP

iter

NP

MG V-cycle

MG V-cycle

MG PGMRES

MG PGMRES

4 (GS) 5 (Ja)

4 (GS) 5 (Ja)

4 (GS) 5 (Ja)

4 (GS) 5 (Ja)

4 (GS) 5 (Ja)

On the square domain, Demko abscissae 3 4

2

Levels

MG V-cycle

NP

7 (GS)

35

31

MG V-cycle

NP

8 (GS)

68

MG V-cycle

NP

8 (GS)

136

297 8 (GS) 21 (Ja)

5

6

MG V-cycle

NP

8 (GS)

327

MG V-cycle 8 (GS)

MG PGMRES

MG PGMRES

MG PGMRES

MG PGMRES

MG PGMRES

4 (GS)

4 (GS)

4 (GS)

4 (GS)

4 (GS)

On the quarter of ring domain, Greville abscissae 3 4

2 NP

150 8 (GS) 21 (Ja)

6 NP

MG PGMRES

iter

8 (GS) 21 (Ja)

5

NP

MG PGMRES

18

70

MG V-cycle

36 8 (GS) 21 (Ja)

NP

iter

MG V-cycle

MG PGMRES

Levels

iter

NP

On the square domain, Greville abscissae 4

19 7 (GS) 20 (Ja)

iter

iter

MG V-cycle

3

MG V-cycle

NP

11 (GS)

59

MG V-cycle

NP

13 (GS)

124

MG V-cycle

NP

13 (GS)

358

5

6

MG V-cycle

NP

12 (GS)

895

MG V-cycle 13 (GS)

MG PGMRES

MG PGMRES

MG PGMRES

MG PGMRES

MG PGMRES

5 (GS)

6 (GS)

5 (GS)

5 (GS)

4 (GS)

as collocation points mainly the tensor-product Greville abscissae (15), and moreover in some tests the Demko abscissae [17] will be employed in order to show the good preconditioner performance with respect to various choices of collocation points. 5.1. BPX preconditioner with GMRES (BPX PGMRES) The IGA collocation problems (16) are solved by the GMRES method with the isogeometric collocation BPX preconditioner (19), with a zero initial guess and as stopping criterion a 10−6 reduction of the relative residual. In the following tests, we study how the optimality of the IGA collocation BPX preconditioner depends on J, p, k and on domain deformations. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

9

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx Table 2 Unpreconditioned (NP) and BPX preconditioners with Jacobi (Ja) smoothers, on the cubic domain: GMRES iteration counts as a function of the number of levels. Fixed p = 3, k = 2. Greville abscissae are considered as collocation points. BPX PGMRES, 3D cubic domain Levels

2

3

4

NP

Ja

NP

Ja

NP

Ja

iter

18

11

35

12

70

11

5.1.1. BPX optimality in levels J with mesh size h J In Table 1, we report the GMRES iteration counts of BPX preconditioners with Jacobi and symmetric Gauss–Seidel smoothers on the unit square with fixed p = 3, k = 2 and the tensor-product Greville abscissae as collocations points, which are similar to those obtained using as collocation points the Demko abscissae [17]. This example indicates that our BPX preconditioners are robust with respect to the choice of collocation points. For this reason, in what follows the Greville abscissae are chosen as collocation points unless otherwise stated. And in Table 1 the optimal convergence of BPX preconditioners is presented on the ring-shaped domain for fixed p = 3, k = 2. Those results are qualitatively similar to those on the square domain. Moreover, these tests show the optimal complexity for BPX preconditioners with Jacobi or symmetric Gauss–Seidel smoothers since iteration counts seem to stay constant independent of the mesh size h J (i.e., the number of levels J ). We also remark that although the iteration count is lower using symmetric Gauss–Seidel as smoothers, the iteration cost in our implementation is higher, thus in our tests the Jacobi smoothers perform better than symmetric Gauss–Seidel ones. In the following BPX tests, we will consider only BPX preconditioners with Jacobi smoothers. The GMRES iteration counts of BPX preconditioners are reported in Table 2 for the reference cubic domain in three spatial dimension. These 3D results confirm the h-independence and the optimality of BPX preconditioners, since the iteration counts are bounded above by a constant independent of J (we could not increase further the number of levels due to memory constraints). We notice by these results for varying J that the proposed IGA collocation BPX preconditioners perform as well as the analogous BPX preconditioners established for IGA Galerkin discretizations in [9]. 5.1.2. BPX dependence on p and k The aim of this test is to investigate how BPX preconditioners work for the polynomial degree p and spline regularity k, compared to the unpreconditioned case. We recall from Section 3.3 that even IGA Galerkin theory does not cover the dependence of BPX preconditioners on p and k. In order to investigate BPX preconditioners for higher values of p, we consider in Table 3 a discrete problem on the ring-shaped domain, with varying p from 3 to 7 with maximum regularity and various levels from 2 to 5. The GMRES iteration counts of the unpreconditioned case grow for increasing p, while our BPX preconditioners are optimal with respect to p with bounded iteration counts for fixed level. Also the GMRES iteration numbers on the unit square are reported in Table 4 where the spline degree p varies from 3 to 6 and the regularity k varies from 2 to 5. The GMRES iteration counts for the unpreconditioned method increase when p grows while they improve when the spline regularity k increases. The iteration count of BPX preconditioner seems to have a good behavior with respect to k for fixed p while for fixed k the iteration counts remain bounded when increasing p. 5.1.3. Robustness of BPX preconditioners with respect to domain deformations The four domains depicted in Fig. 1 are taken into account, deforming the rectangle [0, 0.5] × [0, 2] (A) into increasingly more curved boomerang-shaped domains (B, C, D). All the domains are discretized by IGA NURBS spaces associated with 4 level and 5 level (J = 4, 5), fixed p = 3, k = 2 where both Greville and Demko collocation points are considered. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

10

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

Table 3 Unpreconditioned (NP) and BPX preconditioners with Jacobi (Ja) smoothers (upper row); Multigrid V-cycle and GMRES preconditioned with V-cycle (lower row) with GS smoothers at each fine level; On the quarter of ring domain: GMRES iteration counts as a function of the number of levels and the spline polynomial degree p. Fixed k = p − 1. On the quarter of ring domain, Greville abscissae BPX PGMRES Levels

2 3 4 5

p=3

p=4

p=5

p=6

p=7

NP

Ja

NP

Ja

NP

Ja

NP

Ja

NP

Ja

31 59 124 358

9 12 13 12

37 66 165 459

12 13 14 15

43 73 214 593

12 15 16 16

59 81 265 726

12 16 18 19

80 98 295 860

13 17 19 19

MG PGMRES/MG V-cycle Levels

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

2 3 4 5

5 6 5 5

11 13 13 12

5 5 5 4

10 11 12 12

6 5 5 4

13 12 11 11

7 6 5 4

19 17 15 14

9 8 6 5

29 24 22 19

Table 4 Unpreconditioned (NP) and BPX preconditioners with Jacobi (Ja) smoothers (upper row); Iteration numbers of Multigrid V-cycle and GMRES preconditioned with V-cycle (lower row) with GS smoothers at each fine level; On the unit square domain: GMRES iteration counts as a function of the spline polynomial degree p and the regularity k. 5 level is fixed. On the square domain, fixed 5 level, Greville BPX PGMRES p=3

k k k k

=2 =3 =4 =5

p=4

p=5

p=6

NP

Ja

NP

Ja

NP

Ja

NP

Ja

150 – – –

12 – – –

211 205 – –

14 15 – –

245 233 229 –

15 15 17 –

348 332 308 298

17 17 17 19

MG PGMRES/MG V-cycle

k k k k

=2 =3 =4 =5

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

MG PGMRES

MG V-cycle

4 – – –

8 – – –

5 4 – –

10 8 – –

8 6 5 –

28 15 10 –

17 11 10 7

96 50 25 13

Table 5 shows that the GMRES iteration counts increase (going from A to D) by a factor 2.6 (4 level)/ 3.2 (5 level) for the unpreconditioned solver, while the iteration counts of BPX preconditioners appear to be independent of the domain deformation. We can state that BPX preconditioners are relatively less sensitive to domain deformations than the unpreconditioned solver and also expect that BPX preconditioners are very effective for the problems on more deformed CAD domain. 5.2. Multigrid V-cycle (MG V-cycle) and MG V-cycle preconditioner with GMRES (MG PGMRES) The numerical experiments are performed using stopping criterion

∥rk ∥2 ∥r0 ∥2

< 10−6 , where r0 is the initial residual

vector and rk is the residual of kth iteration. For the multigrid V-cycle (MG V-cycle), we will use Jacobi or symmetric Gauss–Seidel as the pre-smoother and the post-smoother on each fine level and use the exact solver on the coarsest level. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

11

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

Fig. 1. Domains for the 2D boomerang test.

Table 5 2D boomerang test: Iteration counts of unpreconditioned (NP) GMRES and BPX preconditioners with Jacobi (Ja) smoothers (upper row); Iteration counts of unpreconditioned GMRES (NP), multigrid V-cycle (MG V-cycle) and GMRES preconditioned with V-cycle (MG PGMRES) with GS smoothers at each fine level (lower row); On the curved boomerang-shaped domains. Fixed p = 3, k = 2. Both Greville and Demko abscissae are considered as collocation points. On the curved boomerang-shaped domains, p=3, k=2 Domain

BPX PGMRES, Greville 4 levels

BPX PGMRES, Demko 5 levels

4 levels

5 levels

NP

Ja

NP

Ja

NP

Ja

NP

Ja

A B C D

131 148 209 343

15 15 15 15

261 335 461 830

16 16 16 16

126 143 198 335

15 16 16 15

250 320 449 807

16 16 17 16

Domain

MG PGMRES/MG V-cycle, Greville abscissae 4 levels

A B C D

5 levels

NP

MG PGMRES

MG V-cycle

NP

MG PGMRES

MG V-cycle

131 148 209 343

8 9 12 17

20 26 42 90

261 335 461 830

8 9 13 18

20 28 53 117

5.2.1. MG V-cycle/MG PGMRES optimality in levels J with mesh size h J Table 1 displays the comparison of different choices of collocation points: Greville abscissae and Demko ones. The results thus lead to robustness in terms of the selection of collocation points. Also these tests show that our multigrid V-cycle methods are robust in terms of the choice of smoothers (e.g., Jacobi or symmetric Gauss–Seidel ones) at each fine level. Thus in the following multigrid V-cycle tests, a symmetric Gauss–Seidel smoother is adopted as a smoother at each fine level. Moreover, to show the effect of refinement levels on the convergence of MG V-cycle, the iteration counts on both the unit square and the deformed domain are reported in Table 1. It is not surprising to see that almost all iteration counts are same for all refinement levels. Although only the numerical results for the V-cycle methods are listed, we can expect that the W-cycle methods have similar convergence behavior, as was pointed out in Remark 3.3. We will also solve the preconditioned system (5) by using the following V-cycle multigrid as a preconditioner in GMRES method: One step of the standard V-cycle multigrid B : V → V is recursively defined as follows: Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

12

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

Let B0 = A−1 0 , for i > 0 and g ∈ Vi , define Bi g = w3 : (i) Presmoothing : w1 = Ri g; (ii) Correction: w2 = w1 + Bi−1 Q i−1 (g − Ai w1 ); (iii) Postsmoothing: w3 = w2 + Ri (g − Ai w2 ): Set B = B J . As pointed out earlier in Remark 3.2, in order to show the efficiency of GMRES preconditioned with V-cycle, the numerical results for MG PGMRES are given in Table 1. From these tests, we can see that both MG V-cycle and MG PGMRES converge uniformly with respect to the number of levels. 5.2.2. MG V-cycle/MG PGMRES dependence on p and k First we keep in mind that even IGA Galerkin theory does not cover the dependence of V-cycle methods on p and k. As a point of comparison, the GMRES iteration counts of MG V-cycle and MG PGMRES for various p and k are summarized in Tables 3 and 4. Table 3 shows that the GRMES iteration counts of MG PGMRES algorithm are optimal with respect to p with bounded iteration numbers for fixed level, while those of MG V-cycle algorithm grow for increasing p. It is noted that in Table 3, the faster convergence on higher levels can be explained by the fact that the multigrid solution gets closer to the IGA collocation solution as the mesh gets finer. Table 4 reports the GMRES iteration counts, varying p from 3 to 6, k from 2 to 5. The iteration counts of both MG V-cycle methods show a good performance with respect to k for fixed p, while for fixed k they grow when increasing p. In addition, the iteration numbers appear to improve when the regularity k increases. 5.2.3. Robustness of MG V-cycle/MG PGMRES with respect to domain deformations Finally we study how MG V-cycle and MG PGMRES perform in the presence of domain deformations drawn in Fig. 1. Both Greville and Demko collocation points are considered. The results presented in Table 5 show that the V-cycle methods greatly reduce the ill-conditioning of the unpreconditioned problem. The results imply that MG PGMRES is less sensitive to domain deformations than the unpreconditioned solver since the GMRES iteration counts increase (going from A to D) by a factor 2.6 (4 level)/ 3.18(5 level) for the unpreconditioned solver, while only by a factor 2.1 (4 level)/ 2.25 (5 level) for MG PGMRES algorithm. 5.3. Comparison of BPX preconditioners with multigrid V-cycle methods Tables 1 and 2 show that BPX preconditioners and two multigrid V-cycle methods are all very robust with respect to the choice of collocation points as well as with respect to the selection of smoothers and also show the effective iteration counts versus the number of levels of our proposed multilevel preconditioners under h-refinements. From Tables 3 and 4, we observe that for various p and k, BPX preconditioners and two multigridbased preconditioners have much better convergence behavior than the unpreconditioned case. Moreover, Table 5 shows that the preconditioners that we have proposed are less sensitive with respect to domain deformation than the unpreconditioned solver and that comparing with BPX case, the domain deformation has a stronger effect on both multigrid preconditioners. The theory for the stronger effect is beyond the scope of the present paper and will be addressed as part of future work. We would like to remark from the numerical experiments that multigrid V-cycle algorithms outperform BPX preconditioners in almost all cases. 6. Conclusions We have presented optimal multilevel preconditioners for IGA collocation methods. The optimality of the preconditioners under h-refinement has been obtained by several numerical experiments in 2D and 3D, showing the practical performance of the proposed preconditioners. We expect that our preconditioners can be extended to elliptic systems such as compressible elasticity, while p-robust multilevel preconditioners will need further research. Acknowledgments The work of Durkbin Cho was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (2015R1A1A1A05001109) and by the Ministry of Education, South Korea (2018R1D1A1B07048773). This support is gratefully acknowledged. The author is grateful to anonymous referee for his/her comments and suggestions. Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

13

References [1] C. Anitescu, Y. Jia, Y.J. Zhang, T. Rabczuk, An isogeometric collocation method using superconvergent points, Comput. Methods Appl. Mech. Engrg. 284 (2015) 1073–1097. [2] F. Auricchio, F. Calabro, T.J.R. Hughes, A. Reali, G. Sangalli, A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 15–27. [3] F. Auricchio, L. Beirão da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli, Isogeometric collocation methods, Math. Models Methods Appl. Sci. 20 (11) (2010) 2075–2107. [4] F. Auricchio, L. Beirão da Veiga, T.J.R. Hughes, A. Reali, G. Sangalli, Isogeometric collocation for elastostatics and explicit dynamics, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 2–14. [5] F. Auricchio, L. Beirão da Veiga, J. Kiendl, C. Lovadina, A. Reali, Locking-free isogeometric collocation methods for spatial Timoshenko rods, Comput. Methods Appl. Mech. Engrg. 263 (2013) 113–126. [6] Y. Bazilevs, L. Beirão da Veiga, J.A. Cottrell, T.J.R. Hughes, G. Sangalli, Isogeometric analysis: approximation, stability and error estimates for h-refined meshes, Math. Models Methods Appl. Sci. 16 (2006) 1–60. [7] C. de Boor, A Practical Guide to Splines, Springer, 2001. [8] J.H. Bramble, J.E. Pasciak, J. Xu, Parallel multilevel preconditioners, Math. Comp. 55 (1990) 122. [9] A. Buffa, H. Harbrecht, A. Kunoth, G. Sangalli, BPX-preconditioning for isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 265 (2013) 63–70. [10] L. Chen, R. Nochetto, J. Xu, Optimal multilevel methods for graded bisection grids, Numer. Math. 120 (2012) 1–34. [11] D. Cho, R. Vazquez, BPX preconditioners for isogeometric analysis using analysis suitable T-splines, IMA J. Numer. Anal. (2019) http://dx.doi.org/10.1093/imanum/dry032, dry032. [12] N. Collier, L. Dalcin, D. Pardo, V.M. Calo, The cost of continuity: Performance of iterative solvers on isogeometric finite elements, SIAM J. Sci. Comput. 35 (2013) 767–784. [13] N. Collier, D. Pardo, L. Dalcin, M. Paszynsk, V.M. Calo, The cost of continuity: a study of the performance of isogeometric finite elements using direct solvers, Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 353–361. [14] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis. Towards Integration of CAD and FEA, Wiley, 2009. [15] C. De Falco, A. Reali, R. Vazquez, GeoPDEs: A research tool for isogeometric analysis of PDEs, Adv. Eng. Softw. 42 (2011) 1020–1034, 2010. [16] L. De Lorenzis, J.A. Evans, T.J.R. Hughes, A. Reali, Isogeometric collocation: Neumann boundary conditions and contact, Comput. Methods Appl. Mech. Engrg. 284 (2015) 21–54. [17] S. Demko, On the existence of interpolation projectors onto spline spaces, J. Approx. Theory 43 (1985) 151–156. [18] M. Donatelli, C. Garoni, C. Manni, S. Serra-Capizzano, H. Speleers, Robust and optimal multi-iterative techniques for IgA Galerkin linear systems, Comput. Methods Appl. Mech. Engrg. 284 (2015) 230–264. [19] M. Dryja, O.B. Widlund, Domain decomposition algorithms with small overlap, SIAM J. Sci. Comput. 15 (3) (1994) 604–620. [20] G.E. Farin, NURBS Curves and Surfaces: From Projective Geometry to Practical Use, A.K. Peters, 1995. [21] K.P.S. Gahalaut, J.K. Kraus, S.K. Tomar, Multigrid methods for isogeometric discretization, Comput. Methods Appl. Mech. Engrg. 253 (2013) 413–425. [22] K.P.S. Gahalaut, S.K. Tomar, Condition number estimates for matrices arising in the isogeometric discretizations. Tech. Report RICAM, 23, 2012. [23] K.P.S. Gahalaut, S.K. Tomar, J.K. Kraus, Algebraic multilevel preconditioning in isogeometric analysis: Construction and numerical studies , Comput. Methods Appl. Mech. Engrg. 266 (2013) 40–56. [24] H. Gomez, A. Reali, G. Sangalli, Accurate, efficient, and (iso)geometrically flexible collocation methods for phase-field models, J. Comput. Phys. 262 (2014) 153–171. [25] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD,finite elements NURBS exact geometry mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [26] T.J.R. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 199 (2010) 301–313. [27] R. Kettler, Analysis and comparison of relaxation schemes in robust multigrid and preconditioned conjugate gradient methods, in: Multigrid Methods (Cologne, 1981), in: Lecture Notes in Math., vol. 960, Springer, 1982, pp. 502—534. [28] S.K. Kleiss, C. Pechstein, B. Jüttler, S. Tomar, IETI-isogeometric tearing and interconnecting, Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 201–215. [29] M. Montardini, G. Sangalli, L. Tamellini, Optimal-order isogeometric collocation at Galerkin superconvergent points, Comput. Methods Appl. Mech. Engrg. 316 (2017) 741–757. [30] L.F. Pavarino, S. Scacchi, Isogeometric block FETI-DP preconditioners for the Stokes and mixed linear elasticity systems, Comput. Methods Appl. Mech. Engrg. 310 (2016) 694–710. [31] L. Piegl, W. Tiller, The NURBS Book, 2nd Ed., Springer-Verlag, 1997. [32] A. Reali, T.J.R. Hughes, An introduction to isogeometric collocation methods, in: Isogeometric Methods for Numerical Simulation, in: CISM International Centre for Mechanical Sciences, vol. 561, Springer, Vienna, 2015, pp. 173–204. [33] D.F. Rogers, An Introduction to NURBS with Historical Perspective, Academic Press, 2001. [34] D. Schillinger, J.A. Evans, A. Reali, M.A. Scott, T.J.R. Hughes, Isogeometric collocation: cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations, Comput. Methods Appl. Mech. Engrg. 267 (2013) 170–232. [35] B.F. Smith, P. Bjørstad, W.D. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, 1996.

Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.

14

D. Cho / Mathematics and Computers in Simulation xxx (xxxx) xxx

[36] O. Tatebe, Y. Oyanagi, Efficient implementation of the multigrid preconditioned conjugate gradient method on distributed memory machines, Supercomputing (1994) 194–203. [37] A. Toselli, O.B. Widlund, Domain decomposition methods: Algorithms and theory, in: Computational Mathematics, Vol. 34, Springer-Verlag, Berlin, 2004. [38] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Comput. Math. Appl. 72 (3) (2016) 523–554. [39] L. Beirão da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, Overlapping Schwarz methods for isogeometric analysis, SIAM J. Numer. Anal. 50 (2012) 1394–1416. [40] L. Beirão da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, BDDC preconditioners for isogeometric analysis, Math. Models Methods Appl. Sci. 23 (2013) 1099–1142. [41] L. Beirão da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, Isogeometric Schwarz preconditioners for linear elasticity systems, Comput. Methods Appl. Mech. Engrg. 253 (2013) 439–454. [42] L. Beirão da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, Overlapping Schwarz preconditioners for isogeometric collocation methods, Comput. Methods Appl. Mech. Engrg. 278 (2014) 239–253. [43] L. Beirão da Veiga, C. Lovadina, A. Reali, Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods, Comput. Methods Appl. Mech. Engrg. 241–244 (2012) 38–51. [44] L. Beirão da Veiga, L.F. Pavarino, S. Scacchi, O.B. Widlund, S. Zampini, Isogeometric BDDC preconditioners with deluxe scaling, SIAM J. Sci. Comput. 36 (3) (2014) A1118–A1139. [45] L. Beirão da Veiga, L.F. Pavarino, S. Scacchi, O.B. Widlund, S. Zampini, Adaptive selection of primal constraints for isogeometric BDDC deluxe preconditioners, SIAM J. Sci. Comput. 39 (1) (2017) A281–A302. [46] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Rev. 34 (4) (1992) 581–613.

Please cite this article as: D. Cho, Optimal multilevel preconditioners for isogeometric collocation methods, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.08.003.