IMPROVED APPROXIMATION OF EIGENVALUES IN ISOGEOMETRIC METHODS FOR MULTI-PATCH GEOMETRIES AND NEUMANN BOUNDARIES THOMASHORGER,ALESSANDROREALI,BARBARAWOHLMUTH, ANDLINUSWUNDERLICH 7 1 0 Abstract. We present a systematic study on isogeometric approximations 2 of second order elliptic eigenvalue problems. In the case of either Neumann n boundaries or multi-patch geometries and a mortar approach spurious eigen- a values-typicallyreferredtoas“outliers”-occur. Thispollutionofthediscrete J spectrum can be avoided by a hybrid approach. In addition to the imposed weakconsistencywithrespecttotheNeumannboundaryandtotheinterface 3 mortarcouplingcondition,weapplyhigher-orderpenaltytermsforthestrong 2 coupling. The mesh dependent weights are selected such that uniform conti- nuity is guaranteed and optimal a priori estimates are preserved. Numerical ] A resultsshowthegoodbehavioroftheproposedmethodinthecontextofsimple 1D and 2D Laplace problems, but also of more complex multi-patch mapped N geometriesandlinearelasticity. . h Isogeometric analysis, Eigenvalues, Mortar methods, Penalty methods, Multi- t patch geometries, Neumann boundaries a m [ 1. Introduction 1 Isogeometric analysis (IGA), introduced in 2005 by Hughes et al. in [1], is a v family of methods using highly regular basis functions typical of CAD systems, 3 like B-splines or non-uniform rational B-splines (NURBS), to construct numerical 5 approximations of partial differential equations (PDEs). The idea of using spline 3 6 functions for the approximation of PDEs can even be found in earlier works (see, 0 e.g.,[2]andreferencestherein). However,inIGAthecomputationaldomainisalso 1. represented in terms of the same functions used for approximation, in an isopara- 0 metric fashion, with the goal of vastly simplifying the mesh generation and refine- 7 ment processes with respect to standard finite element analysis (FEA), possibly 1 bridging the gap between CAD and analysis. : v i X (T. Horger, B. Wohlmuth, L. Wunderlich) Institute for Numerical Mathematics, Tech- nischeUniversita¨tMu¨nchen,Boltzmannstraße3,85748Garchingb. Mu¨nchen,Germany r a (A. Reali) Department of Civil Engineering and Architecture, University of Pavia, via Ferrata 3, 27100 Pavia, Italy (A.Reali)Institute for Advanced Study, Technische Universita¨t Mu¨nchen, Lichten- bergstraße 2a, 85748 Garching b. Mu¨nchen, Germany E-mail addresses: (T. Horger) [email protected], (A. Reali) [email protected], (B.Wohlmuth) [email protected], (L.Wunderlich, Correspondingauthor) [email protected]. TH, BW, and LW would like to gratefully acknowledge the funds provided by the Deutsche Forschungsgemeinschaftunderthecontract/grantnumbers: WO671/11-1,WO671/13-2andWO 671/15-1 (within the Priority Programme SPP 1748, ”Reliable Simulation Techniques in Solid Mechanics. DevelopmentofNon-standardDiscretisationMethods,MechanicalandMathematical Analysis”. AR would like to gratefully acknowledge the funds provided by Fondazione Cariplo - Regione Lombardia through the project “Verso nuovi strumenti di simulazione super veloci ed accuratibasatisull’analisiisogeometrica”,withintheprogramRST-rafforzamento. 1 2 IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION Thankstoitshigh-regularity,IGAhasalsoshowntopossesssuperiorapproxima- tionproperties,onaper-degree-of-freedombasis,withrespecttoclassicalFEA,and has been successfully applied in many fields of computational mechanics, including the approximation of eigenvalue problems. Eigenvalue analysis arises in many im- portant applications in science and engineering, where the amount of eigenvalues of interest can be quite different from application to application. For example in vibro-acoustics one is typically only interested in the first part of the spectrum, while for explicit dynamics the approximation quality of the entire spectrum is rel- evant. ComparedtoFEA,itwasobservedthatIGApossessessuperiorapproxima- tionofeigenvalues, see[3,4,5,6]. Suchstudiesarehoweverlimitedtosingle-patch geometries and Dirichlet boundary conditions. In general, when dealing with non-trivial engineering applications, the isogeo- metric computational domain is constituted by the assembly of multiple patches and thus efficient algorithmic techniques to couple different patches are required. To retain the flexibility of the meshes at the interfaces, mortar methods are a very attractive option, originally introduced for the coupling of non-matching meshes in spectral and finite element methods [7, 8]. Several alternative but equivalent discrete mortar formulations exist [9]. Most implementations are based on the unconstrained saddle point approach since it allows for a local assembling. While mortarfiniteelementformulationsarequiteoftenmotivatedbytheflexibilityofdo- main decomposition techniques or by the robustness with respect to non-matching meshesindynamicapplications,IGAleadsinanaturalwaytoamulti-patchsitua- tion in case of complex geometries, and mortar-type techniques are mandatory for theweakcouplingsee,e.g.,[10,11,12]. Amathematicalstabilityandapriorianal- ysis enlightening the use of different dual spaces can be found in [13]. Higher-order couplingsrecentlygainedattention. Wereferto[14,15]forKirchoff-Loveshells. A discussion of strong C1 couplings in multi-patch settings is given in [16], whereas weak continuity of the normal stress is realized in [17]. Alternative higher-order coupling methods based on least-squares techniques were proposed in [18]. Here, we investigate the influence of multi-patch coupling on the isogeometric eigenvalueapproximation. AscomparedtoIGAwithmaximalregularity,bothweak and strong continuity couplings introduce outlier eigenvalues. This observation motivates us to propose a new hybrid mortar variant. In addition to the weak continuity satisfied in terms of a Lagrange multiplier, we apply a penalty approach forthejumptermsinthehigherderivatives. Theweightsinthepenalizationterms areselectedsuchthatboththeconditionnumbergrowthrateofthealgebraicsystem and the optimal a priori convergence rate are preserved. The same approach is also successfully employed to cure the outliers produced when Neumann boundary conditions are imposed. Thepaperisstructuredasfollows: InSection2,webrieflyreviewtheisogeomet- ric mortar discretization, formulate the eigenvalue problem within the saddle point framework and introduce our hybrid mortar variant. Numerical results illustrate in Section 3 the influence of the penalization in the jump terms. In particular, the higher part of the spectrum can be significantly improved. Finally, in Section 4, conclusions are given. 2. Problem setting and hybrid mortar formulation Let Ω ⊂ Rd, d = 1,2,3, be a bounded domain with Γ¯ ∪Γ¯ = ∂Ω and Γ ∩ D N D Γ =∅. We consider the following Laplace eigenvalue problem with Dirichlet and N IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION 3 Neumann boundary conditions on Γ and Γ , respectively: D N −∆u=λu, in Ω, u=0, on Γ , D ∂ u=0, on Γ . n N In the following, we briefly recapture isogeometric mortar methods, and for more details we refer to [13]. For the ease of presentation, we restrict ourselves to the two dimensional case. The generalization to one or three dimensions follows the same lines. 2.1. Mortar saddle point formulation. Let the domain Ω be decomposed into K non-overlapping subdomains Ω , i.e., k K (cid:91) Ω= Ω , and Ω ∩Ω =∅,i(cid:54)=j. k i j k=1 Here, we limit our presentation to the basic isogeometric concepts and notations used throughout the paper and refer to the literature [5, 19, 20, 21] for more de- tails. Each of the subdomains is a NURBS geometry, i.e., there exists a NURBS parametrization F based on an open knot vector Ξ and a degree p such that Ω k k k is the image of the parametric space Ω(cid:98) =(0,1)d by Fk. We set Np(Ξ ) as the multivariate NURBS space (corresponding to Ω ) in the k k parametric domain, with the standard NURBS basis functions N(cid:98)p . For a set of k,i control points C ∈Rd, i∈I, we define a parametrization of a NURBS surface as k,i a linear combination of the basis functions and control points (cid:88) Fk(ζ)= Ck,iN(cid:98)kp,i(ζ), i∈I and assume the regularity stated in [22, Assumption 3.1]. For 1≤k <k ≤K, we define the interface as the interior of the intersection of 1 2 the boundaries, i.e., γ =∂Ω ∩∂Ω , where γ is open. Let the non-empty k1k2 k1 k2 k1k2 interfaces be enumerated by γ , l=1, ..., L. l For each Ω , we introduce H1(Ω ) = {v ∈ H1(Ω ),v = 0}, where we k ∗ k k k k|∂Ω∩∂Ωk use standard Sobolev spaces, as defined in [23], endowed with their usual norms. InordertosetaglobalfunctionalframeworkonΩ, weconsiderthebrokenSobolev space V =(cid:81)K H1(Ω ), endowed with the broken norm (cid:107)v(cid:107)2 =(cid:80)K (cid:107)v(cid:107)2 . k=1 ∗ k V k=1 H1(Ωk) For any interface γ ⊂ ∂Ω , we define by H−1/2(γ ) the dual space of H1/2(γ ), l k l 00 l whichisthespaceofallfunctionsthatcanbetriviallyextendedon∂Ω \γ byzero k l to an element of H1/2(∂Ω ). k In the following, we set our non-conforming approximation framework. On each subdomain Ω , based on the NURBS parametrization, we introduce the approxi- k mation space V = {v = v ◦F−1,v ∈ Np(Ξ )}. On Ω, we define the discrete k k (cid:98)k k (cid:98)k k product space V = (cid:81)K V ⊂ V, which forms a H1(Ω) non-conforming space h k=1 k discontinuous over the interfaces. On the skeleton Γ=(cid:83)L γ , we define the discrete Lagrange multiplier product l=1 l space M as M = (cid:81)L M . Based on the interface knot vector of one of the h h l=1 l adjacent subdomains, M is a spline space of degree p defined on the interface γ . l l An appropriate local degree reduction performed at the crosspoints guarantees the inf-sup stability of the mortar coupling, see [13] for more details, while preserving optimal order error decay rates. 4 IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION We define the bilinear forms a: V ×V →R and m: V ×V →R, such that K (cid:90) K (cid:90) (cid:88) (cid:88) a(u,v)= ∇u·∇v dx, m(u,v)= uv dx k=1 Ωk k=1 Ωk as well as L (cid:90) (cid:88) b(τ,v)= τ[v] dσ, l l=1 γl where[·] denotesthejumpoverγ . Thesaddlepointformulationoftheisogeometric l l mortareigenvalueproblemreadsasfollows: Find(u ,τ )∈V ×M , λh ∈R, such h (cid:98)h h h that a(u ,v )+b(τ ,v )=λhm(u ,v ), v ∈V , h h (cid:98)h h h h h h b(τ ,u )=0, τ ∈M , h h h h with the normalization m(u ,u ) = 1. The saddle point formulation introduces h h 2dimM spuriouseigenvaluestothespectrum. Werestrictourselvestothephysical h relevant eigenpairs (λh,u ). These are characterized by the fact that they are also h eigenpairs of the constrained mortar formulation, i.e., they satisfy a(u ,v )=λhm(u ,v ), v ∈{V :b(τ ,v )=0,τ ∈M }. h h h h h h h h h h 2.5 3.5 IGA, smooth IGA, smooth IGA, C0-line IGA, C0-line IGA Mortar 3 IGA Mortar 2 2.5 n n h / n h / n 2 1.5 1.5 1 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n/N n/N Figure 1. Normalized spectrum for p = 2. Left: Dirichlet BC. Right: Neumann BC. To get a better feeling for the influence of the mortar coupling, we compare in Figure 1 three different IGA approaches for two different boundary conditions. As reference, we use the classical IGA approach with maximal regularity. Since we set p=2, we have C1-regularity. Additionally, we take into account the classical IGA approachwithreducedregularityatthelinex=1,i.e.,C0-regularity(see[5,Chap- ter3.5]),andthestandardmortarIGAapproach,i.e.,continuityisonlyguaranteed in a weak sense, see Figure 5 for the used initial meshes. On the left, we show the normalized spectrum for Dirichlet boundary conditions on ∂Ω. It can be seen that only the IGA approach with maximal regularity exhibits moderate outliers while bothalternativeshavemoredominantoutliers. Thequalityoftheapproximationis notsignificantlyinfluencedbythechoiceofthecontinuitycondition. Bothoptions, strong or weak continuity, yield quantitatively the same results. However, the ac- curacy of the higher part in the spectrum is considerably improved if additionally strong continuity in the normal derivatives across the line x = 1 is imposed. On the right of Figure 1, we consider Neumann boundary conditions. Here all three approaches show severe outliers. We recall that in contrast to Dirichlet boundary IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION 5 conditions, Neumann boundary conditions are typically not satisfied strongly but only weakly in a variationally consistent way. These preliminary observations motivate us to impose higher regularity terms in the formulation. While on matching meshes strong C1-continuity can be easily guaranteed, the situation is different for non-matching multi-patch meshes. Thus we impose higher regularity by additional penalty terms. As already mentioned weak C1-continuity can be imposed in terms of an additional Lagrange multiplier approach. However this requires a careful selection of the space to guarantee uni- form stability and involves a careful handling of the resulting algebraic system. Here we propose an alternative approach. 2.2. Hybrid mortar formulation. Toimprovetheglobalcontinuity,wepenalize the jump in the normal derivatives across the interfaces of a multi-patch geometry andthefirstnormalderivativeonaNeumannboundarypart. Toavoidlocking, we takethelocalL2-projectionπ0 ontothepiecewiseconstantfunctionsonthe(slave) h boundary mesh. The extra penalty term is given as L p−1 (cid:88) (cid:88) c (u ,v )=C c (u ,v )+ Cmcm(u ,v )+C c (u ,v ) h h h N N h h l l h h CP CP h h l=1m=1 with appropriate penalty constants Cm,C ,C ≥ 0. The smooth interface cou- l N CP pling (cid:90) cm(u ,v )= h2m−1π0([∂mu ] ) π0([∂mv ] ) dσ, l h h s h n h l h n h l γl as well as the Neumann penalty term (cid:90) c (u ,v )= hπ0(∂ u) π0(∂ v) dσ, N h h h n h n ΓN are properly weighted with the local mesh-size h on the slave side, respectively h s on the boundary. The index CP in the bilinear form c refers to contributions CP from the crosspoints. At the end points of each interface and each corner of the Neumann boundary, we introduce additional point evaluations. More precisely, for each interface γ , we add the term l p−1 (cid:88) (cid:88) h2m[∂mu ] (x¯)[∂mv ] (x¯). s n h l n h l x¯∈∂γlm=1 On each corner x¯ of the Neumann boundary, the gradient is zero, so we add the term h2∇u (x¯)∇v (x¯), h h and for each point x¯, where an interface and the Neumann boundary intersect, we add the term h2∂ u (x¯)∂ v (x¯), n h n h where n indicates the outward normal. We then consider the penalized eigenvalue problem a (u ,v )+b(τ ,v )=λhm(u ,v ), v ∈V , h h h (cid:98)h h h h h h b(τ ,u )=0, τ ∈M , h h h h with a (u ,v )=a(u ,v )+c (u ,v ). h h h h h h h h We now show broken H1 continuity of the penalized bilinear form a and for h |Γ | > 0 ellipticity on the kernel of the mortar coupling. The ellipticity trivially D follows from the ellipticity of a, since a (u ,u )=a(u ,u )+c (u ,u )≥a(u ,u )≥α(cid:107)u (cid:107)2 . h h h h h h h h h h h Vh 6 IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION To show continuity, it remains to prove |c (u ,v )| ≤ C(cid:107)u (cid:107) (cid:107)v (cid:107) . For the h h h h Vh h Vh integral terms, using standard estimates and stability of the L2-projection, this reduces to an inverse inequality (cid:107)∂mv (cid:107) ≤Ch1/2−m(cid:107)v (cid:107) for m<p. n h L2(γl) h H1(Ωk) For simplicity, we restrict ourselves to a single domain and one boundary part. Standard trace and inverse inequalities, see [19, Theorem 4.2], yield (cid:107)∂mv (cid:107)2 ≤C(cid:107)v (cid:107) (cid:107)v (cid:107) n h L2(γl) h Hm(Ωk) h Hm+1(Ωk) ≤Ch1−m(cid:107)v (cid:107) h1−(m+1)(cid:107)v (cid:107) . h H1(Ωk) h H1(Ωk) Noting that we have for polynomials an inverse inequality between L∞- and L2- norms, the point evaluations can be handled analogously. Remark 1. We note that the a priori analysis of the new hybrid mortar approach can be easily worked out within the abstract framework of non-conforming finite el- ement techniques. The hybrid form allows us to use the Lemmata of Strang to show optimal order a priori bounds for right hand side problems and then further extend this to the eigenvalue theory analogously to [24] . We point out that the weights in the penalty term are selected such that the condition number of the algebraic system is still O(h−2). 3. Numerical results In this section, we present the eigenvalue approximation in different situations. In particular, after a systematical investigation in simple 1D and 2D settings, we study more complex 2D geometries composed of several patches, and we finally consider also a non-trivial multi-patch example in the framework of linear elas- ticity. The numerical simulations are based on the isogeometric Matlab toolbox GeoPDEs [25, 26]. We compare globally smooth spaces, C0-couplings and the previously introduced higher-order penalty couplings. All results reported in the following are in terms of the normalized discrete eigenvalue λh/λ, which directly relates to the relative error in the eigenvalue since (λh −λ)/λ = λh/λ−1. Note that in the case of pure Neumann boundary conditions, the first eigenvalue is zero, so we exclude it from the spectrum. The numerically obtained eigenvalues can be grouped into physical relevant eigenvaluesandspuriouseigenvaluesduetothecoupling. Weneglecttheseunphysi- calmodesandonlyshowthephysicalpartoftheresultingspectrum. Thenumerical criterion to distinguish between these parts is a relative criterion λh /λh > 100. n+1 n We note that the spurious eigenvalues of the saddle point formulation are infinite. For simple domains, analytical solutions are known for the eigenpairs, otherwise reference solutions need to be computed. On the unit line, the set of eigenvectors √ with pure Dirichlet conditions is given by u (x) = 2sin(nπx), with the corre- n sponding eigenvalue λ = n2π2, n = 1,2,..., see [6]. In the Neumann setting n we also have zero as eigenvalue. The two-dimensional eigenpairs on rectangles are obtained by a tensor product approach. We consider the two-dimensional domain (0,2)×(0,1), with eigenvalues (n2+ n(cid:48)2/4)π2, for n,n(cid:48) = 1,2,... in the pure Dirichlet case and n,n(cid:48) = 0,1,... in the pure Neu- mann case. 3.1. Single-patch 2D example with Neumann boundaries. Aswehaveseen, Neumann boundary conditions induce more significant outliers compared to the Dirichlet case. Thus, we first consider the discrete spectrum of a pure Neumann problem on a single-patch unit square with penalty. In Figure 2, the standard results for p = 2 are compared to those including a penalty for the Neumann IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION 7 3.5 IGA smooth 3 IGA smooth, Neumann penalty n2.5 / hn 2 1.5 1 0 0.2 0.4 0.6 0.8 1 n/N Figure 2. Single-patch 2D example with Neumann boundaries. Normalized discrete spectra for p = 2 with and without the Neu- mann boundary penalty. conditions. We observe a considerable improvement, revealing the potential of the new hybrid variant. By adding a penalty term measuring the violation of the Neumann boundary conditioninastrongform, thecompletespectrumcanbeapproximatedasgoodas in the Dirichlet case. 3.2. One-dimensional results. In this subsection, we report on one-dimensional results obtained with several higher-order penalty couplings. C1 C2 C3 C4 κpenalty/κkink,p=2 6.7e+02 - - - κpenalty/κkink,p=3 1.1e+03 1.1e+03 - - κpenalty/κkink,p=4 1.6e+03 1.6e+03 2.6e+03 - κpenalty/κkink,p=5 2.1e+03 2.8e+03 7.7e+03 7.7e+03 Table1. Relativeconditionnumbersfordifferentansatzandcou- pling degrees. Ratio of the condition number κpenalty for the penalty coupling and the condition number κkink for the C0- interface. Penalty constant chosen as Cm =103−m. l In contrast to a pure penalty approach which typically leads to extremely poor conditionnumbersofthealgebraicsystem,ourhybridvariantisratherrobustwith respect to p. As already noted the condition number has the same growth rate as in the standard mortar setting. This results form the proper mesh dependent weighting. In Table 1, we consider the influence of p on the condition number. As referenceweuseanIGAapproachwithoneC0 point. Theratiodependssensitively on the choice of our penalty parameters but does not deteriorate as the mesh-size tends to zero. Furthermore it should be noted that by decreasing the penalty constant also the condition number decreases. However, also the magnitude of the spurious modes decreases and they get more close to the physical ones. Figure3depictsthedifferentdiscretespectraofthepureDirichletproblemusing splines of degrees 2 and 3. A zoom into the higher frequency part of the spectra is also shown. As it is well-known (see, e.g., [4]), finite elements fail to approximate the higher part of the spectrum while IGA with maximal regularity allows to obtain good 8 IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION 2.2 1.2 IGA, smooth IGA, smooth 2 FEM FEM IGA, one C0 point 1.15 IGA, one C0 point 1.8 IGA, with penalty IGA, with penalty λn λn hλ / n1.6 hλ / n 1.1 1.4 1.05 1.2 1 1 0 0.2 0.4 0.6 0.8 1 0.8 0.85 0.9 0.95 1 n/N n/N 3 1.5 IGA, smooth FEM 2.5 IGA, one C0 point 1.48 IGA, one C1 point hλλ / nn 2 IIGGAA,, wwiitthh CC 12 ppeennaallttyy hλλ / nn11..4446 IFGEAM, smooth IGA, one C0 point 1.5 IGA, one C1 point 1.42 IGA, with C 1 penalty IGA, with C 2 penalty 1 1.4 0 0.2 0.4 0.6 0.8 1 0.99 0.992 0.994 0.996 0.998 1 n/N n/N Figure 3. 1D Dirichlet problem: Top: p = 2, Bottom: p = 3. Left: Entire normalized discrete spectrum. Right: Zoom of the last part of the normalized discrete spectrum. results. However IGA with reduced regularity, e.g., introduced by a C0-line or by a weak mortar coupling across multiple-patches, also introduces outliers at the high frequency end of the spectrum. By our new hybrid variant which additionally penalizes jumps in the normal derivatives, the number of these outliers can be significantly reduced. We obtain quantitatively approximation quality as in the single-patch smooth case. The results for the pure Neumann problem are shown in Figure 4, where we also see outliers for the smooth space. However, the penalty can reduce even these outliers and we end up with better results than with the original smooth spline space. The results for degrees 4 and 5 are similar to those shown for degrees 2 and 3 and, therefore, are not reported here. 3.3. Rectangular domain. In this subsection, we show two-dimensional studies on two different mesh cases depicted in Figure 5, case (a), and Figure 6, case (b). Figure 7 depicts the resulting normalized discrete spectra for the Dirichlet and theNeumannproblemsonthemeshesofcase(a). TheIGAwithmaximalregularity and the IGA with a C0-line are computed on the conforming mesh, whereas the weaklycoupledmortarIGAapproachesarecomputedonthenonconformingmesh. We test the cases p=2,3,4,5 and show only p=2,5. In all cases, the lower frequency part of the spectrum is well approximated. In the higher frequency part of the spectrum (above 75%) we observe, as in the one dimensionalcase,apollutionoftheapproximationduetolow-continuitycouplings. In contrast to the one-dimensional case, these couplings do not just pollute a fixed number of modes, but a higher range of the spectrum due to the multivariate structure of the problem. Nevertheless it is observed that the outliers introduced by the weak or strong C0-coupling can be avoided by including a penalization for IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION 9 2.2 2.5 IGA, smooth IGA, smooth 2 FEM FEM IGA, one C0 point IGA, one C0 point 1.8 IGA, with penalty 2 IGA, with penalty λn λn hλ / n1.6 hλ / n 1.4 1.5 1.2 1 1 0 0.2 0.4 0.6 0.8 1 0.99 0.995 1 n/N n/N 3 5.5 IGA, smooth IGA, smooth 5 FEM FEM 2.5 IGA, one C0 point 4.5 IGA, one C0 point IGA, one C1 point 4 IGA, one C1 point hλλ / nn 2 IIGGAA,, wwiitthh CC 12 ppeennaallttyy hλλ / nn3.53 IIGGAA,, wwiitthh CC 12 ppeennaallttyy 2.5 1.5 2 1.5 1 1 0 0.2 0.4 0.6 0.8 1 0.99 0.995 1 n/N n/N Figure 4. 1D Neumann problem: Top: p = 2, Bottom: p = 3. Left: Entire normalized discrete spectrum. Right: Zoom of the last part of the normalized discrete spectrum. 1 1 0.51 0.5 0.49 0 0 0 0.5 1 1.5 2 0 0.5 1 1.5 2 Figure 5. Mesh case (a). Conforming and nonconforming meshes. 1 1 1 0.7 0.4 0.5 0.5 0.3 0 0 0 0 0.5 1 1.5 2 0 0.5 1 1.3 1.7 2 0 0.5 1 1.5 2 Figure 6. Mesh case (b). Left: First conforming mesh, Middle: Second conforming mesh, Right: Mortar mesh. the jump terms in the normal derivatives. Only very few outliers which remain in some cases. For the Neumann problem, we observe the expected improvement of the outliers even in comparison to the IGA approach with maximal regularity, similar to the single-patch setting considered in Section 3.1. For mesh case (b), the results depicted in Figure 8 are similar to mesh case (a). Furthermore, we observe an additional negative effect of mesh perturbation on the spectrum. Remark 2. The accurate evaluation of the mesh dependent bilinear form c re- h quires, as the accurate evaluation of the mortar coupling bilinear form b, a suitable quadrature formula on a merged mesh. Here we use the same quadrature formula 10 IMPROVED ISOGEOMETRIC EIGENVALUE APPROXIMATION 2.5 10 IGA, smooth IGA, smooth IGA, C0-line 9 IGA, C0-line IGA Mortar IGA Mortar IGA Mortar with C 1 penalty 8 IGA Mortar with C 1 penalty 2 7 IIGGAA MMoorrttaarr wwiitthh CC 23 ppeennaallttyy n n 6 IGA Mortar with C 4 penalty h / n h / n 5 1.5 4 3 2 1 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n/N n/N 3.5 20 IGA, smooth IGA, smooth IGA, C0-line 18 IGA, C0-line 3 IIGGAA MMoorrttaarr with C 1 penalty 16 IIGGAA MMoorrttaarr with C 1 penalty 14 IGA Mortar with C 2 penalty IGA Mortar with C 3 penalty n2.5 n12 IGA Mortar with C 4 penalty h / n h / n10 2 8 6 1.5 4 2 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n/N n/N Figure 7. Normalized discrete spectrum for mesh case (a). Left: p = 2, Right: p = 5. Top: Dirichlet problem, Bottom: Neumann problem. 4 IGA 1, smooth 16 IGA 1, smooth IGA 1 C0-line IGA 1 C0-line 3.5 IGA 2, smooth 14 IGA 2, smooth IGA 2 C0-line IGA 2 C0-line IGA Mortar 12 IGA Mortar 3 IGA Mortar with C 1 penalty IGA Mortar with C 1 penalty h / nn2.5 h / nn108 IIIGGGAAA MMMooorrrtttaaarrr wwwiiittthhh CCC 234 pppeeennnaaallltttyyy 2 6 1.5 4 2 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n/N n/N 5 IGA 1, smooth 30 IGA 1, smooth IGA 1 C0-line IGA 1 C0-line 4.5 IGA 2, smooth 25 IGA 2, smooth IGA 2 C0-line IGA 2 C0-line 4 IGA Mortar IGA Mortar IGA Mortar with C 1 penalty 20 IGA Mortar with C 1 penalty h / nn3.53 h / nn15 IIIGGGAAA MMMooorrrtttaaarrr wwwiiittthhh CCC 234 pppeeennnaaallltttyyy 2.5 10 2 1.5 5 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n/N n/N Figure 8. Normalized discrete spectrum for mesh case (b). Left: p = 2, Right: p = 5. Top: Dirichlet problem, Bottom: Neumann problem. for both bilinear forms. More precisely a Gauss-quadrature with p+1 nodes per ele- ment is used on the merged mesh combining the boundary knots from the slave side