Fast multigrid solvers for conforming and nonconforming multipatch Isogeometric Analysis
Abstract
Isogeometric Analysis is a highorder discretization method for boundary value problems that uses a number of degrees of freedom which is as small as for a loworder method. Standard isogeometric discretizations require a global parameterization of the computational domain. In nontrivial cases, the domain is decomposed into patches having separate parameterizations and separate discretization spaces. If the discretization spaces agree on the interfaces between the patches, the coupling can be done in a conforming way. Otherwise, nonconforming discretizations (utilizing discontinuous Galerkin approaches) are required. The author and his coworkers have previously introduced multigrid solvers for Isogeometric Analysis for the conforming case. In the present paper, these results are extended to the nonconforming case. Moreover, it is shown that the multigrid solves get even more powerful if the proposed smoother is combined with a (standard) GaussSeidel smoother.
keywords:
Isogeometric Analysis, Multipatch domains, Symmetric interior penalty discontinuous GalerkinremarkRemark \newproofproofProof
1 Introduction
Isogeometric Analysis (IgA), see Hughes:2005 (), is a spline based approach for approximating the solution of a boundary value problem (BVP). One of the big strengths of IgA is that it has the approximation power of a highorder method while the number of degrees of freedom behaves basically like for a loworder method. To obtain this behavior, we have to be able to increase the spline degree while we simultaneously increase the smoothness. In IgA, this is typically called refinement and leads to spline based discretizations.
In IgA, the computational domain (usually called physical domain) is parameterized using spline or NURBS functions. Since it might be too restrictive to parameterize the whole computational domain using just one global (smooth) geometry function, one typically represents the physical domain as the union of subdomains, in IgA called patches. Then, each of the patches is parameterized by its own geometry function (multipatch IgA).
On each patch, a space of trial and test functions is introduced. The simplest approach to set up such a function space is to use tensorproduct Bsplines on the unit square or unit cube (parameter domain) and to use the geometry function to map them onto the physical domain or, in the multipatch case, onto one of the patches. If we set up the function spaces such that the basis functions on the interfaces between the patches agree, we can use conforming discretizations. Approximation errors (cf. Hughes:2005 (); Bazilevs:BeiraoDaVeiga:Cottrell:Hughes:Sangalli:2006 (); BeiraoDaVeiga:Buffa:Rivas:Sangalli:2011 (); BeiraoDaVeiga:Buffa:Sangalli:Vazquez:2014 (); Takacs:Takacs:2015 (); Floater:Sande:2017 (); Sande:Manni:Speleers:2018 (); deLaRiva:Rodrigo:Gaspar:2018 (); Tielen:Moller:Goddeke:Vuik:2019 () and many others) and multigrid solvers (cf. Gahalaut:Kraus:Tomar:2013 (); Speleers:2015 (); Hofreither:Takacs:Zulehner:2017 (); Hofreither:Takacs:2017 () and others) for such discretizations have been previously discussed. Since we are interested in refinement, we need results that are explicit in the spline degree. For the singlepatch case, such error estimates have originally been given in Takacs:Takacs:2015 () and later improved in Sande:Manni:Speleers:2018 (). In Hofreither:Takacs:2017 (), a robust singlepatch multigrid solver has been proposed and analyzed based on the error estimates from Takacs:Takacs:2015 (). In Takacs:2018 (), both the approximation error estimates and the multigrid solver have been extended to the conforming multipatch case. These results are the foundation of the present paper.
If conforming discretizations are not feasible, discontinuous Galerkin (dG) approaches are possible. One standard dG approach is the Symmetric Interior Penalty discontinuous Galerkin (SIPG) method, see Arnold:1982 (); Arnold:Brezzi:Cockburn:Marini:2002 (). Already in LMMT:2015 (); LT:2015 (), it has been proposed to utilize these approaches to couple patches in IgA. Recently, also the dependence of the approximation error on the spline degree has been analyzed, see Takacs:2019 (). It was not possible to show that the approximation error is robust in the spline degree but it could be proven that it only grows logarithmically.
(Robust) multigrid solvers for such nonconforming discretizations are not known so far. In the present paper, it is shown how the multigrid solver from Takacs:2018 () can be extended to SIPG discretizations; we observe – as in Takacs:2018 () – that the numerical experiments show both robustness in the grid size and the spline degree. For completeness, we also show how to extend the convergence analysis from Takacs:2018 () to SIPG discretizations. It is worth noting that there are alternative solvers for multipatch IgA, like FETItype approaches, cf. KleissEtAl:2012 (); Hofer:2017 () and others, overlapping Schwarz type methods, cf. BeiraoDaVeiga:Cho:Pavarino:Scacchi:2012 (), or BDDC methods, cf. BeiraoDaVeiga:Cho:Pavarino:Scacchi:2013 (); most of them, however, have not been worked out for the nonconforming case.
Note that the idea behind the proposed subspace corrected mass smoother is that the boundary value problem on the physical domain (on one patch) can be well approximated by a boundary value problem on the parameter domain. Thus, the tensorproduct structure on the parameter domain can be used. This is true if the geometry function is not too distorted. Otherwise, the convergence behavior suffers significantly. The same behavior can be observed by other fast solvers that are based on the same idea, cf. the fast diagonalization method Sangalli:2016 (). Here, the authors have improved their method by incorporating the geometry information into the preconditioner, cf. Montardini:Sangalli:Tani:2018 (). For the multigrid setting, it has turned out that one can overcome these problems quite well if the subspace corrected mass smother is combined with a GaussSeidel smoother (hybrid smoother) since both approaches have strengths that seem to be somewhat orthogonal to each other (robustness in spline degree vs. robustness in the geometry), cf. also Sogn:Takacs:2018 ().
In the present paper, we illustrate our findings with numerical experiments. All presented numerical experiments are available in the G+Smo library gismoweb ().
This paper is organized as follows. We give the model problem and a conforming discretization in Section 2. Then, in Section 3, we discuss why a nonconforming discretization might be of interest. Moreover, we propose a discontinuous Galerkin approach that fits our needs. We proceed to multigrid solvers: In Section 4, we discuss GaussSeidel smoothers and their performance. Motivated by that section, we introduce a subspace corrected mass smoother in Section 5 and finally a hybrid smoother in Section 6. In Section 7, we conclude and give some outlook. The Appending finally contains the proofs of the theorems stated in the paper.
2 Model problem and standard Galerkin discretization
Let with be an open and simply connected Lipschitz domain. Most of the numerical experiments are done for the twodimensional domains shown in Figure 1.
The first domain is an Lshaped domain consisting of three quadrilaterals. Here, the geometry function is just the identity or a translation. On the coarsest grid level , each patch consists only of one element, i.e., the local basis functions are polynomials only. The second domain is the Yeti footprint which consists of the 21 patches depicted in Figure 1. Here, the grid on the coarsest grid level is as follows. The five patches at the bottom consist of two elements each, which are constructed by subdividing the patches on their longer sides. The remaining patches consist only of one element each. The grid levels are obtained by uniform refinement.
Consider the following Poisson model problem with Dirichlet boundary conditions. Find such that
(1) 
where and are given functions. Here and in what follows, , and are the standard Lebesgue and Sobolev spaces.
The experiments are performed for the choice and ; note that is the exact solution of the problem.
After homogenization (, ), the problem reads in variational form as follows. Find such that
(2) 
The computational domain is a standard multipatch domain. Thus, we assume that is composed of nonoverlapping patches :
(3) 
where each patch is represented by a sufficiently smooth bijective geometry function
(4) 
which can be continuously extended to , the closure of . Moreover, we assume that the mesh introduced by the patches satisfies the following condition.
Assumption 1
For any , the intersection is either (a) empty, (b) one common vertex, (c) the closure of one common edge, or – for – (d) the closure of one common face.
For each of the patches, we assume to have a hierarchy of grids with levels obtained by uniform refinement, which we denote by
(5) 
where is the space of tensorproduct splines of degree , smoothness (or, equivalently, ) and grid size on the parameter domain . Note that the grid can be nonuniform and both the spline degree and the grid size can depend on the spatial direction and of the patch number; for simplicity, we do not write down this dependence explicitly. Note that the grid needs to be quasiuniform, i.e., there needs to be a constant such that all knot spans on grid level are bounded from below by . The functions are assumed to form a (standard) Bspline or NURBS basis of .
To be able to set up a conforming discretization, we need to assume that the function spaces are fully matching on the interfaces, cf. (Takacs:2018, , Assumption 2.4). For tensorproduct Bspline bases, the following assumption characterizes fully matching discretizations.
Assumption 2
On each interface between two patches, the geometry functions, the knot vector in tangential direction, and the spline degree in tangential direction agree.
Assuming a fully matching discretization, we define the conforming discretization space by
(6) 
A basis for this space is visualized in Figure 2 (left), where all basis functions are represented by their Greville point. The support of the basis functions with Greville point in the interior of a patch is completely contained in that patch. The basis functions with Greville points on the interfaces are combinations of the matching patchlocal basis functions. Their support extends to the vertices if and only if the Greville point is located on the vertex.
The conforming discretization of the model problem is obtained using the standard Galerkin principle: Find such that
(7) 
Using the abovementioned basis for the space , we obtain a standard matrixvector problem: Find such that
(8) 
where is the stiffness matrix, the vector is the representation of with respect to the chosen basis and the load vector is obtained by testing the function with the basis functions.
3 Symmetric interior penalty discontinuous Galerkin (SIPG) discretization
Following LMMT:2015 (); LT:2015 (); Takacs:2019 (), we use a conforming isogeometric discretization for each patch and couple the patches using discontinuous Galerkin. We assume that the domain is again subdivided into patches such that (3), (4) and Assumption 1 are satisfied. We assume again to have patchlocal spaces as in (5), which are combined in a nonconforming (i.e., discontinuous) way, i.e., we just define
(9) 
This allows us to drop Assumption 2. Note that we strongly enforce the Dirichlet boundary conditions in our example; alternatively, one could use the SIPG method also to enforce the Dirichlet boundary conditions.
Since we have a discontinuous function space, we can visualize the degrees of freedom by tearing apart the individual patches, cf. Figure 3. Here, neither the Greville points nor the basis functions need to agree on the interfaces; the support of each basis function is contained in one single patch.
Since , it is not feasible to use the standard Galerkin principle for discretization. Thus, we couple the patches using the Symmetric Interior Penalty discontinuous Galerkin (SIPG) method. First, we define
to be the set of interfaceindices. For each interface with , we define the following symbols.

n is the outer normal vector of . (Thus, is the outer normal vector of .)

is the jump operator: .

is the averaging operator: .
Now, we can formulate the SIPG discretization as follows: Find such that
(10) 
where we define
(11)  
There is some independent of the grid size, the spline degree and the number of patches such that for all , the bilinear form is coercive, cf. (Takacs:2019, , Theorems 8 and 9). Thus, for , the Theorem of Lax Milgram states that the problem (10) has exactly one solution. The combination of Ceá’s Lemma and a naive approximation error estimate yields a discretization error estimate of the form
cf. Takacs:2019 (). By doing a more careful analysis, we obtain estimates of the form
see (Takacs:2019, , eq. (15)). This significantly decreases the influence of the spline degree.
Note that the penalization term has the form
i.e., it depends on the grid size on the finest grid . This follows the ideas from Gopalakrishnan:Kanschat:2003 (). The idea behind that is that
(12) 
hold, i.e., we obtain a multigrid solver with conforming coarsegrid correction. This means that – on the coarse grid levels – the discretization is over penalized by a factor of , i.e.,
where is the canonical parameter and is the chosen one. We will see that this does not cause any problems for the examples we consider; following Gopalakrishnan:Kanschat:2003 (), convergence theory only holds if the number of smoothing steps is sufficiently increased for the coarser grid levels, cf. Remark 5.
Using a basis for the space , we obtain a standard matrixvector problem: Find such that
(13) 
4 Multigrid solvers with GaussSeidel smoothers
In this and the following sections, we discuss several possible choices of multigrid smoothers, illustrate their convergence behavior with numerical experiments, and comment on the convergence theory.
We consider conforming discretizations and nonconforming discretizations which are set up as discussed in the last two sections. As we have nested spaces in all cases, the matrix is always the canonical embedding from into and the restriction matrix is its transpose. The chosen method is presented as pseudocode as Algorithm 1, where we choose for the Vcycle method or for the Wcycle method.
In the finite element world, GaussSeidel smoothers are known to be very efficient smoothers; thus, as a first attempt, we consider such a smoother. One forward GaussSeidel sweep can be represented by
where is a lowertriangular matrix containing the coefficients of the stiffness matrix , i.e., it is given by
To be able to use our multigrid solver as preconditioner for a conjugate gradient solver, the postsmoothing procedure uses the transposed matrix , which represents a backwardGaussSeidel sweep.
All tables show the number of iterations required until the stopping criterion
is satisfied.
As usual, the convergence behavior of the overall solver can be improved if the multigrid method is not just used directly as a solver, but as a preconditioner within a preconditioned conjugate gradient (PCG) solver. Thus, we present results for both possibilities; in the following sections we will restrict ourselves to the more efficient PCG solver. Since the Vcycle and the Wcycle methods yield comparable iteration counts, we present the results for the more efficient Vcycle only. The number of smoothing steps is chosen as in all cases.
The multigrid solver was implemented in C++ based on the G+Smo library gismoweb (). The tables shown in the remainder of this section are obtained with the following command line instructions, where the values and are substituted accordingly.
The results for the PCG experiments, presented in Tables 1 (b) and 2 (b), are obtained by replacing the option i d by the option i cg.
(a) Direct – Conforming  (b) PCG – Conforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
4  9  24  74  –  –  –  –  8  15  28  53  –  –  – 
5  9  24  73  –  –  –  –  8  15  28  52  –  –  – 
6  9  24  72  –  –  –  –  8  15  28  53  –  –  – 
7  10  24  72  –  –  –  –  8  15  28  54  –  –  – 
8  10  24  72  –  –  –  –  8  15  28  54  –  –  – 
(a) Direct – Conforming  (b) PCG – Conforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
3  13  25  75  –  –  –  –  10  15  28  54  –  –  – 
4  14  25  74  –  –  –  –  10  15  28  53  –  –  – 
5  15  25  74  –  –  –  –  11  16  28  54  –  –  – 
6  15  25  72  –  –  –  –  11  16  29  54  –  –  – 
7  17  25  73  –  –  –  –  12  16  29  55  –  –  – 
In Table 1, we observe that the multigrid solver is certainly robust in the grid size. While this approach is very efficient for small numbers of spline degrees, we observe that the convergence rates deteriorate significantly if the spline degree is increased. On the right side of the table, one can see the iteration counts for a preconditioned conjugate gradient method where one Vcycle of the mentioned multigrid method is used as a preconditioner. We observe that the iteration counts are significantly smaller than the iteration counts obtained by directly applying the multigrid solver. However, we simultaneously observe that we do not observe any qualitative improvement. In Table 2, we give the iteration counts for the Yeti footprint. We observe that – despite the fact that the geometry function is now nontrivial – the iteration counts are very similar to those of the Lshaped domain.
When one turns to the nonconforming discretizations, it immediately turns out that the multigrid solver utilizing the GaussSeidel smoother does not converge well at all.
One can show using standard arguments that the multigrid method converges with rates that are independent of the grid size and of the number of patches. The convergence analysis (for the conforming case) employs estimates that increase exponentially in the spline degree, cf. Gahalaut:Kraus:Tomar:2013 (). The numerical experiments show that this is not only a matter of the proof.
Since we did not obtain convincing results, we are interested in more advanced smoothers that work well also for SIPG discretizations and which do not deteriorate if the spline degree is increased.
5 Multigrid with subspace corrected mass smoother
In this section, we employ the subspace corrected mass smoother as introduced in Hofreither:Takacs:2017 (). That smoother requires that the spline space has a tensorproduct structure; in our examples, we have such a structure on each patch but not on the overall domain. The extension of that smoother to conforming discretizations has been discussed in Takacs:2018 () based on a domaindecomposition approach. The key idea was to decompose all degrees of freedom on a perpiece bases. Pieces are the patchinteriors and the interface pieces. In two dimensions, the interface pieces are the edges and the vertices of each of the patches. In three dimensions, the interface pieces are the faces, the edges and the vertices of each of the patches. The decomposition of the degrees of freedom is depicted in Figure 4 (left).
The piecelocal smoothers, which we denote by , are defined as follows. For being a patchinterior, we choose to be the subspace corrected mass smoother as proposed in Hofreither:Takacs:2017 (). We choose the scaling parameter (which was called in Hofreither:Takacs:2017 ()) to be for some suitable chosen parameter . If is an interface piece, we choose
where the matrix represents the embedding of the piece in the whole space. The symbol refers to be the application of a direct solver. Applying a direct solver on the interfaces is feasible since the interfaces have much smaller numbers of degrees of freedom than the interiors of the patches. The overall smoother is just an additive composition of the piecelocal smoothers, i.e., we choose
(14) 
where the sum is taken over all pieces . Here, is some damping parameter to be chosen.
The convergence theory from Takacs:2018 () can be summarized by the following theorem.
Theorem 3
Assume that is such that full elliptic regularity is satisfied (cf. (Takacs:2018, , Assumption 3.1)). Consider the conforming discretization and a multigrid solver with the smoother (14). There are constants , and which are independent of , , and (but may depend particularly on the geometry functions and the maximum number of neighbors of a patch) such that for
(15) 
the Wcycle multigrid method converges with a convergence rate .
Note that the terms and imply that the convergence degrades if too small values of those parameters are chosen. Thus, it is of interest to choose these parameters in a rather optimal way.
Note that the convergence theorem requires full elliptic regularity (cf. (Takacs:2018, , Assumption 3.1)). Thus, it is applicable to the Yeti footprint but it is not directly applicable to the Lshaped domain. Convergence results for the case with full elliptic regularity often carry over in practice to cases where that regularity assumption does not hold. The same behavior can be observed for the numerical experiments we have considered.
Observe that the convergence theorem suggests that the number of smoothing steps should increase with . As already outlined in Takacs:2018 (), this seems to be too pessimistic since the numerical experiments have shown that in all cases yields good convergence rates.
The numerical experiments are again applied within the same setup as in the last section. We set up a Vcycle multigrid method with 1+1 smoothing steps of the proposed smoother (on all grid levels). The damping parameter is chosen as indicated with the option MG.Damping and the scaling parameter is chosen as indicated with the option MG.Scaling below. The tables for the conforming case shown in this section are obtained with the following code, where the values and are substituted accordingly:
⬇ > ./multiGrid_example g domain2d/ldomain.xml r p s scms MG.Damping 1 MG.Scaling .12 i d // Table 3 (a) > ./multiGrid_example g domain2d/yeti_mp2.xml r p s scms MG.Damping .25 MG.Scaling .2 i cg // Table 5 (a)The results for the PCG experiments presented Table 3 (b) are obtained by replacing the option i d by the option i cg.
(a) Direct – Conforming  (b) PCG – Conforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
4  27  23  22  20  17  16  15  16  14  13  12  11  11  10 
5  29  27  27  26  24  24  22  17  16  15  15  14  14  14 
6  30  30  28  27  27  27  26  17  17  16  16  15  16  15 
7  31  30  29  28  28  27  27  17  17  17  16  16  16  16 
8  32  31  30  29  28  28  28  18  17  17  17  16  16  16 
All numerical experiments show that the proposed method is robust both in the grid size and the spline degree. However, when comparing the results for the Yeti footprint from Table 5 (a) with the corresponding results for the Lshaped domain from Table 3 (b), we see that the multigrid solver suffers from distorted geometry functions.
The numbers for the Yeti footprint seem not to be completely robust in the grid size. Since we have given a convergence theory, we know that the convergence numbers are bounded uniformly. Thus, the observed behavior is preasymptotic. The reason for this is that on coarser grid levels, the geometry is not resolved exactly. Let be the original stiffness matrix and be the simplified stiffness matrix obtained by neglecting the geometry function. Then, we have
which yields
(16)  
which can be finally bounded by a constant times some power of the quantity . Of none of these relations is satisfied by equality. In such a case, the iteration counts are likely to increase if the grid gets refined. For more on this topic, see (Sogn:2018, , Section 7.4).
As a next step, we turn towards the nonconforming discretizations. Here, each degree of freedom is assigned to exactly one patch. So, it would be tempting to set up a patchwise splitting of the degrees of freedom. Unfortunately, numerical experiments have shown that this approach does not work well. So, we follow the approach from Takacs:2018 () also in the nonconforming case and split the degrees of freedom again into pieces . This means that we avoid breaking the coupling which was enforced by the penalty term. So, the degrees of freedom belonging to one edge (face, vertex) are considered to be one piece, even if the degrees of freedom belong to different patches, see Figure 4 (right).
For this choice, we can give the following convergence theorem.
Theorem 4
Assume that is such that full elliptic regularity holds (cf. (Takacs:2019, , Assumption 4)) and assume that the geometry functions (but not necessarily the discretizations) agree on the interfaces (cf. (Takacs:2019, , Assumption 2)). Consider the SIPG discretization and a multigrid solver with the smoother (14). There are constants , and which are independent of , , and (but may depend particularly on the geometry functions and the maximum number of neighbors of a patch) such that for
(17) 
the Wcycle multigrid method converges with a convergence rate .
We give the proof of this theorem in the Appendix; the proof is based on the error estimates from Takacs:2019 ().
One might observe that the number of smoothing steps required by this convergence theorem increases like . This follows the approach suggested in Gopalakrishnan:Kanschat:2003 () and is related to the chosen overpenalization discussed in Section 3. {remark} Note that the number of degrees of freedom on the coarser grid levels is smaller by a factor of . So, also using these additional smoothing steps, the overall complexity of the multigrid solver is still linear in the number of unknowns on the finest grid level if (a) or (b) the Vcycle is considered. If we consider and the Wcycle, the choice (17) yields that the computational complexity grows like , where is the number of unknowns on the finest grid level. In the numerical experiments, we did not observe that increasing the number of smoothing steps has been required. Analogously to the conforming case, also the stated dependence on is too pessimistic; thus, we again choose on all grid levels.
Now, we provide numerical experiments for the SIPG discretization. Theoretically, we could just use exactly the discretization that has been chosen for the conforming case. This, however, yields a (particularly uninteresting) special case since Assumption 2 holds. In this special case, we have and the SIPG formulation converges to the conforming discretization for . Instead, we are interested in a discretization such that Assumption 2 does not hold: We modify the setup of the spaces. For one third of the patches, we use the original spline space . For one third of the patches, we use the spline space . For the last third of the patches, we use the spline space . This particular setting is obtained with the command line option NonMatching. In this way, we obtain a setup where a conforming discretization is not possible.
The tables for the nonconforming case shown in this section are obtained with the following code, where the values and are substituted accordingly:
⬇ > ./multiGrid_example g domain2d/ldomain.xml r p DG NonMatching s scms MG.Damping .9 MG.Scaling .12 i d // Table 4 (a) > ./multiGrid_example g domain2d/yeti_mp2.xml r p DG NonMatching s scms MG.Damping .25 MG.Scaling .2 i cg // Table 5 (b)The results for the PCG experiments presented in Table 4 (b) are obtained by replacing the option i d by the option i cg.
(a) Direct – Nonconforming  (b) PCG – Nonconforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
4  22  28  34  33  23  35  23  17  16  15  14  13  12  12 
5  71  48  45  69  35  32  57  19  19  18  17  17  17  16 
6  73  71  70  46  69  57  145  21  20  20  19  19  20  19 
7  100  106  86  71  92  67  61  22  22  21  21  21  22  22 
8  90  94  127  98  291  106  73  23  23  22  22  22  22  23 
(a) PCG – Conforming  (b) PCG – Nonconforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
3  44  42  41  39  37  36  34  40  38  36  35  34  33  31 
4  48  47  45  43  43  40  41  44  44  42  42  40  40  39 
5  51  49  48  47  45  45  44  49  47  47  46  46  45  44 
6  52  51  49  48  47  46  45  58  57  57  56  55  54  53 
7  54  53  51  50  49  48  47  74  73  72  71  71  72  70 
The PCG discretizations are presented in Tables 4 (b) and 5 (b); we again obtain robustness in the grid size and the spline degree. Here, for the Yeti footprint, we have again iteration counts that are slightly increasing with the grid size; again, this observation can be explained by the fact that finer grids allow to resolve the geometry function better, cf. (16).
In principle, the method works also if the multigrid solver is applied directly, cf. Table 4 (a). Here, we suffer from numerical instabilities which are amplified with an increasing number of levels. One can avoid these instabilities, e.g., by increasing the number of pre and post smoothing steps. However, using the multigrid method as a preconditioner within a PCG solver is obviously the more efficient approach.
6 Multigrid with hybrid smoother
We have observed that a multigrid method with the subspace corrected mass smoother is robust in the grid size and the spline degree and works well for both conforming and discontinuous Galerkin discretizations. We have also observed that this approach suffers from nonsimple geometry functions since it is based on the close connection between the stiffness matrix and the simplified stiffness matrix . The results for the GaussSeidel smoother are different: the multigrid solver works badly both for large spline degrees and for discontinuous Galerkin discretizations. However, by comparing Table 1 with Table 2, we observe that the method behaves quite robust in the geometry function.
Since the behavior of the two smothers is somewhat orthogonal, we can hope for a good method if we combine them. Our idea is to use one forward GaussSeidel sweep followed by the subspace corrected mass smoother for presmoothing and the subspace corrected mass smoother followed by one backward GaussSeidel sweep for postsmoothing. The overall method is presented as Algorithm 2.
The convergence analysis from Section 5 can be easily carried over to the hybrid smoother. The iteration matrix for the (V or W cycle) multigrid method with the hybrid smoother is given by
where is the iteration matrix of the (V or W cycle, respectively) multigrid method with the subspace corrected mass smoother. Since the GaussSeidel iteration is stable in the energy norm, we obtain
So, we have using the results from the last section the convergence of the Wcycle multigrid method with hybrid smoother. Thus, we obtain as follows.
Corollary 5
(a) Direct – Conforming  (b) PCG – Conforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
3  12  15  22  29  35  40  47  9  11  14  17  20  21  22 
4  13  16  23  32  38  45  50  10  11  15  19  21  23  25 
5  14  16  24  35  41  47  53  11  11  16  19  22  24  26 
6  15  16  26  37  45  52  55  11  12  16  20  23  26  27 
7  16  17  26  38  47  54  57  12  12  17  21  24  26  28 
(a) Direct – Nonconforming  (b) PCG – Nonconforming  
2  3  4  5  6  7  8  2  3  4  5  6  7  8  
3  23  18  35  34  37  24  28  15  17  18  19  20  20  21 
4  19  23  31  42  32  59  32  16  19  21  22  24  26  26 
5  20  25  34  50  46  67  62  17  20  22  25  27  29  31 
6  22  27  30  47  48  52  69  17  20  23  25  27  30  32 
7  22  29  33  40  58  52  65  17  21  23  26  29  31  33 
The tables for the experiments with the hybrid smoother are obtained with the following code, where the values and are substituted accordingly:
⬇ > ./multiGrid_example g domain2d/yeti_mp2.xml r p s hyb MG.Damping .25 MG.Scaling .1 i d // Table 6 (a) > ./multiGrid_example g domain2d/yeti_mp2.xml r p DG NonMatching s hyb MG.Damping .25 MG.Scaling .1 i d // Table 7 (a)The results for the PCG experiments, presented in Tables 6 (b) and 7 (b), are obtained by replacing the option i d by the option i cg.
For both cases, we obtain that the iteration counts are quite robust in the grid size (even if the maximum number of iterations is not reached for the coarser grid levels). We observe that the number of iterations increases with the spline degree. This is indeed due to the fact that for small values of , the GaussSeidel smoother yields very fast convergence and that convergence behavior is carried over to the hybrid smoother. For larger spline degrees, the hybrid smoother’s convergence behavior degrades mildly; this is due to the fact that the GaussSeidel smoother is not completely capable to capture all effects perfectly. Still, keeping in mind that the condition number of the stiffness matrix grows exponentially with the spline degree, the observed behavior is still very satisfactory.
Compared to applying the subspace corrected mass smoother only, the hybrid smoother pays of in most cases. Certainly, applying the hybrid smoother with means basically that pre and postsmoothing steps are applied. Since the GaussSeidel smoother is slightly cheaper than the subspace corrected mass smoother, the costs for one such cycle are smaller than the costs of two multigrid cycles with the subspace corrected mass smoother only.
Besides the twodimensional examples considered so far, the proposed methods can be directly extended to three dimensional problems (even if the details of the convergence theory have not been worked out for these cases). We consider the two domains depicted in Figure 5: the Fichera corner and a variant of that domain with nontrivial geometry function, which we call twisted Fichera corner.
The tables for the three dimensional domains are obtained with the following code, where the values and are substituted accordingly:
⬇ > ./multiGrid_example g domain2d/fichera.xml r p s hyb MG.Scaling .12 MG.Damping 1 i cg // Table 8 (a) > ./multiGrid_example g domain2d/twisted_fichera.xml r p s hyb MG.Scaling .12 MG.Damping .25 i cg // Table 9 (a)The results for the DG experiments, presented in the Tables 8 (b) and 9 (b), are obtained by adding the command line options DG NonMatching.
(a) PCG – Conforming  (b) PCG – Nonconforming  
2  3  4  5  6  2  3  4  5  6  
2  6  6  7  8  8  12  13  15  16  19 
3  6  7  7  8  8  13  15  16  18  20 
4  6  8  9  9  10  14  16  18  19  20 
5  6  8  9  10  11  14  16  17  19  20 
(a) PCG – Conforming  (b) PCG – Nonconforming  
2  3  4  5  6  2  3  4  5  6  
2  10  13  14  15  17  22  25  28  30  31 
3  13  15  18  20  22  29  31  33  36  38 
4  14  17  19  22  25  31  34  37  40  42 
5  16  17  20  23  26  31  36  40  44  47 
Similar to the results for the Yeti footprint, Tables 8 and 9 again show small iteration counts. For the twisted Fichera corner, we observe that the number of iterations increases mildly when the grid gets refined; this is again to be explained by the better resolution of the geometry function. Moreover, we observe a very mild dependence on the spline degree.
7 Conclusions and outlook
We have presented robust multigrid solvers for multipatch IgA with conforming and nonconforming discretizations. We have given convergence results that exactly state the robustness of the solvers in the grid size. Concerning the dependence on the spline degree, the statements seem to be too pessimistic since the solvers have been completely or (at least) rather robust in practice.
We have addressed another issue, which causes problems for all solvers that use the tensorproduct structure on the parameter domain: the dependence on the geometry function. We have proposed a hybrid smoother between our subspace corrected mass smoother and the GaussSeidel smoother which seems to reduce the effect on the geometry function. Finding approaches to better incorporate the geometry function into the smoother itself seems to be an interesting topic for further research.
Appendix
In the appendix, we give a proof of Theorem 4 and of some auxiliary results.
Every constant used within the appendix is assumed to be independent of the grid size, the grid level, the spline degree and the number of patches, but it may depend on the geometry function (cf. (Takacs:2019, , Assumption 3)), the number of neighbors of a patch (cf. (Takacs:2018, , Assumption 2.3)), the constant in the elliptic regularity assumption (cf. (Takacs:2018, , Assumption 3.1)) and the quasiuniformity of the grid, i.e., the ratio between the largest and the smallest knot span of the knot vectors on one level. We write if and only if there is a constant such that and we write if and only if and .
First, we show the following lemma, which is basically a trace inequality.
Lemma 6
Let be the space of tensorproduct splines of degree on a quasiuniform grid with size on the parameter domain . Then, the estimate
holds for all and all .
For , let be the eigenfunctions of , i.e., such that
where is the Kronecker delta and are the corresponding eigenvalues. Using coercivity of and a standard inverse estimate, cf. (Schwab:1998, , Corollary 3.94), we obtain
We define level sets
for , where
is the number of level sets. Note that by construction every eigenvalue belongs to exactly one level set. Every function can be represented as
A standard trace estimate, cf. (Takacs:2018, , Lemma 4.4), yields
(18)  
where
for . Since and since all eigenvalues are in or , respectively, we obtain
Using (18), we obtain further
Finally, the CauchySchwarz inequality and orthogonality of the basis functions (both in and ) yield
which finishes the proof. ∎
Now, we give bounds on the smoother which allow to show the smoothing property.
Lemma 7
The proof of this Lemma requires the notation from Takacs:2018 (), i.e., we denote the set of all patchinteriors by , the set of all edges by and the set of all vertices by .
Observe that we have
(19) 
for all , where is defined by (11) and (12). For , this statement directly follows from (Takacs:2019, , Theorem 8). Since (Takacs:2019, , Theorem 8) also holds in cases of overpenalization, we can apply that theorem also to the case and obtain (19) also in that case.
As first step, we bound from below. (The following arguments are analogous to (Takacs:2018, , Lemma 4.3).) The triangle inequality yields
(20) 
For , (Hofreither:Takacs:2017, , Lemma 8) and (19) yield . For , we have by definition . Thus, we obtain from (20)
and for all with small enough further
which shows the first part of the desired inequality.
Now, we bound from above. We use the decomposition
cf. (11). Using (19), we obtain
for all and, therefore,
(21) 
where
and
Completely analogous to (Takacs:2018, , Lemma 4.7), we obtain
(22) 
So, it remains to bound from above. Since the restriction of to any patchinterior vanishes, the same arguments as in the proof of (Takacs:2018, , Lemma 4.7) and the triangle inequality yield