Isogeometric spectral approximation for elliptic differential operators

G Model

ARTICLE IN PRESS

JOCS-879; No. of Pages 7

Journal of Computational Science xxx (2018) xxx–xxx

Contents lists available at ScienceDirect

Journal of Computational Science journal homepage: www.elsevier.com/locate/jocs

Isogeometric spectral approximation for elliptic differential operators Quanling Deng a,b,∗ , Vladimir Puzyrev a,b , Victor Calo a,b,c a

Curtin Institute for Computation, Curtin University, Kent Street, Bentley, Perth, WA 6102, Australia Department of Applied Geology, Curtin University, Kent Street, Bentley, Perth, WA 6102, Australia c Mineral Resources, Commonwealth Scientific and Industrial Research Organisation (CSIRO), Kensington, Perth, WA 6152, Australia b

a r t i c l e

i n f o

Article history: Received 2 February 2018 Received in revised form 11 May 2018 Accepted 14 May 2018 Available online xxx Keywords: Differential operator Spectral approximation Isogeometric analysis Optimally-blended quadratures Schrödinger operator

a b s t r a c t We study the spectral approximation of a second-order elliptic differential eigenvalue problem that arises from structural vibration problems using isogeometric analysis. In this paper, we generalize recent work in this direction. We present optimally-blended quadrature rules for the isogeometric spectral approximation of a diffusion-reaction operator with both Dirichlet and Neumann boundary conditions. The blended rules improve the accuracy of the isogeometric approximation. In particular, the optimal blending rules minimize the dispersion error and lead to two extra orders of super-convergence in the eigenvalue error. Various numerical examples (including the Schrödinger operator for quantum mechanics) in one and three spatial dimensions demonstrate the performance of the blended rules. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Differential eigenvalue problems arise in a wide range of applications, such as the vibration of elastic bodies in structural mechanics and the multi-group diffusion in nuclear reactors [42]. In general, analytical solutions for these problems are impossible and numerical methods are used. Numerical methods for approximating these differential eigenvalue problems lead to a generalized matrix eigenvalue problem, which is then solved numerically. Different numerical methods result in different matrices and the widely-used methods include finite elements [3,11,16,17,34,35,37,42], isogeometric elements [10,13,15,24,25,32,39,40], discontinuous Galerkin (DG) [2,29], hybridizable discontinuous Galerkin (HDG) [30], and a recently developed hybrid high-order (HHO) method [14]. Early work [11,37,42] used conforming finite elements on simplicial meshes and the method demonstrated convergence rates of order h2p for the eigenvalues and of order hp in the energy norm for the eigenfunctions provided that the eigenfunctions are smooth enough in the sense of function regularity (mainly derivative continuities, see [18]). Herein, h is the mesh size. For pth order finite element approximations, the convergence rate of 2p of eigenvalue

∗ Corresponding author at: Curtin Institute for Computation, Curtin University, Kent Street, Bentley, Perth, WA 6102, Australia. E-mail addresses: [email protected] (Q. Deng), [email protected] (V. Puzyrev), [email protected] (V. Calo).

2p

error means that the error decreases to (roughly) h when the mesh size h is decreased to h . Any rate greater than 2p is referred to as a superconvergent rate (see, for example, [18,26,42]). Sharp and optimal estimates of the numerical eigenfunctions and eigenvalues of finite element analysis are established in [4–6]. In [16,34,35], the authors studied the spectral approximation of elliptic operators by mixed and mixed-hybrid methods and optimal error estimates were established. Similar results were obtained more recently in [2,14,29,30] using the DG, HDG, and HHO methods, which are numerical methods using discontinuous basis functions. The spectral approximation by the HDG method leads to a convergence of order h2p+1 for the eigenvalues. A non-trivial post-processing, which utilizes a Rayleigh quotient, is also examined numerically in [30] which leads to an improved convergence of order h2p+2 for p ≥ 1. The HHO approximation delivers a convergence of order h2p+2 for the eigenvalue errors for all polynomial degrees (p ≥ 0). Isogeometric analysis is a numerical method introduced in 2005 [20,31]. The spectral approximation of the elliptic operators arising in structural vibrations were investigated using isogeometric analysis in [21,33] and the method shows improved spectral approximations over the classical finite elements [20]. In [33], a duality principle, which induces a bijective map from spectral analysis to dispersion analysis, was established, which unifies the spectral analysis for structural vibrations (eigenvalue problems) and the dispersion analysis for wave propagations. Further advantages of the method on spectral approximation properties are investigated in [32].

https://doi.org/10.1016/j.jocs.2018.05.009 1877-7503/© 2018 Elsevier B.V. All rights reserved.

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009

G Model

ARTICLE IN PRESS

JOCS-879; No. of Pages 7

Q. Deng et al. / Journal of Computational Science xxx (2018) xxx–xxx

2

The recent work in [15,39] studies both theoretically and numerically the optimally blended quadrature rules [1] for the isogeometric analysis of the Laplace eigenvalue problem. The blended quadrature is a combination of two classical quadratures (see Section 3 below). In [15], the authors establish for p = 1, 2, . . ., 7 the super-convergence of order h2p+2 for the eigenvalue errors while maintaining optimal convergence of orders hp and hp+1 for the eigenfunction errors in the H1 -seminorm and in the L2 -norm, respectively. The work [25] introduces the dispersion-minimized mass for isogeometric analysis and generalizes the results to arbitrary B-spline of degree p. In [40], the authors study the optimally blended quadratures for isogeometric analysis with variable continuity. To reduce the computational costs, [24] describes new quadrature rules to replace the optimal blending rules. For the source problems, optimal (Gaussian) quadrature rules were proposed for isogeometric analysis in [7–9]. In this work, we generalize the work in [15,39] to include the reaction effects in the differential operator as well as to consider different boundary conditions. We study the optimal blending quadratures numerically for the generalized differential operator with both Dirichlet and Neumann boundary conditions. We apply the blending rules to approximate the spectrum of the Schrödinger operator. The outline of the remainder of this paper is as follows. We first describe the model problem and the isogeometric discretization in Section 2. We introduce classical and blended quadrature rules in Section 3. A brief dispersion error estimations is presented in this section. Numerical examples are given in Section 4. Finally, Section 5 summarizes our findings and describes future research directions. 2. Problem statement We consider the second-order differential eigenvalue problem: Find the eigenpair (, u) such that −u + u

= u in ,

u

= 0 on ∂,

(2.2)

with an associated set of orthonormal eigenfunctions uj



uj (x)uk (x) dx = ıjk ,

(2.3)



where ıjk is the Kronecker delta. The set of all the eigenvalues is the spectrum of the operator. We normalize the eigenfunctions in the L2 space and hence the eigenfunctions are orthonormal with each other under the scalar inner product. Now, let us define two bilinear forms



∇ w · ∇ v + wv dx and b(w, v) 



= (w, v) =

wv dx,

∀w, v ∈ H01 (),

(2.4)



where H01 () is the Sobolev space with functions vanishing at the boundary ∂. These two inner products are associated with the following energy and L2 norms w E =



a(w, w),

w = w L2 () =



(w, w).

(2.5)

(2.6)

where we have used the integration by parts on (2.1) and selected the weighting functions to be eigenfunctions. Let H01 () be the standard Hilbert space where both functions and their weak derivatives are in L2 () with functions vanishing ¯ = [0, 1]d into a uniat the boundary. We discretize the domain  form tensor-product grid with the grid nodes located at (h1 n1 , . . ., hd nd ), where hk = N1 , k = 1, . . ., d, is the size of the kth dimenk

sion and nk = 0, 1, . . ., Nk and Nk is the number of intervals in the kth dimension. We denote each element as K and their collection ¯ = ∪K ∈ T K. Due to the tensor-product structure of as Th such that  h the discretization, the element size is h =

 d

h2 . At the contink=1 k

uous level, the weak formulation for the eigenvalue problem (2.1) is: Find all eigenvalues  ∈ R and eigenfunctions u ∈ V ⊂ H01 () such that, a(w, u) = b(w, u),

∀ w ∈ H01 (),

(2.7)

while at the discrete level, the isogeometric analysis for the eigenvalue problem (2.1) is: Find all eigenvalues h ∈ R and eigenfunctions uh ∈ Vh such that, a(wh , uh ) = h b(wh , uh ),

∀ wh ∈ Vh (),

(2.8)

where Vh ⊂ H01 () is the solution and test space, which is spanned by the B-spline or non-uniform rational basis spline (NURBS) basis functions. In this paper, we focus on B-spline basis functions. Following [23,38], the definition of the B-spline basis functions with maximum continuity in one dimension is as follows. Let X = {x0 , x1 , . . ., xm } be a knot vector with knots xj , that is, a nondecreasing sequence of real numbers called knots. The jth B-spline j basis function of degree p, denoted as p (x), is defined as

 j

0 < 1 < 2 ≤ · · · ≤ j ≤ · · ·

a(w, v) =

a(uj , uk ) = (−uj + uj , uk ) = (j uj , uk ) = j (uj , uk ) = j ıjk ,

(2.1)

where  = ∇ 2 is the Laplacian,  = (x) ∈ L2 () and is non-negative function, and  ⊂ Rd , d = 1, 2, 3 is a bounded open domain with Lipschitz boundary. This problem is a Sturm-Liouville eigenvalue problem (see, for example, [27,42]) which has a countable infinite set of eigenvalues j ∈ R

(uj , uk ) =

Using this notation, the eigenfunctions are also orthogonal with each other with respect to the energy inner product, i.e.,

0 (x) j p (x)

1,

if xj ≤ x < xj+1

= 0, otherwise x − xj xj+p+1 − x j j+1 =  (x) +  (x). xj+p − xj p−1 xj+p+1 − xj+1 p−1

(2.9)

In this paper, we use the B-splines on uniform meshes with non-repeating knots, that is, we use B-splines with maximum continuity. We approximate the eigenfunction as a linear combination of the B-spline basis functions. Using linearity and substituting all the B-spline basis functions for wh in (2.8) leads to the matrix eigenvalue problem KU = h MU,

(2.10)

where Kab = a(a , b ), Mab = b(a , b ), and U contains the coefficients of the eigenfunction uh , for its representation in the linear combination of the B-spline basis functions. Herein, we refer to [20] for the treatment on the boundary elements. For simplicity, the matrix K (although it contains a scaled mass) is referred as the stiffness matrix while the matrix M is referred as the mass matrix, and (h , uh ) is the unknown eigenpair. 3. Blending quadratures and dispersion errors In this section, we present the quadratures as well as their optimal blendings. Following earlier work [15,39] and its recent generalization [25], we omit the details to briefly give the dispersion errors for the quadrature rules. Herein, the dispersion error is the error of approximate frequency (characterized by eigenvalues in the following context; see (3.23)) to the exact one and we refer

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009

G Model

ARTICLE IN PRESS

JOCS-879; No. of Pages 7

Q. Deng et al. / Journal of Computational Science xxx (2018) xxx–xxx

to [1] for its formal definition. The blending rules are optimal in the sense of delivering minimal dispersion error. 3.1. Quadrature rules In computational applications, we evaluate the integrals involved in a(wh , uhj ) and b(wh , uhj ) numerically, that is, approximated by quadrature rules. On a reference element Kˆ = [0, 1]d , where d is the number of dimensions, a quadrature rule is of the form



fˆ (ˆx) dˆx = Kˆ

n 

Remark 3.1. For multidimensional problems on tensor-product grids, the stiffness and mass matrices can be expressed as Kronecker products of 1D matrices [28]. For example, in the 2D case, assume j that  is a constant. We define  =  2D = 2 1D and let pi (x)p (y) and pk (x)pl (y) be two 2D basis functions. Using the definition (2.4), we calculate j

Kijkl = a(pi (x)p (y), pk (x)pl (y))





ˆ l fˆ (nˆ l ) + Eˆ n ,



(3.11)

=

f (x) dx = K

l,K f (nl,K ) + En ,

(3.12)

l=1

ˆ l and nl,K = (nˆ l ). We refer to [41,26,12] for where l,K = det(JK ) details of how quadrature rules work on general elements. The quadrature rule is exact for a given function f(x) when the remainder En is exactly zero. For simplicity and notation purpose, we denote by Gm the m-point Gauss-Legendre quadrature rule, by Lm the m-point Gauss-Lobatto quadrature rule, and by Op the optimal blending scheme for the pth order isogeometric analysis with maximum continuity. In one dimension, Gm and Lm fully integrate polynomials of order 2m − 1 and 2m − 3, respectively (see, for example, [8,41]). We will use them in the following notations when referring to a particular quadrature rule. Applying the quadrature rules to (2.8), we have the approximated form ˜ h bh (wh , u˜ h ), ah (wh , u˜ h ) = 

∀ wh ∈ Vh ,

(3.13)

ah (w, v) =



n 1 

K ∈ Th n2 

+

l=1



pi (x)pk (x) dx

X

Y

(2) (2) (2) (2) (nl,K )w(nl,K )v(nl,K ) l,K

∂pj (y) ∂pl (y) ∂y

∂y

(3.17)

1D Herein, M1D ij and Kij are the mass and stiffness matrices of the 1D problem with  =  1D in (2.1). We refer the reader to [22] for the description of the summation rules.

3.2. Blended quadratures (1)

(2) (2) n2 { l,K , nl,K } , l=1

the blended quadrature rule, denoted as Q , is

Q = Q1 + (1 − )Q2 ,

(3.19)

where  is referred to as the blending parameter. A blended quadrature rule, which leads to minimized discrete dispersion error, is referred to as an optimally-blended quadrature rule (see [39,1]). For  = 0, 1, it reduces to a standard quadrature rule. Thus, we assume / 0, 1. Applying the blended rule Q for the intein this section  = gration of a function f, we have

n 1 



Q (1) (1) f (nl,K ) + En11 l,K

l=1

+(1 − )

n 2 

Q (2) (2) f (nl,K ) + En22 l,K

l=1



n2 

(j)

(3.15)

(j)

where { l,K , nl,K } with j = 1, 2, 3 specifies three (possibly different) quadrature rules. Here, we assume that we apply the same quadrature rules for the L2 inner products in (3.14) and (3.15). With these quadrature rules, we can rewrite (with slight abuse of the notation) the matrix eigenvalue problem (2.10) as ˜ h MU, ˜ = ˜ KU

(1) (1) f (nl,K ) + (1 − ) l,K

l=1

K ∈ Th l=1

(3.16)

˜ is the corresponding where Kab = ah (a , b ), Mab = bh (a , b ), and U representation of the eigenvector as the coefficients of the basis functions.

and Q2 =

defined as

= (3) (3) (3) w(nl,K )v(nl,K ), l,K

(1) n1 l=1

Given two quadrature rules Q1 = { l,K , nl,K }

n1

n3 

dy

(3.18)

(3.14)

and

 j

+ 1D p (y)pl (y)

1D Mijkl = M1D ik ⊗ Mjl .

l=1

bh (w, v) =

Y

where X and Y specify the intervals of each dimension in  and ⊗ refers to the tensor-product. Similarly, we obtain

K



j

p (y)pl (y) dy

dx

1D 1D 1D = K1D ik ⊗ Mjl + Mik ⊗ Kjl ,

f (x) dx =  (1) (1) ∇ w(n(1) ) · ∇ v(nl,K ) l,K l,K



+ 1D pi (x)pk (x)

∂x



+



where for w, v ∈ Vh



∂pi (x) ∂pk (x) ∂x

X

ˆ l are the weights, nˆ l are the coordinates of the quadrature where nodes, fˆ is the integrand f transformed into this reference element ˆ n is the number of quadrature points, and Eˆ n is the error of the K. quadrature rule. For each element K, we assume that there is an ˆ which determines the correinvertible map such that K = (K), ˆ Assuming JK is the spondence between the functions on K and K. Jacobian of the mapping, (3.11) induces a quadrature rule over the element K given by n 

∇ (pi (x)pj (y)) · ∇ (pk (x)pl (y)) + 2D pi (x)pj (y)pk (x)pl (y) dx

=

l=1



3



Q

Q

+ En11 + (1 − )En22

(2) (2) f (nl,K ) l,K

l=1

,

(3.20)

where EnQjj is the quadrature error of Qj for j = 1, 2, respectively. Thus, the error for the blending rule is the same as blending of the errors, that is, E Q = EnQ11 + (1 − )EnQ22 .

(3.21)

Assuming that Q1 and Q2 integrate polynomials up to order k1 and k2 , respectively, (3.21) shows that the blending rule integrates polynomials up to order min(k1 , k2 ). For example, in one dimension, the blending rule Q = Gm + (1 − )Lm

(3.22)

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009

G Model JOCS-879; No. of Pages 7

ARTICLE IN PRESS Q. Deng et al. / Journal of Computational Science xxx (2018) xxx–xxx

4

integrates polynomials up to order 2m − 3. For the dispersion analysis on the Helmholtz equation ( = 0,  = ω2 in (2.1)), the blended rule shows smaller dispersion errors. In fact, the optimal blending of spectral elements and finite elements, which is realized by optimally blended quadratures, leads to two extra order of super-convergence on the dispersion error; see [1]. This fact motivates the work (see [15,39,25]) of finding the optimal blending rules for the isogeometric analysis for differential eigenvalue problems (2.1) with  = 0. In the following section, we present the quadrature blending rules which minimize the dispersion errors of the discrete approximations. 3.3. Dispersion errors and optimal blending quadratures Following earlier work [15,25], based on the dual principle in [33], the dispersion errors of the isogeometric elements using quadratures can be characterized by the eigenvalue errors. For simplicity, we assume that  = 0. For C1 quadratic isogeometric elements (see [39] for classical linear finite elements), which use quadratic B-spline basis functions with C1 continuity (see [20] for details and for other higher order isogeometric elements), the relative errors are G

h 3 − 

=

 L

h3 − 

1

4 + O( 6 ), 720



=

2 − 3 4

+ O( 6 ). 1440

(3.24)

For  = 2/3, we obtain the two extra orders in the error representation and we call this case the optimal blending. The error representation of the optimal blending is O

h p −  

=

11

6 + O( 8 ). 60480

(3.25)

For C2 cubic elements, the optimal blending parameter is  =−3/2 and we refer to [15] for p ≤ 7 and [25] for the general case. The convergence rate for eigenpairs computed using isogeometric elements is O( 2p ) as shown in [21]. The optimal blending leads to a O( 2p+2 ) convergence rate for the relative eigenvalue errors. For a constant , we redefine the eigenvalue problem (2.1) as −u u

ˆ = u

in ,

= 0 on ∂,

(3.26)

ˆ =  − . Once the eigenvalue problem (3.26) is solved where  using isogeometric elements with optimal blending rules, we postprocess the approximated eigenvalue of (2.1) as ˆ h + . h = 

compared to the corresponding exact eigenvalues j . Figure 4.1 compares the approximation relative eigenvalue |− | error, which is defined as  h , of C1 quadratic isogeometric elements using the standard Gaussian quadrature and the optimal rule for problem (2.1) with homogeneous Neumann boundary conditions. The use of the optimal quadrature leads to more accurate results. The optimal ratio of blending of the Lobatto and Gauss quadrature rules in this case is 2:1 ( = 2/3), which in this particular case coincides with the ratio proposed by Ainsworth and Wajid [1] for finite-spectral elements of the same polynomial order. This ratio is different for higher order isogeometric elements [15]. Figure 4.2 shows the convergence of the errors in the eigenvalue approximation with C1 quadratic isogeometric elements. The optimal quadrature rule has two extra orders of convergence in the eigenvalue errors compared to the standard fully-integrated isogeometric elements. Not only the convergence rate, but also the errors are significantly lower for the optimal rule.

4.2. 3D results

where = ωh with ω2 =  being the exact eigenvalue and Q h denotes the approximate eigenvalue while using the quadrature rule Q with Q = G3 , L3 . The blending of these two rules, that is, Q = G3 + (1 − )L3 , leads to the error representation  Q − h

The 1D elliptic eigenvalue problem (2.1) with  = 0 and homogeneous Neumann boundary √ conditions (both ends) has the exact eigenpairs j = j2 2 , uj = 2 cos(jx), j = 1, 2, . . .. The approximate eigenvalues hj are sorted in ascending order and are

(3.23)

1 =−

4 + O( 6 ), 1440



4.1. 1D results

(3.27)

4. Numerical examples To show how the optimal quadrature rules reduce the errors in isogeometric analysis, we present one and three dimensional studies of equation (2.1) in this section, for the case in which  is both constant and variable.

Next, we continue our study with the dispersion properties of the three-dimensional eigenvalue problem (2.1) on tensor-product meshes. Optimal methods for multidimensional problems with constant coefficients and affine mappings can be formed by tensor product of the 1D mass and stiffness matrices (3.18). The exact eigenvalues and eigenfunctions of the 3D eigenvalue problem are given by klm = (k2 + l2 + m2 )2 ,

uklm = 2 sin(kx) sin(ly) sin(mz), (4.28)

for k, l, m = 1, 2, . . .. Figure 4.3 shows the relative eigenvalue errors with C1 quadratic isogeometric elements. Similar to the 1D case, the optimal scheme has two extra orders of convergence in the eigenvalue errors. Figure 4.4 compares the eigenvalue errors of the standard Gauss rule using C1 quadratic elements with the optimal scheme ( = 2/3). The latter has significantly better approximation properties in the entire domain. These results demonstrate that the use of optimal quadratures in isogeometric analysis significantly improves the accuracy of the discrete approximations compared to the fully-integrated Gaussbased method.

4.3. Spectral approximation of Schrödinger operator Following the analytical work on Schrödinger operators in [19], we study their numerical approximations in this subsection. We consider the 1D Schrödinger equation of a quantum particle trapped by the Pöschl-Teller potential



2

˛(˛ + 1) ˇ(ˇ + 1) + cos2 (y) sin2 (y)

 ,

0 < y <

 , ˛, ˇ > 0. 2

(4.29)

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009

G Model

ARTICLE IN PRESS

JOCS-879; No. of Pages 7

Q. Deng et al. / Journal of Computational Science xxx (2018) xxx–xxx

5

Fig. 4.1. Relative eigenvalue errors for C1 quadratic isogeometric elements with standard Gauss quadrature rule and optimal rule on linear (left) and logarithmic scales (right). The total number of degrees of freedom (discrete modes) is N = 1000 and I = 1, 2, . . ., N.

Fig. 4.2. Convergence of the errors in the eigenvalue approximation using C1 quadratic isogeometric elements with standard and optimal quadratures. The second (left), fourth (middle) and eighth (right) eigenvalues are shown.

Fig. 4.3. Convergence of the relative errors in the eigenvalue approximation using C1 quadratic isogeometric elements with standard and optimal quadratures. The second (left), tenth (middle) and sixteenth (right) eigenvalues are shown.

Applying the scaling x = y, the eigenvalue problem reads: Find the eigenpair (, u) such that

where we choose ˛ = ˇ = 1 for simplicity. This eigenvalue problem has the true eigenvalues (see for example [19]) 2

 = (4 + 2j) , −

d2 u + dx2



˛(˛ + 1) ˇ(ˇ + 1) + cos2 x sin2 x

 u

= u,

u(0) = 0,  u( ) = 0, 2

0
 , 2 (4.30)

j = 0, 1, 2, . . ..

(4.31)

Table 1 shows the relative eigenvalue errors for the first, second and fourth eigenmodes. We present the errors while using both the Gauss rule and optimally blended rule. Here, since the Pöschl-Teller potential blows up at the points x = 0, 2 and the Lobatto rules utilize the interval element end knots as quadrature points, we use the Gp+1 and Gp optimally blended rules (alternatively, one can use the

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009

G Model

ARTICLE IN PRESS

JOCS-879; No. of Pages 7

Q. Deng et al. / Journal of Computational Science xxx (2018) xxx–xxx

6

Fig. 4.4. Relative eigenvalue errors for C1 quadratic isogeometric elements with standard Gauss (left) and optimal quadrature rule (right). Color represents the absolute value of the relative error. Isosurfaces show 0.2% and 1.0% levels of the relative error. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 1 Relative eigenvalue (EV) errors for C0 linear and C1 quadratic isogeometric elements with Gauss rule Gp+1 and optimally-blended rules Qp . Set

|h1 − 1 |/1

|h2 − 2 |/2

|h4 − 4 |/4

p

N

Gp+1

Op

Gp+1

Op

Gp+1

Op

1

40 80 160

3.19e−3 7.41e−4 1.78e−4

6.60e−4 8.43e−5 1.06e−5

1.06e−2 2.49e−3 6.04e−4

1.65e−3 2.19e−4 2.80e−5

3.95e−2 9.33e−3 2.27e−3

3.81e−3 5.97e−4 8.07e−5

10 20 40

2.08 1.63e−3 7.94e−5 4.62e−6

2.98 2.65e−4 2.39e−6 1.11e−7

2.07 1.68e−2 6.68e−4 3.61e−5

2.94 4.29e−3 6.54e−5 5.24e−7

2.06 1.02e+0 9.07e−3 4.07e−4

2.78 2.73e−1 1.95e−3 2.83e−5

4.23

5.61

4.43

6.50

5.64

6.62

1 2 2

equivalent nonstandard quadratures; see [15,24] for details). The table shows that the eigenvalue errors converge in an order of 2p when using the (p + 1)-point Gauss rule while the error converges in an order of 2p + 1 and 2p + 2 when using the optimal rule Op for p = 1 and p = 2, respectively. The optimal rules were developed for operators with constant coefficients. It is still an open question to develop optimal rules for the operators with variable coefficients. Herein, for the Schrödinger operator with variable Pöschl-Teller potential, the optimal rules improve the eigenvalue errors significantly but the two-extra orders of convergence are not ensured. 5. Conclusions and future outlook We apply the optimally-blended quadrature rules to approximate the spectrum of a general elliptic differential operator where we account for reaction effects. We show that the optimally blended rules lead to two extra orders of convergence in the eigenvalue errors for both 1D and 3D examples. For a more general problem of (2.1) where the Laplacian − is replaced with an elliptic operator − ∇ · ((x) ∇ u) with diffusion coefficient (x), our numerical tests show that the optimally-blended quadrature rules reduce the eigenvalue errors significantly but do not lead to two extra orders of convergence in multiple dimensions. This is because we lose the Kronecker product structure as in (3.17) when (x) is variable. One future direction is the study on the non-uniform meshes and non-constant coefficient differential eigenvalue problems. Another future direction is to develop the blending techniques for

the isogeometric analysis of multiscale problems (c.f., [36]). The study with variable continuity of the B-spline basis functions is also of interest. We will study the dispersion properties of variable continuity in the basis functions on isogeometric elements and study how the dispersion can be minimized by designing goal-oriented quadrature rules.

Acknowledgement This publication was made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of the Commonwealth Scientific Industrial Research Organisation, CSIRO, of Australia. Additional support was provided by the European Union’s Horizon 2020 Research and Innovation Program of the Marie Skłodowska-Curie grant agreement No. 644202, the Megagrant of the Russian Federation Government (N 14.Y26.31.0013) and the Curtin Institute for Computation. The J. Tinsley Oden Faculty Fellowship Research Program at the Institute for Computational Engineering and Sciences (ICES) of the University of Texas at Austin has partially supported the visits of VMC to ICES.

References [1] M. Ainsworth, H.A. Wajid, Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration, SIAM J. Numer. Anal. 48 (1) (2010) 346–371. [2] P.F. Antonietti, A. Buffa, I. Perugia, Discontinuous Galerkin approximation of the Laplace eigenproblem, Comput. Methods Appl. Mech. Eng. 195 (25) (2006) 3483–3503. [3] I. Babuˇska, J. Osborn, Eigenvalue problems. Handbook of Numerical Analysis, vol. II, 1991, pp. 641–787, North-Holland, Amsterdam. [4] U. Banerjee, A note on the effect of numerical quadrature in finite element eigenvalue approximation, Numer. Math. 61 (1) (1992) 145–152. [5] U. Banerjee, J.E. Osborn, Estimation of the effect of numerical integration in finite element eigenvalue approximation, Numer. Math. 56 (8) (1989) 735–762. [6] U. Banerjee, Analysis of numerical integration in p-version finite element eigenvalue approximation, Numer. Methods Partial Differ. Equ. 8 (4) (1992) 381–394. ˇ V.M. Calo, Gaussian quadrature for splines via homotopy [7] M. Barton, continuation: rules for C2 cubic splines, J. Comput. Appl. Math. 296 (2016) 709–723. ˇ V.M. Calo, Optimal quadrature rules for odd-degree spline spaces [8] M. Barton, and their application to tensor-product-based isogeometric analysis, Comput. Methods Appl. Mech. Eng. 305 (2016) 217–240.

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009

G Model JOCS-879; No. of Pages 7

ARTICLE IN PRESS Q. Deng et al. / Journal of Computational Science xxx (2018) xxx–xxx

ˇ V.M. Calo, Gauss-Galerkin quadrature rules for quadratic and cubic [9] M. Barton, spline spaces and their application to isogeometric analysis, Comput.-Aided Des. 82 (2017) 57–67. ˇ V.M. Calo, Q. Deng, V. Puzyrev, Generalization of the Pythagorean [10] M. Barton, Eigenvalue Error Theorem and its Application to Isogeometric Analysis, 2018 (in press). [11] J.H. Bramble, J.E. Osborn, Rate of convergence estimates for nonselfadjoint eigenvalue approximations, Math. Comput. 27 (123) (1973) 525–549. [12] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Volume 15 of Texts in Applied Mathematics, 3rd ed., Springer, New York, 2008. [13] V. Calo, Q. Deng, V. Puzyrev, Quadrature blending for isogeometric analysis, Proc. Comput. Sci. 108 (2017) 798–807. [14] V.M. Calo, M. Cicuttin, Q. Deng, A. Ern, Spectral Approximation of Elliptic Operators by the Hybrid High-Order Method, 2017 arXiv:1711.01135. [15] V.M. Calo, Q. Deng, V. Puzyrev, Dispersion Optimized Quadratures for Isogeometric Analysis, 2017 arXiv:1702.04540. [16] C. Canuto, Eigenvalue approximations by mixed methods, RAIRO Anal. Numér 12 (1) (1978) 27–50. [17] F. Chatelin, Spectral Approximation of Linear Operators, 1983. [18] P.G. Ciarlet, Finite Element Method for Elliptic Problems, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002. [19] H. Ciftci, R.L. Hall, N. Saad, Construction of exact solutions to eigenvalue problems by the asymptotic iteration method, J. Phys. A: Math. Gen. 38 (5) (2005) 1147. [20] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009. [21] J.A. Cottrell, A. Reali, Y. Bazilevs, T.J.R. Hughes, Isogeometric analysis of structural vibrations, Comput. Methods Appl. Mech. Eng. 195 (41) (2006) 5257–5296. [22] J.D. De Basabe, M.K. Sen, Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations, Geophysics 72 (6) (2007) T81–T95. [23] C. De Boor, A Practical Guide to Splines, vol. 27, Springer-Verlag New York, 1978. ˇ V. Puzyrev, V. Calo, Dispersion-minimizing quadrature [24] Q. Deng, M. Barton, rules for C1 quadratic isogeometric analysis, Comput. Methods Appl. Mech. Eng. 328 (2018) 554–564. [25] Q. Deng, V. Calo, Dispersion-Minimized Mass for Isogeometric Analysis, 2017 arXiv:1711.02979. [26] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, vol. 159, Springer Science & Business Media, 2013. [27] L.C. Evans, Partial Differential Equations, Volume 19 of Graduate Studies in Mathematics, 2nd ed., American Mathematical Society, Providence, RI, 2010. [28] L. Gao, V.M. Calo, Fast isogeometric solvers for explicit dynamics, Comput. Methods Appl. Mech. Eng. 274 (2014) 19–41. [29] S. Giani, hp-Adaptive composite discontinuous Galerkin methods for elliptic eigenvalue problems on complicated domains, Appl. Math. Comput. 267 (2015) 604–617. [30] J. Gopalakrishnan, F. Li, N.-C. Nguyen, J. Peraire, Spectral approximations by the HDG method, Math. Comput. 84 (293) (2015) 1037–1059. [31] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Eng. 194 (39) (2005) 4135–4195. [32] T.J.R. Hughes, J.A. Evans, A. Reali, Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems, Comput. Methods Appl. Mech. Eng. 272 (2014) 290–320. [33] T.J.R. Hughes, A. Reali, G. Sangalli, Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of p-method finite elements with k-method NURBS, Comput. Methods Appl. Mech. Eng. 197 (49) (2008) 4104–4124. [34] B. Mercier, J.E. Osborn, J. Rappaz, P.-A. Raviart, Eigenvalue approximation by mixed and hybrid methods, Math. Comput. 36 (154) (1981) 427–453. [35] B. Mercier, J. Rappaz, Eigenvalue approximation via non-conforming and hybrid finite element methods, Publ. Sém. Math. Inform. Rennes 1978 (S4) (1978) 1–16, Available at http://www.numdam.org/item?id=PSMIR 1978 S4 A10 0.

7

[36] H. Nguyen-Xuan, T. Hoang, V.P. Nguyen, An isogeometric analysis for elliptic homogenization problems, Comput. Math. Appl. 67 (9) (2014) 1722–1741. [37] J.E. Osborn, Spectral approximation for compact operators, Math. Comput. 29 (131) (1975) 712–725. [38] L. Piegl, W. Tiller, The NURBS Book, Springer Science & Business Media, 1997. [39] V. Puzyrev, Q. Deng, V.M. Calo, Dispersion-optimized quadrature rules for isogeometric analysis: modified inner products, their dispersion properties, and optimally blended schemes, Comput. Methods Appl. Mech. Eng. 320 (2017) 421–443. [40] V. Puzyrev, Q. Deng, V.M. Calo, Spectral approximation properties of isogeometric analysis with variable continuity, Comput. Methods Appl. Mech. Eng. 334 (2018) 22–39. [41] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, vol. 12, Springer Science & Business Media, 2013. [42] G. Strang, G.J. Fix, An Analysis of the Finite Element Method, vol. 212, Prentice-Hall Englewood Cliffs, NJ, 1973. Quanling Deng is a researcher at Curtin University, Australia. His research interests include numerical methods for PDEs, scientific computing, and flow in porous media. He has authored more than 10 peer-reviewed journal papers and book chapters, and has given more than 20 talks in international conferences, workshops and seminars.

Vladimir Puzyrev is a Research Fellow at Curtin University, Australia. His research interests include numerical methods for PDEs, inverse problems in geosciences, sparse linear solvers, and high-performance computing. He has published 20 journal papers and book chapters, and has given more than 40 talks in international conferences, workshops and seminars.

Victor M. Calo was appointed to the position of Professor and CSIRO Professorial Chair in Computational Geoscience, a joint appointment with CSIRO and Curtin University in February 2016. Before working in Western Australia, he taught at the King Abdullah University of Science and Technology in Makkah, Saudi Arabia from 2009–2016. His research interests include modelling and simulation of geomechanics, fluid dynamics, flow in porous media, phase separation, fluid-structure interactions, solid mechanics and high-performance computing. The technologies Professor Calo is currently developing in collaboration with other researchers, aim to simplify the use of state-of-the-art techniques for high-performance computing. These methods can model geological processes that are relevant to the resource extraction industry. Professor Calo is a Highly-Cited Researcher in the Academic Ranking of World Universities of the Shanghai Jiao Tong University and Thomson Reuters and has authored more than 200 peer-reviewed publications. In the past two years, he has given more than 20 invited presentations and keynotes at conferences and seminars, and organized 10 mini symposia at international conferences.

Please cite this article in press as: Q. Deng, et al., Isogeometric spectral approximation for elliptic differential operators, J. Comput. Sci. (2018), https://doi.org/10.1016/j.jocs.2018.05.009