Fast multigrid solvers for conforming and non-conforming multi-patch Isogeometric Analysis

TL;DR

This paper extends multigrid solvers for Isogeometric Analysis to handle both conforming and non-conforming multi-patch discretizations, demonstrating improved efficiency when combined with Gauss-Seidel smoothing.

Contribution

It introduces multigrid solvers for non-conforming multi-patch Isogeometric Analysis and shows enhanced performance with combined smoothing techniques.

Findings

Multigrid solvers effectively handle non-conforming discretizations.

Combining the proposed smoother with Gauss-Seidel improves convergence.

Extended multigrid methods outperform previous conforming-only approaches.

Abstract

Isogeometric Analysis is a high-order discretization method for boundary value problems that uses a number of degrees of freedom which is as small as for a low-order method. Standard isogeometric discretizations require a global parameterization of the computational domain. In non-trivial 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, non-conforming 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 non-conforming case. Moreover, it is shown that the multigrid solves get even more powerful if the proposed smoother is combined with a (standard) Gauss-Seidel smoother.

Figures

Tables

Table 1. Table 1: V-cycle with Gauss-Seidel smoother for the L-shaped domain

	(a) Direct – Conforming							(b) PCG – Conforming
$L\backslash p$	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	–	–	–

Table 2. Table 2: V-cycle with Gauss-Seidel smoother for the Yeti footprint

	(a) Direct – Conforming							(b) PCG – Conforming
$L\backslash p$	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	–	–	–

Table 3. Table 3: V-cycle with subspace corrected mass smoother for the L-shaped domain

	(a) Direct – Conforming							(b) PCG – Conforming
$L\backslash p$	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

Table 4. Table 4: V-cycle with subspace corrected mass smoother for the L-shaped domain

	(a) Direct – Non-conforming							(b) PCG – Non-conforming
$L\backslash p$	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

Table 5. Table 5: V-cycle with subspace corrected mass smoother for the Yeti footprint

	(a) PCG – Conforming							(b) PCG – Non-conforming
$L\backslash p$	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

Table 6. Table 6: V-cycle with hybrid smoothing strategy for the Yeti footprint

	(a) Direct – Conforming							(b) PCG – Conforming
$L\backslash p$	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

Table 7. Table 7: V-cycle with hybrid smoothing strategy for the Yeti footprint

	(a) Direct – Non-conforming							(b) PCG – Non-conforming
$L\backslash p$	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

Table 8. Table 8: V-cycle with hybrid smoothing strategy for the Fichera corner

	(a) PCG – Conforming					(b) PCG – Non-conforming
$L\backslash p$	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

Table 9. Table 9: V-cycle with hybrid smoothing strategy for the twisted Fichera corner

	(a) PCG – Conforming					(b) PCG – Non-conforming
$L\backslash p$	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

Equations

$$-\Delta \mathrm{u}=\mathrm{f}\text{mbox}in\Omega ,\mathrm{u}=\mathrm{g}\text{mbox}on\partial \Omega ,$$

$$(\nabla u,\nabla v)_{L_{2}(\Omega)}=(f,v)_{L_{2}(\Omega)}\qquad\mbox{for all $v\in H^{1}_{0}(\Omega)$}.$$

$$\overline{\Omega }=\bigcup _{k=1}^K\overline{{\Omega }_k}\text{mbox}with{\Omega }_k\cap {\Omega }_l=\varnothing \text{mbox}fork\ne l,$$

$$G_k:\hat{\Omega }:=(0,1{)}^d\to {\Omega }_k:=G_k(\hat{\Omega })\subset {\mathbb{R}}^d$$

$$V_{k,\ell }:=\{v\in L_2({\Omega }_k):v\circ G_k\in \underset{\delta =1}{\overset{d}{⨂}}{S}_{p,h_{\ell }}\}=\text{mbox}span\{{\phi }_{k,\ell }^{(i)}{\}}_{i=1}^{N_{k,\ell }},$$

$$V_{\ell }^c:=\{v\in V:v{|}_{{\Omega }_k}\in V_{k,\ell }\text{mbox}fork=1,\dots ,K\}.$$

$$(u_{\ell},v_{\ell})_{H^{1}(\Omega)}=(f,v_{\ell})_{L_{2}(\Omega)}\qquad\mbox{for all $v_{\ell}\in V_{\ell}^{c}$}.$$

$$A_{\ell }{\underline{u}}_{\ell }={\underline{f}}_{\ell },$$

$$V_{\ell }^n:=\{v\in L_2(\Omega ):v{|}_{{\Omega }_k}\in V_{k,\ell }\text{mbox}fork=1,\dots ,K\text{mbox}andv{|}_{\partial \Omega }=0\}.$$

$$\mathcal{N}:=\{(k,j):k<j\text{mbox}with{\Omega }_k\text{mbox}and{\Omega }_l\text{mbox}haveacommonedge.\}$$

$$(u_{\ell},v_{\ell})_{A_{\ell}}=(f,v_{\ell})_{L_{2}(\Omega)}\qquad\mbox{for all $v_{\ell}\in V_{\ell}^{n}$},$$

$$(u,v{)}_{A_{\ell }}$$

$$(u,v{)}_{K_{\ell }}$$

$$(u,v{)}_{B_{\ell }}$$

$$|u-u_L{|}_{Q_L}^2\le c{p}^2{h}_L^2|u{|}_{H^2(\Omega )}^2,$$

$$|u-u_L{|}_{Q_L}^2\le c(\log p{)}^4{h}_L^2|u{|}_{H^2(\Omega )}^2,$$

$$\frac{\sigma p^2}{h_L},$$

$$(u_{\ell },v_{\ell }{)}_{A_{\ell }}=(u_{\ell },v_{\ell }{)}_{A_{\ell +1}}\text{mbox}and(u_{\ell },v_{\ell }{)}_{Q_{\ell }}=(u_{\ell },v_{\ell }{)}_{Q_{\ell +1}}$$

$$\underset{{\tilde{\Sigma }}_{\ell }:=}{\underbrace{\frac{\sigma p^2}{h_L}}}=2^{L-\ell }\underset{{\Sigma }_{\ell }:=}{\underbrace{\frac{\sigma p^2}{h_{\ell }}}},$$

$$A_{\ell }{\underline{u}}_{\ell }={\underline{f}}_{\ell }.$$

$${\underline{u}}_{\ell }\leftarrow {\underline{u}}_{\ell }+L_{\ell }^{-1}\left({\underline{f}}_{\ell }-A_{\ell }{\underline{u}}_{\ell }\right),$$

$$(L_{\ell})_{i,j}=\left\{\begin{array}[]{ll}(A_{\ell})_{i,j}&\mbox{ if }i\geq j\\ 0&\mbox{ if }i<j\\ \end{array}\right..$$

$$\frac{\Vert A_L{\underline{u}}_L-{\underline{f}}_L{\Vert }_{{\ell }^2}}{\Vert {\underline{f}}_L{\Vert }_{{\ell }^2}}\le ϵ:=10^{-8}$$

$$L_{\ell ,T}:=P_T^{\top }{A}_{\ell }{P}_T,$$

$$L_{\ell }^{-1}:=\tau \sum _T{P}_{\ell ,T}{L}_{\ell ,T}^{-1}{P}_{\ell ,T}^{\top },$$

$$\tau \in (0,{\tau }^{\ast }),\delta \in (0,{\delta }^{\ast })\text{mbox}and{\nu }_{\ell }>{\nu }_{\ell }^{\ast }:=p\frac{{\tau }^{\ast }}{\tau }\frac{{\delta }^{\ast }}{\delta }\theta ,$$

$$\kappa ({\hat{A}}_{\ell }^{-1}{A}_{\ell })=\underset{v_{\ell }\in V_{\ell }}{\sup }\frac{|v_{\ell }{|}_{H^1(\Omega )}^2}{{\sum }_{k=1}^K|v_{\ell }\circ G_k{|}_{H^1(\hat{\Omega })}^2}\underset{v_{\ell }\in V_{\ell }}{\sup }\frac{{\sum }_{k=1}^K|v_{\ell }\circ G_k{|}_{H^1(\hat{\Omega })}^2}{|v_{\ell }{|}_{H^1(\Omega )}^2},$$

$$\kappa ({\hat{A}}_0^{-1}{A}_0)$$

$$\le \underset{v\in V}{\sup }\frac{|v{|}_{H^1(\Omega )}^2}{{\sum }_{k=1}^K|v\circ G_k{|}_{H^1(\hat{\Omega })}^2}\underset{v\in V}{\sup }\frac{{\sum }_{k=1}^K|v\circ G_k{|}_{H^1(\hat{\Omega })}^2}{|v{|}_{H^1(\Omega )}^2},$$

$$\tau \in (0,{\tau }^{\ast }),\delta \in (0,{\delta }^{\ast })\text{mbox}and{\nu }_{\ell }>{\nu }_{\ell }^{\ast }:=2^{L-\ell }(1+L-\ell {)}^{2}p(\log p{)}^4\frac{{\tau }^{\ast }}{\tau }\frac{{\delta }^{\ast }}{\delta }\theta ,$$

$${\tilde{W}}_{\ell }:=(I-(L_{\ell }^{GS}{)}^{-\top }{A}_{\ell })W_{\ell }(I-(L_{\ell }^{GS}{)}^{-1}{A}_{\ell }),$$

$$\Vert {\tilde{W}}_{\ell }{\Vert }_{A_{\ell }}\le \Vert I-(L_{\ell }^{GS}{)}^{-\top }{A}_{\ell }{\Vert }_{A_{\ell }}\Vert W_{\ell }{\Vert }_{A_{\ell }}\Vert I-(L_{\ell }^{GS}{)}^{-1}{A}_{\ell }{\Vert }_{A_{\ell }}\le \Vert W_{\ell }{\Vert }_{A_{\ell }}.$$

Full text

††journal: arXiv

1 Introduction

Isogeometric Analysis (IgA), see [18], 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 high-order method while the number of degrees of freedom behaves basically like for a low-order 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 $k$-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 (multi-patch 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 tensor-product B-splines 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 multi-patch 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. [18, 4, 5, 6, 31, 11, 24, 9, 32] and many others) and multigrid solvers (cf. [12, 10, 17, 16] and others) for such discretizations have been previously discussed. Since we are interested in $k$-refinement, we need results that are explicit in the spline degree. For the single-patch case, such error estimates have originally been given in [31] and later improved in [24]. In [16], a robust single-patch multigrid solver has been proposed and analyzed based on the error estimates from [31]. In [29], both the approximation error estimates and the multigrid solver have been extended to the conforming multi-patch 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 [2, 3]. Already in [20, 21], 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 [30]. 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 non-conforming discretizations are not known so far. In the present paper, it is shown how the multigrid solver from [29] can be extended to SIPG discretizations; we observe – as in [29] – 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 [29] to SIPG discretizations. It is worth noting that there are alternative solvers for multi-patch IgA, like FETI-type approaches, cf. [19, 15] and others, overlapping Schwarz type methods, cf. [7], or BDDC methods, cf. [8]; most of them, however, have not been worked out for the non-conforming 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 tensor-product 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 [25]. Here, the authors have improved their method by incorporating the geometry information into the preconditioner, cf. [23]. 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 Gauss-Seidel 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 [28].

In the present paper, we illustrate our findings with numerical experiments. All presented numerical experiments are available in the G+Smo library [22].

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 non-conforming 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 Gauss-Seidel 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 $\Omega\subset\mathbb{R}^{d}$ with $d\in\{2,3\}$ be an open and simply connected Lipschitz domain. Most of the numerical experiments are done for the two-dimensional domains shown in Figure 1.

The first domain is an L-shaped domain consisting of three quadrilaterals. Here, the geometry function is just the identity or a translation. On the coarsest grid level $\ell=0$, 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 $\ell=1,2,\ldots,$ are obtained by uniform refinement.

Consider the following Poisson model problem with Dirichlet boundary conditions. Find $\mathrm{u}\in H^{1}(\Omega)$ such that

[TABLE]

where $\mathrm{f}\in L_{2}(\Omega)$ and $\mathrm{g}\in H^{2}(\Omega)$ are given functions. Here and in what follows, $L_{2}(\Omega)$, $H^{r}(\Omega)$ and $H^{r}_{0}(\Omega)$ are the standard Lebesgue and Sobolev spaces.

The experiments are performed for the choice $\mathrm{g}(x,y):=\sin(\pi x)\,\sin(\pi y)$ and $\mathrm{f}:=-\Delta\mathrm{g}$; note that $\mathrm{g}$ is the exact solution of the problem.

After homogenization ($u:=\mathrm{u}-\mathrm{g}$, $f:=\mathrm{f}+\Delta\mathrm{g}$), the problem reads in variational form as follows. Find $u\in V:=H^{1}_{0}(\Omega)$ such that

[TABLE]

The computational domain $\Omega$ is a standard multi-patch domain. Thus, we assume that $\Omega$ is composed of $K$ non-overlapping patches $\Omega_{k}$:

[TABLE]

where each patch is represented by a sufficiently smooth bijective geometry function

[TABLE]

which can be continuously extended to $\overline{\widehat{\Omega}}$, the closure of $\widehat{\Omega}$. Moreover, we assume that the mesh introduced by the patches satisfies the following condition.

Assumption 1

For any $k\not=l$, the intersection $\overline{\Omega_{k}}\cap\overline{\Omega_{l}}$ is either (a) empty, (b) one common vertex, (c) the closure of one common edge, or – for $d=3$ – (d) the closure of one common face.

For each of the patches, we assume to have a hierarchy of grids with levels $\ell=0,1,\ldots,L$ obtained by uniform refinement, which we denote by

[TABLE]

where $\bigotimes_{\delta=1}^{d}S_{p,h_{\ell}}$ is the space of tensor-product splines of degree $p$, smoothness $H^{p}(\widehat{\Omega})$ (or, equivalently, $C^{p-1}(\widehat{\Omega})$) and grid size $h_{\ell}=2^{\ell}h_{0}$ on the parameter domain $\widehat{\Omega}$. Note that the grid can be non-uniform 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 quasi-uniform, i.e., there needs to be a constant $c>0$ such that all knot spans on grid level $\ell$ are bounded from below by $c\,h_{\ell}$. The functions $\varphi_{k,\ell}^{(i)}$ are assumed to form a (standard) B-spline or NURBS basis of $V_{k,\ell}$.

To be able to set up a conforming discretization, we need to assume that the function spaces are fully matching on the interfaces, cf. [29, Assumption 2.4]. For tensor-product B-spline 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

[TABLE]

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 patch-local 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 $u_{\ell}\in V_{\ell}^{c}$ such that

[TABLE]

Using the abovementioned basis for the space $V_{\ell}^{c}$, we obtain a standard matrix-vector problem: Find $\underline{u}_{\ell}\in\mathbb{R}^{N_{\ell}}$ such that

[TABLE]

where $A_{\ell}$ is the stiffness matrix, the vector $\underline{u}_{\ell}$ is the representation of $u_{\ell}$ with respect to the chosen basis and the load vector $\underline{f}_{\ell}$ is obtained by testing the function $f$ with the basis functions.

3 Symmetric interior penalty discontinuous Galerkin (SIPG) discretization

Following [20, 21, 30], we use a conforming isogeometric discretization for each patch and couple the patches using discontinuous Galerkin. We assume that the domain $\Omega$ is again subdivided into patches such that (3), (4) and Assumption 1 are satisfied. We assume again to have patch-local spaces $V_{k,\ell}$ as in (5), which are combined in a non-conforming (i.e., discontinuous) way, i.e., we just define

[TABLE]

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 $V_{\ell}^{n}\not\subset V$, 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

[TABLE]

to be the set of interface-indices. For each interface $I_{k,j}$ with $(k,l)\in\mathcal{N}$, we define the following symbols.

n is the outer normal vector of $\Omega_{k}$. (Thus, $-\textnormal{{n}}$ is the outer normal vector of $\Omega_{l}$.)

• 2.

$\llbracket\cdot\rrbracket$ is the jump operator: $\llbracket u\rrbracket:=u|_{\Omega_{k}}-u|_{\Omega_{l}}$.

• 3.

$\{\cdot\}$ is the averaging operator: $\{u\}:=\tfrac{1}{2}(u|_{\Omega_{k}}+u|_{\Omega_{l}})$.

Now, we can formulate the SIPG discretization as follows: Find $u_{\ell}\in V_{\ell}^{n}$ such that

[TABLE]

where we define

[TABLE]

There is some $\sigma_{0}>0$ independent of the grid size, the spline degree and the number of patches such that for all $\sigma\geq\sigma_{0}$, the bilinear form $(\cdot,\cdot)_{A_{\ell}}$ is coercive, cf. [30, Theorems 8 and 9]. Thus, for $\sigma\geq\sigma_{0}$, 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

[TABLE]

cf. [30]. By doing a more careful analysis, we obtain estimates of the form

[TABLE]

see [30, eq. (15)]. This significantly decreases the influence of the spline degree.

Note that the penalization term has the form

[TABLE]

i.e., it depends on the grid size on the finest grid $h_{L}$. This follows the ideas from [13]. The idea behind that is that

[TABLE]

hold, i.e., we obtain a multigrid solver with conforming coarse-grid correction. This means that – on the coarse grid levels – the discretization is over penalized by a factor of $2^{L-\ell}$, i.e.,

[TABLE]

where $\Sigma_{\ell}$ is the canonical parameter and $\widetilde{\Sigma}_{\ell}$ is the chosen one. We will see that this does not cause any problems for the examples we consider; following [13], convergence theory only holds if the number of smoothing steps is sufficiently increased for the coarser grid levels, cf. Remark 1.

Using a basis for the space $V_{\ell}^{n}$, we obtain a standard matrix-vector problem: Find $\underline{u}_{\ell}\in\mathbb{R}^{N_{\ell}}$ such that

[TABLE]

4 Multigrid solvers with Gauss-Seidel 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 non-conforming discretizations which are set up as discussed in the last two sections. As we have nested spaces in all cases, the matrix $I_{\ell-1}^{\ell}$ is always the canonical embedding from $V_{\ell-1}$ into $V_{\ell}$ and the restriction matrix $I_{\ell}^{\ell-1}$ is its transpose. The chosen method is presented as pseudo-code as Algorithm 1, where we choose $\mu=1$ for the V-cycle method or $\mu=2$ for the W-cycle method.

In the finite element world, Gauss-Seidel smoothers are known to be very efficient smoothers; thus, as a first attempt, we consider such a smoother. One forward Gauss-Seidel sweep can be represented by

[TABLE]

where $L_{\ell}$ is a lower-triangular matrix containing the coefficients of the stiffness matrix $A_{\ell}$, i.e., it is given by

[TABLE]

To be able to use our multigrid solver as preconditioner for a conjugate gradient solver, the post-smoothing procedure uses the transposed matrix $L_{\ell}^{\top}$, which represents a backward-Gauss-Seidel sweep.

All tables show the number of iterations required until the stopping criterion

[TABLE]

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 V-cycle and the W-cycle methods yield comparable iteration counts, we present the results for the more efficient V-cycle only. The number of smoothing steps is chosen as $\nu_{\ell}:=1$ in all cases.

The multigrid solver was implemented in C++ based on the G+Smo library [22]. The tables shown in the remainder of this section are obtained with the following command line instructions, where the values $L$ and $p$ are substituted accordingly.

⬇

git clone https://github.com/gismo/gismo.git

cd gismo

make

cd build/bin

./multiGrid_example -g domain2d/ldomain.xml -r $L$ -p $p$

  -s gs -i d                                         // Table [1](#S4.T1) (a)

./multiGrid_example -g domain2d/yeti_mp2.xml -r $L$ -p $p$

  -s gs -i d                                         // Table [2](#S4.T2) (a)

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.

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 V-cycle 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 non-trivial – the iteration counts are very similar to those of the L-shaped domain.

When one turns to the non-conforming discretizations, it immediately turns out that the multigrid solver utilizing the Gauss-Seidel 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. [12]. 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 [16]. That smoother requires that the spline space has a tensor-product 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 [29] based on a domain-decomposition approach. The key idea was to decompose all degrees of freedom on a per-piece bases. Pieces are the patch-interiors 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 piece-local smoothers, which we denote by $L_{\ell,T}^{-1}$, are defined as follows. For $T$ being a patch-interior, we choose $L_{\ell,T}^{-1}$ to be the subspace corrected mass smoother as proposed in [16]. We choose the scaling parameter (which was called $\sigma$ in [16]) to be $\delta^{-1}h_{\ell}^{-2}$ for some suitable chosen parameter $\delta>0$. If $T$ is an interface piece, we choose

[TABLE]

where the matrix $P_{\ell,T}$ represents the embedding of the piece $T$ in the whole space. The symbol $L_{\ell,T}^{-1}$ 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 piece-local smoothers, i.e., we choose

[TABLE]

where the sum is taken over all pieces $T$. Here, $\tau>0$ is some damping parameter to be chosen.

The convergence theory from [29] can be summarized by the following theorem.

Theorem 3

Assume that $\Omega\subset\mathbb{R}^{2}$ is such that full elliptic regularity is satisfied (cf. [29, Assumption 3.1]). Consider the conforming discretization and a multigrid solver with the smoother (14). There are constants $\tau^{*}$, $\delta^{*}$ and $\theta$ which are independent of $K$, $h$, $L$ and $p$ (but may depend particularly on the geometry functions and the maximum number of neighbors of a patch) such that for

[TABLE]

the W-cycle multigrid method converges with a convergence rate $q\leq\max_{\ell}\nu_{\ell}^{*}/\nu_{\ell}$.

Note that the terms $\frac{\tau^{*}}{\tau}$ and $\frac{\delta^{*}}{\delta}$ 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. [29, Assumption 3.1]). Thus, it is applicable to the Yeti footprint but it is not directly applicable to the L-shaped 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 $p$. As already outlined in [29], this seems to be too pessimistic since the numerical experiments have shown that $\nu_{\ell}:=1$ 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 V-cycle multigrid method with 1+1 smoothing steps of the proposed smoother $L_{\ell}$ (on all grid levels). The damping parameter $\tau$ is chosen as indicated with the option --MG.Damping and the scaling parameter $\delta$ 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 $L$ and $p$ are substituted accordingly:

⬇

./multiGrid_example -g domain2d/ldomain.xml -r $L$ -p $p$

  -s scms --MG.Damping 1 --MG.Scaling .12 -i d        // Table [3](#S5.T3) (a)

./multiGrid_example -g domain2d/yeti_mp2.xml -r $L$ -p $p$

  -s scms --MG.Damping .25 --MG.Scaling .2 -i cg     // Table [5](#S5.T5) (a)

The results for the PCG experiments presented Table 3 (b) are obtained by replacing the option -i d by the option -i cg.

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 L-shaped 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 pre-asymptotic. The reason for this is that on coarser grid levels, the geometry is not resolved exactly. Let $A_{\ell}$ be the original stiffness matrix and $\widehat{A}_{\ell}$ be the simplified stiffness matrix obtained by neglecting the geometry function. Then, we have

[TABLE]

which yields

[TABLE]

which can be finally bounded by a constant times some power of the quantity $\|\nabla G\|_{L_{\infty}(\widehat{\Omega})}\|(\nabla G)^{-1}\|_{L_{\infty}(\widehat{\Omega})}$. 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 [27, Section 7.4].

As a next step, we turn towards the non-conforming discretizations. Here, each degree of freedom is assigned to exactly one patch. So, it would be tempting to set up a patch-wise splitting of the degrees of freedom. Unfortunately, numerical experiments have shown that this approach does not work well. So, we follow the approach from [29] also in the non-conforming case and split the degrees of freedom again into pieces $T$. 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 $\Omega\subset\mathbb{R}^{2}$ is such that full elliptic regularity holds (cf. [30, Assumption 4]) and assume that the geometry functions (but not necessarily the discretizations) agree on the interfaces (cf. [30, Assumption 2]). Consider the SIPG discretization and a multigrid solver with the smoother (14). There are constants $\tau^{*}$, $\delta^{*}$ and $\theta$ which are independent of $K$, $h$, $L$ and $p$ (but may depend particularly on the geometry functions and the maximum number of neighbors of a patch) such that for

[TABLE]

the W-cycle multigrid method converges with a convergence rate $q\leq\max_{\ell}\nu_{\ell}^{*}/\nu_{\ell}$.

We give the proof of this theorem in the Appendix; the proof is based on the error estimates from [30].

One might observe that the number of smoothing steps required by this convergence theorem increases like $(1+L-\ell)2^{L-\ell}$. This follows the approach suggested in [13] and is related to the chosen over-penalization discussed in Section 3.

Remark 1

Note that the number of degrees of freedom on the coarser grid levels is smaller by a factor of $2^{d(L-\ell)}$. 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) $d\geq 3$ or (b) the V-cycle is considered. If we consider $d=2$ and the W-cycle, the choice (17) yields that the computational complexity grows like $N_{L}L^{3}$, where $N_{L}$ 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 $p$ is too pessimistic; thus, we again choose $\nu_{\ell}:=1$ 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 $V_{\ell}^{c}\subset V_{\ell}^{n}$ and the SIPG formulation converges to the conforming discretization for $\sigma\rightarrow\infty$. 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 $S_{p,h_{\ell}}(\widehat{\Omega})$. For one third of the patches, we use the spline space $S_{p+1,2h_{\ell}}(\widehat{\Omega})$. For the last third of the patches, we use the spline space $S_{p,2h_{\ell}}(\widehat{\Omega})$. 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 non-conforming case shown in this section are obtained with the following code, where the values $L$ and $p$ are substituted accordingly:

⬇

./multiGrid_example -g domain2d/ldomain.xml -r $L$ -p $p$ --DG

  --NonMatching -s scms --MG.Damping .9 --MG.Scaling .12

  -i d                                            // Table [4](#S5.T4) (a)

./multiGrid_example -g domain2d/yeti_mp2.xml -r $L$ -p $p$ --DG

  --NonMatching -s scms --MG.Damping .25 --MG.Scaling .2

  -i cg                                           // Table [5](#S5.T5) (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.

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 non-simple geometry functions since it is based on the close connection between the stiffness matrix $A_{\ell}$ and the simplified stiffness matrix $\widehat{A}_{\ell}$. The results for the Gauss-Seidel 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 Gauss-Seidel sweep followed by the subspace corrected mass smoother for pre-smoothing and the subspace corrected mass smoother followed by one backward Gauss-Seidel sweep for post-smoothing. 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

[TABLE]

where $W_{\ell}$ is the iteration matrix of the (V or W cycle, respectively) multigrid method with the subspace corrected mass smoother. Since the Gauss-Seidel iteration is stable in the energy norm, we obtain

[TABLE]

So, we have using the results from the last section the convergence of the W-cycle multigrid method with hybrid smoother. Thus, we obtain as follows.

Corollary 5

Consider the multigrid solver with the hybrid smoother. Under the assumptions of Theorem 3 or 4, respectively, the W-cycle multigrid method converges with a convergence rate $q\leq\max_{\ell}\nu_{\ell}^{*}/\nu_{\ell}$.

The tables for the experiments with the hybrid smoother are obtained with the following code, where the values $L$ and $p$ are substituted accordingly:

⬇

./multiGrid_example -g domain2d/yeti_mp2.xml -r $L$ -p $p$

  -s hyb --MG.Damping .25 --MG.Scaling .1 -i d       // Table [6](#S6.T6) (a)

./multiGrid_example -g domain2d/yeti_mp2.xml -r $L$ -p $p$ --DG

  --NonMatching -s hyb --MG.Damping .25 --MG.Scaling .1

  -i d                                               // Table [7](#S6.T7) (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 $p$, the Gauss-Seidel 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 Gauss-Seidel 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 $\nu_{\ell}:=1$ means basically that $2$ pre- and $2$ post-smoothing steps are applied. Since the Gauss-Seidel 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 two-dimensional 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 non-trivial geometry function, which we call twisted Fichera corner.

The tables for the three dimensional domains are obtained with the following code, where the values $L$ and $p$ are substituted accordingly:

⬇

./multiGrid_example -g domain2d/fichera.xml -r $L$ -p $p$

  -s hyb --MG.Scaling .12 --MG.Damping 1 -i cg        // Table [8](#S6.T8) (a)

./multiGrid_example -g domain2d/twisted_fichera.xml -r $L$ -p $p$

  -s hyb --MG.Scaling .12 --MG.Damping .25 -i cg     // Table [9](#S6.T9) (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.

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 multi-patch IgA with conforming and non-conforming 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 tensor-product 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 Gauss-Seidel 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 $c$ 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. [30, Assumption 3]), the number of neighbors of a patch (cf. [29, Assumption 2.3]), the constant in the elliptic regularity assumption (cf. [29, Assumption 3.1]) and the quasi-uniformity of the grid, i.e., the ratio between the largest and the smallest knot span of the knot vectors on one level. We write $A\lesssim B$ if and only if there is a constant $c$ such that $A\leq c\;B$ and we write $A\eqsim B$ if and only if $A\lesssim B$ and $B\lesssim A$.

First, we show the following lemma, which is basically a trace inequality.

Lemma 6

Let $S:=S^{(1)}\otimes S^{(2)}$ be the space of tensor-product splines of degree $p$ on a quasi-uniform grid with size $h$ on the parameter domain $\widehat{\Omega}:=(0,1)^{2}$. Then, the estimate

[TABLE]

holds for all $u\in S$ and all $\theta\geq 1$.

Proof 1

For $\nu=1,2$, let $(\psi_{\nu,i})_{i=1}^{N_{\nu}}$ be the eigenfunctions of $S^{(\nu)}$, i.e., such that

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta and $\lambda_{\nu,1}\leq\lambda_{\nu,2}\leq\cdots\leq\lambda_{\nu,N_{\nu}}$ are the corresponding eigenvalues. Using coercivity of $(\cdot^{\prime},\cdot^{\prime})_{L_{2}(0,1)}$ and a standard inverse estimate, cf. [26, Corollary 3.94], we obtain

[TABLE]

We define level sets

[TABLE]

for $m\in\{1,2,3,\ldots,M\}$, where

[TABLE]

is the number of level sets. Note that by construction every eigenvalue belongs to exactly one level set. Every function $u\in S$ can be represented as

[TABLE]

A standard trace estimate, cf. [29, Lemma 4.4], yields

[TABLE]

where

[TABLE]

for $a,b\in\mathbb{N}_{0}$. Since $\theta\geq 1$ and since all eigenvalues are in $I_{1,m}$ or $I_{2,n}$, respectively, we obtain

[TABLE]

Using (18), we obtain further

[TABLE]

Finally, the Cauchy-Schwarz inequality and orthogonality of the basis functions (both in $L_{2}$ and $H^{1}$) yield

[TABLE]

which finishes the proof. \qed

Now, we give bounds on the smoother which allow to show the smoothing property.

Lemma 7

Provided the assumptions of Theorem 4, the estimate

[TABLE]

holds, where $\widetilde{L}_{\ell}:=Q_{\ell}+(1+2^{\ell-L}h_{\ell}^{-2})M_{\ell}$ and $M_{\ell}$ is the standard mass matrix.

Proof 2

The proof of this Lemma requires the notation from [29], i.e., we denote the set of all patch-interiors by $\mathbb{K}$, the set of all edges by $\mathbb{E}$ and the set of all vertices by $\mathbb{V}$.

Observe that we have

[TABLE]

for all $\ell=0,2,\ldots,L$, where $Q_{\ell}$ is defined by (11) and (12). For $\ell=L$, this statement directly follows from [30, Theorem 8]. Since [30, Theorem 8] also holds in cases of over-penalization, we can apply that theorem also to the case $\ell<L$ and obtain (19) also in that case.

As first step, we bound $L_{\ell}$ from below. (The following arguments are analogous to [29, Lemma 4.3].) The triangle inequality yields

[TABLE]

For $T\in\mathbb{K}$, [16, Lemma 8] and (19) yield $L_{\ell,T}\gtrsim P_{\ell,T}^{\top}Q_{\ell}P_{\ell,T}^{\top}\gtrsim P_{\ell,T}^{\top}A_{\ell}P_{\ell,T}^{\top}$. For $T\in\mathbb{E}\cup\mathbb{V}$, we have by definition $L_{\ell,T}=P_{\ell,T}^{\top}A_{\ell}P_{\ell,T}^{\top}$. Thus, we obtain from (20)

[TABLE]

and for all $\tau\in(0,\tau^{*})$ with $\tau^{*}$ small enough further

[TABLE]

which shows the first part of the desired inequality.

Now, we bound $L_{\ell}$ from above. We use the decomposition

[TABLE]

cf. (11). Using (19), we obtain

[TABLE]

for all $T\in\mathbb{E}\cup\mathbb{V}$ and, therefore,

[TABLE]

where

[TABLE]

and

[TABLE]

Completely analogous to [29, Lemma 4.7], we obtain

[TABLE]

So, it remains to bound $\widetilde{J}_{\ell}$ from above. Since the restriction of $J_{\ell}$ to any patch-interior vanishes, the same arguments as in the proof of [29, Lemma 4.7] and the triangle inequality yield

[TABLE]

and therefore

[TABLE]

Note that $J_{\ell}$ models jumps and note that these jumps can be bounded from above using the triangle inequality with the function values on both sides. Thus, we obtain

[TABLE]

where $u_{\ell}|_{\Omega_{k}}$ is the restriction of $u_{\ell}$ to the patch $\Omega_{k}$ and $u_{\ell}|_{\Omega_{k}}|_{T}$ is the evaluation of the continuous extension of that function to the vertex $T$ at that vertex and $\psi(x)=\max\{1-x/h_{\ell},0\}^{p}$ is the corresponding basis function. Using $\|\psi\|_{L_{2}(0,1)}^{2}\eqsim p^{-1}h_{\ell}$, cf. [29, Eq. (4.16)], Lemma 6 (with $\theta:=(1+h_{\ell}^{-2}2^{\ell-L})^{1/2}$), and $h_{\ell}\leq 1$, we further obtain

[TABLE]

Using [30, Lemma 6], we obtain

[TABLE]

and thus

[TABLE]

This shows together with (23) and $h_{L}\eqsim 2^{\ell-L}h_{\ell}$

[TABLE]

Since $\delta\in(0,\delta^{*})$ and since $\tau^{*}\eqsim 1$, we obtain

[TABLE]

which finishes together with (21) and (22) the proof. \qed

Lemma 8

Let $\|v\|_{Q_{\ell}^{+}}^{2}:=|v|_{H^{1}(\Omega)}^{2}+\sigma^{-2}p^{-4}4^{\ell-L}h_{\ell}^{2}|v|_{H^{2}(\Omega)}^{2}$. The estimate

[TABLE]

holds for all $u\in H^{2}(\Omega)$.

Proof 3

Let

[TABLE]

be the set of all globally continuous functions which are linear locally. Observe that $W\subseteq V_{\ell}^{n}$. Using $u$ and $w$ being continuous, we obtain

[TABLE]

For the choice $w\in H^{1}(\Omega)$ with $w|_{\Omega_{k}}:=w_{k}=\widehat{w}_{k}\circ G_{k}^{-1}$, where

[TABLE]

we further obtain using standard approximation error estimates and [30, Lemma 6]

[TABLE]

Using $h_{\ell}\leq 1$, $\sigma\geq 1$, $p\geq 2$, and $L\geq\ell$, we obtain the desired result.

Lemma 9

Provided the assumptions of Theorem 4, the estimate

[TABLE]

holds, where $\widetilde{L}_{\ell}$ is as in Lemma 7.

Proof 4

Let $u_{\ell}\in V_{\ell}^{n}$ be arbitrary but fixed. Let $f_{\ell}\in V_{\ell}^{n}$ be such that

[TABLE]

Let $u_{\ell-1}\in V_{\ell-1}$ and $u\in V$ be such that

[TABLE]

Using $f_{\ell}\in L_{2}(\Omega)$ and full elliptic regularity, cf. [29, Assumption 3.1], we obtain $u\in H^{2}(\Omega)$ and

[TABLE]

[30, Theorems 12 and 13]** and Lemma 8 yield

[TABLE]

where $\sigma_{\ell}=2^{L-\ell}p^{2}\sigma.$ Using the triangle inequality, $\sigma\eqsim 1$, $\log(ab)\lesssim\log a\log b$, $h_{\ell}\eqsim h_{\ell-1}$, $p\geq 2$, and we obtain further

[TABLE]

Using (19), (24) and the definition of $u_{\ell-1}$, we obtain further

[TABLE]

This yields

[TABLE]

and thus

[TABLE]

Using this estimate and the stability of the $A_{\ell}$-orthogonal projection, we obtain

[TABLE]

and further

[TABLE]

Using the identity $\|A^{\top}A\|\leq\|A\|^{2}$, we finally obtain the desired result. \qed

Finally, we can show Theorem 4. Here, we follow the classical approach as introduced by Hackbusch, cf. [14].

Proof 5** (of Theorem 4)**

Lemma 7 yields $A_{\ell}\leq L_{\ell}$. Using standard arguments, cf. [17, Lemma 2] or [14], the smoothing property

[TABLE]

follows. Using Lemma 7, we obtain further

[TABLE]

which shows together Lemma 9

[TABLE]

This statement shows convergence of the two-grid method if $\nu$ is large enough. Standard arguments, cf. [14], allow to extend the analysis to the W-cycle multigrid method.

Acknowledgments

The author was supported by the Austrian Science Fund (FWF): grant P31048, and by the bilateral project DNTS-Austria 01/3/2017 (WTZ BG 03/2017), funded by Bulgarian National Science Fund and OeAD (Austria).