# Multigrid on unstructured meshes with regions of low quality cells

Yuxuan Chen<EMAIL_ADDRESS>https://orcid.org/0000-0002-4152-2011 School
of Mathematical Sciences, Soochow University, Suzhou, China  Garth N. Wells
<EMAIL_ADDRESS>https://orcid.org/0000-0001-5291-7951 Department of
Engineering, University of Cambridge, Cambridge, United Kingdom

###### Abstract

The convergence of multigrid methods degrades significantly if a small number
of low quality cells are present in a finite element mesh, and this can be a
barrier to the efficient and robust application of multigrid on complicated
geometric domains. The degraded performance is observed also if intermediate
levels in a non-nested geometric multigrid problem have low quality cells,
even when the fine grid is high quality. It is demonstrated for geometric
multigrid methods that the poor convergence is due to the local failure of
smoothers to eliminate parts of error around cells of low quality. To overcome
this, a global–local combined smoother is developed to maintain effective
relaxation in the presence of a small number of poor quality cells. The
smoother involves the application of a standard smoother on the whole domain,
followed by local corrections for small subdomains with low quality cells.
Two- and three-dimensional numerical experiments demonstrate that the degraded
convergence of multigrid for low quality meshes can be restored to the high
quality mesh reference case using the proposed smoother. The effect is
particularly pronounced for higher-order finite elements. The results provide
a basis for developing efficient, non-nested geometric multigrid methods for
complicated engineering geometries.
Keywords: multigrid, multigrid smoother, finite element methods, cell quality,
domain decomposition.

##  1 Introduction

Multigrid methods have the potential to be optimal solvers ($O(n)$ complexity,
where $n$ is the number of unknowns) for systems arising in the solution of
elliptic partial differential equations [6, 15, 37], and lend themselves to
efficient parallel implementations. This makes tractable high fidelity
simulations of complex engineering components and enables the computation of
engineering problems at a system level.

Finite element simulations on complex geometries are invariably performed on
unstructured grids, and cell quality will vary with position. This is
particularly the case when representing geometrically complex shapes with
modest cell counts. It is, however, well recognised that cell quality can have
a significant impact on the performance of iterative solvers [33, 13, 19].
When disappointing performance of multigrid solvers is observed in engineering
practice (slow performance with high iteration counts, or failure to
converge), low cell quality is a common cause [31]. Even very small regions
with just a few low quality cells can cause convergence of a multigrid
preconditioned iterative solver to stall. On the other hand, it has been shown
that if using linear Lagrange finite elements for solving Poisson equation,
the standard _a priori_ estimate may remain valid on a mesh with low quality
regions [10], but the associated finite element linear system is poorly
conditioned.

Ideally, computational grids would be of sufficient quality to not degrade
performance of the linear solver. However, we consider two contexts in which
this is not universally feasible. The first is that the generation of meshes
for highly complex geometries and for which ensuring that all cells are of
high quality may be very difficult. Considerable research has been devoted to
improving mesh quality, e.g. [21, 22, 20]. However, even in cases where it is
technically possible to create a high quality mesh, from a workflow
perspective it may be more efficient (less costly in overall time) if the
solver can deal sub-optimal grids efficiently. The second case is geometric
(non-nested) multigrid, in which the quality of intermediate coarse grids may
be compromised by under-resolution with respect to the domain shape
complexity, i.e. a low-resolution triangulation of a complex domain will
necessarily compromise on cell quality. Ideally the fine grid will be of high
quality to provide good approximation properties, but the coarse grids are
primarily a vehicle for constructing a (hopefully) fast solver. Robust
performance of a geometric multigrid solver with respect to the quality of
intermediate grids is appealing for complex engineering geometries.

We examine the performance of a geometric multigrid method for the finite
element method on non-nested unstructured meshes in the presence of a small
number of low quality cells. It is shown, through examples, that standard
multigrid smoothers fail to eliminate parts of the error in highly localised
regions around small clusters of low quality cells. The component of the
residual that is not reduced by the smoother can be ‘lost’ in the restriction
operation. Building on the observation of how smoothers perform around low
quality cells, we construct a global–local combined smoother in which (i) a
standard multigrid smoother is applied over the entire grid, followed by (ii)
a local correction on small regions with low cell quality using a direct
solver. The combined smoother is effectively a Schwarz-type domain
decomposition method with full overlap [34, 9]. We assume that the number of
unknowns in each low quality region is small such that the cost of applying a
direct solver on these regions is small relative to other operations in the
solver. In this work, we assume conforming boundaries between geometric
multigrid levels. Examples of handling non-conforming boundaries in geometric
multigrid include [28, 35, 8].

The remainder of this paper is structured as follows. In section 2, we give a
brief overview of multigrid. An explanation of why the performance of
multigrid degrades in the presence of low quality cells is presented in
section 3, supported by numerical examples. A global–local combined smoother
is introduced in section 4, followed by a brief detour in section 5 to discuss
a specific issue for Chebyshev smoothers. Numerical examples are presented in
section 6 for Poisson and elasticity problems. Conclusions are drawn in
section 7.

##  2 Multigrid background

Consider a domain $\Omega\subset\mathbb{R}^{\dim}$, where $\dim=1,2,3$, and
triangulations (grids) of the domain $\Omega_{l}$, $l=1,2,\ldots,L$. With
increasing index $l$ the grids become coarser, i.e. the number of cells
reduces. We wish to solve a finite element problem on the finest grid
$\Omega_{1}$. The finite element method generates the discrete operator
$A_{1}\in\mathbb{R}^{n\times n}$, and the task is to solve the linear system

$A_{1}u_{1}=b_{1},$ (1)

where $b_{1}\in\mathbb{R}^{n}$ is the right-hand side vector and
$u_{1}\in\mathbb{R}^{n}$ is the vector of degrees-of-freedom.

If $V_{l}$ is a finite element space on $\Omega_{l}$, in the case of nested
grids and using the same element type on each level we have $V_{l+1}\subset
V_{l}$. With a view to complex geometries, we consider in this work non-nested
meshes, i.e. $V_{l+1}\not\subset V_{l}$. A _prolongation_ operator
$P_{l}:V_{l+1}\rightarrow V_{l}$ projects a finite element function on grid
$\Omega_{l+1}$ onto the next finest grid $\Omega_{l}$. We consider continuous
Lagrange finite element spaces and define the prolongation operator using
interpolation, in which case the components of $P_{l}$ are given by

$[P_{l}]_{ij}=\varphi^{(l+1)}_{j}(x_{i}),$ (2)

where $\varphi^{(l+1)}_{j}$ is the basis function associated with the $j$th
degree-of-freedom on level $l+1$ and $x_{i}$ is the interpolation coordinate
for the $i$th basis function on the finer $\Omega_{l}$ grid. The _restriction_
operator $R_{l}$ maps a function on grid $l$ to the coarser $l+1$ grid.
Following the Galerkin approach, we set $R_{l}=P_{l}^{T}$, leading to

$A_{l+1}:=P_{l}^{T}A_{l}P_{l}.$ (3)

for the operator on level $l+1$, $A_{l+1}$, and

$b_{l+1}:=P_{l}^{T}\left(b_{l}-A_{l}u_{l}\right)$ (4)

for the residual vector on level $l+1$, $b_{l+1}$.

A _smoother_ $S$ provides an approximate solution $u_{l}$ to the problem
$A_{l}u_{l}=b_{l}$:

$u_{l}\leftarrow S^{\nu}(A_{l},b_{l},u_{l}),$ (5)

where $\nu$ denotes the number of applications of the smoother. Commonly used
smoothers include (weighted) Jacobi, (symmetric) Gauss–Seidel and Chebyshev
iterations. For the coarsest grid, $\Omega_{L}$, a direct solver is employed.
Given the fine grid operator $A_{1}$, the fine grid right-hand side vector
$b_{1}$ and the prolongation operators, a V-cycle Galerkin multigrid algorithm
with $L$ levels for solving eq. 1 is summarised in algorithm 1.

Algorithm 1 Multigrid V-cycle of $L$ levels to solve $A_{1}u_{1}=b_{1}$.

1:procedure $u_{l}\leftarrow\textrm{Vcycle}(A_{l},b_{l},u_{l},l,\nu)$

2: for $l=1,2,\cdots,L$ do

3: if $l\neq L$ then

4: Pre-smoothing $u_{l}\leftarrow S^{\nu}(A_{l},b_{l},u_{l})$.

5: Coarse grid construction $A_{l+1}=P_{l}^{T}A_{l}P_{l}$,
$b_{l+1}=P_{l}^{T}(b_{l}-A_{l}u_{l})$ and $u_{l+1}=P_{l}^{T}u_{l}$.

6: else

7: Direct solver on the coarsest grid for $A_{L}u_{L}=b_{L}$.

8: for $l=L-1,L-2,\cdots,1$ do

9: Updating current solution $u_{l}\leftarrow u_{l}+P_{l}u_{l+1}$.

10: Post-smoothing $u_{l}\leftarrow S^{\nu}(A_{l},b_{l},u_{l})$.

We consider a hierarchy of non-nested grids, but restricted to the case where
all grids conform to the same boundary. A natural extension would be to non-
conforming boundaries, which would require some additional considerations [28,
35, 8].

##  3 Why multigrid converges slowly with low quality meshes

We explore how multigrid behaves with non-nested unstructured grids when
levels have small regions of low quality cells. There is no universal measure
of cell quality, and no sharp distinction between ‘low’ and ‘high’ quality
cells. Some cell quality measures are discussed in [30, 23]. We consider
simplex cells and select the radius ratio as a measure of cell quality; the
radius ratio is the ratio of a cell’s inscribed sphere radius ($R_{I}$) to its
circumscribing sphere radius ($R_{C}$). We define a normalised radius ratio
$\gamma$ as:

$\gamma:=\frac{R_{I}}{\gamma^{*}R_{C}},$ (6)

where $\gamma^{*}$ is the optimal radius ratio, which is equal to the
geometric dimension. The normalised radius ratio is in the range $(0,1]$, with
ideal cells having a measure of $1$, and degenerate cells having the measure
$0$.

###  3.1 Model problem

We solve the homogeneous Poisson problem on the unit square
$\Omega=[0,1]^{2}$,

$\displaystyle-\nabla^{2}u$ $\displaystyle=0\quad\text{in}\ \Omega,$ (7)
$\displaystyle u$ $\displaystyle=0\quad\text{on}\ \partial\Omega,$

using the finite element method with degree one Lagrange basis functions and a
two-level multigrid solver (V-cycle) with unstructured, non-nested meshes. The
coarse grid contains $68$ cells and the fine grid contains $272$ cells. The
initial solution guess interpolates $u^{(0)}=\sin(10\pi x)\sin(10\pi y)$.

Two cases are considered: the first uses a high quality fine grid, and the
second uses a perturbation of the high quality fine grid to create a small
region of low quality cells. The perturbed fine grid is shown in fig. 1. The
minimum angle over all cells in the perturbed fine grid is approximately
$\pi/180$, corresponding to a normalised radius ratio of $10^{-3}$ for the
cell. There are three ‘low’ quality cells with a normalised radius ratio of
less than $0.1$.

Figure 1: Unstructured grid of the unit square with a poor quality region near
the centre.

For a symmetric Gauss–Seidel smoother (one iteration at each application of
the smoother, pre- and post-smoothing), the relative residual after each cycle
is recorded and presented in fig. 2. The convergence rate for the low quality
fine grid case is dramatically slower than the high quality case.

Figure 2: Relative residual at the end of each multigrid cycle for high and
low quality meshes. Symmetric Gauss–Seidel is used as the smoother.

###  3.2 Local failure of smoothers

We investigate the performance of a smoother only for the model problem on the
grid shown in fig. 1. The interpolation of the initial guess $u^{(0)}$ is
shown in fig. 3(a), and the absolute value of the error after five symmetric
Gauss–Seidel iterations is shown in fig. 3(b). After applications of the
smoother, in most of the domain the oscillatory error has been eliminated.
However, a localised error persists around the low quality cells.

(a) Initial guess $u^{(0)}$.

(b) Absolute value of error on each vertex.

Figure 3: A (a) high frequency Fourier mode as initial guess, and (b) absolute
value of the error at each vertex after five iterations of symmetric
Gauss–Seidel.

We also apply a Jacobi–preconditioned Chebyshev smoother [1] to the model
problem. The largest eigenvalue used in the Chebyshev method is estimated by
the Krylov–Schur method [36, 18] (tolerance set as $10^{-8}$), and the
smallest eigenvalue is set as one tenth of the largest one. The absolute value
of the error at each vertex after five applications of the smoother is shown
in fig. 4. As with the Gauss–Seidel smoother, the localised error persists in
the region of low cell quality.

Figure 4: Absolute value of error on each vertex after five iterations of
Jacobi-preconditioned Chebyshev method on the low quality unit square model
problem.

The observation from the model problem is that for a mesh with a small region
of low quality cells, standard multigrid smoothers are ineffective in the
region around the low quality cells. The error is not eliminated and may
subsequently be lost in the restriction operator to the coarser grid, and this
slows the rate of convergence of multigrid. The localised nature of the
residual error motivates the consideration of local smoothers/solvers that can
provide effective smoothing in small regions.

##  4 Global–local correction smoother

We consider a global–local combined smoother that overcomes the issues
demonstrated in the previous section. We consider the linear system $Au=b$,
where $A\in\mathbb{R}^{n\times n}$. Let $u^{(k)}$ be the approximate solution
obtained after $k$ iterations of some iterative method, and
$r^{(k)}:=b-Au^{(k)}$ be the corresponding residual.

###  4.1 Local residual correction

Suppose there are $D$ small, closed subdomains containing low quality cells,
which we denote by $\Omega_{B}^{d}\subset\Omega,\ d=1,2,\ldots,D$. The
subdomains are defined to be disjoint,
$\bigcap_{d}\bar{\Omega}_{B}^{d}=\emptyset$. Let $B_{d}$ be the set of all
degrees-of-freedom (DOFs) $\beta_{i}$ associated with the closure of the
subdomain $\bar{\Omega}_{B}^{d}$, i.e.,
$B_{d}=\\{\beta_{i}:i=1,2,\ldots,n_{B}^{d}\\}$, where $n_{B}^{d}:=\lvert
B_{d}\rvert$. Let $I_{d}\in\mathbb{R}^{n\times n_{B}^{d}}$ be the natural
inclusion mapping degrees-of-freedom in $\bar{\Omega}_{B}^{d}$ to the the
whole domain $\Omega$, given by:

$[I_{d}]_{ij}:=\begin{cases}1\quad\quad i=\beta_{j},\,\beta_{j}\in B_{d},\\\
0\quad\quad\text{otherwise}.\end{cases}$ (8)

Consider a restriction of the matrix $A$ to subdomain $\bar{\Omega}_{B}^{d}$,

$A^{[d]}:=I_{d}^{T}AI_{d},$ (9)

which has size $n_{B}^{d}\times n_{B}^{d}$. Similarly for residual,

$r^{[d]}:=I_{d}^{T}r^{(k)}.$ (10)

A local residual correction system on subdomain $\bar{\Omega}_{B}^{d}$ is then
given by:

$A^{[d]}e^{[d]}=r^{[d]},$ (11)

where $e^{[d]}$ is the local error correction. If the original system $A$ is
symmetric positive-definite, then the local system $A^{[d]}$ is symmetric
positive-definite [38, Lemma 3.1]. We assume that the local correction systems
are small and can be solved efficiently using a direct solver:

$e^{[d]}=(I_{d}^{T}AI_{d})^{-1}I_{d}^{T}r^{(k)}.$ (12)

The local error correction can be mapped back to the global domain by applying
$I_{d}$ and the approximate solution $u^{(k)}$ corrected,

$u^{(k+1)}=u^{(k)}+\sum\limits_{d=1}^{D}I_{d}e^{[d]}.$ (13)

The local residual correction can be written in a preconditioner form as

$u^{(k+1)}=u^{(k)}+S_{c}\left(b-Au^{(k)}\right),$ (14)

where

$S_{c}:=\sum\limits_{d=1}^{D}I_{d}(I_{d}^{T}AI_{d})^{-1}I_{d}^{T}.$ (15)

In the case of an interior region of low quality cells, the local correction
is essentially the solution of a local Dirichlet problem.

The local residual correction procedure is summarised in algorithm 2.

Algorithm 2 Local residual correction $S_{c}$

procedure $u^{(k+1)}=S_{c}(A,b,u^{(k)})$

Identify low quality regions $\Omega_{B}^{d}$ and the corresponding DOF sets
$B_{d}$, $d=1,2,\cdots,D$.

for $d=1,2,\cdots,D$ do

Construct the local residual correction systems $A^{[d]}$, $r^{[d]}$ via eq. 9
and eq. 10.

Solve the local residual correction system $A^{[d]}e^{[d]}=r^{[d]}$ by a
direct method.

Correct $u^{(k)}$ by adding local errors via eq. 13.

It is straightforward to implement and the additional computational cost is
small if the poor quality regions are few and small. It is noted that the
local residual correction is similar to subspace correction methods, including
the parallel subspace correction and the successive subspace correction,
proposed by Xu [39, 38]. In particular, the low quality regions can be viewed
as (local) subdomains. The local correction $S_{c}$ corrects the error on the
local subdomains, and can be interpreted as a parallel subspace correction on
the subdomains with poor quality cells.

###  4.2 A global–local combined smoother for multigrid

We consider a global–local combined multigrid smoother, with the local
smoother applied in a number of small regions to overcome the deleterious
effect of any low quality cells. The smoother has a ‘sandwich’ form, with the
local correction smoother $S_{c}$ applied on subdomains $\Omega_{B}^{d}$,
followed by a standard global smoother $S_{g}$, e.g. a symmetric Gauss–Seidel,
on the whole domain/level, followed by another application of the local
smoother $S_{c}$. The smoother involves three steps:

$\displaystyle u^{(k+1/3)}$
$\displaystyle=u^{(k)}+S_{c}\left(b-Au^{(k)}\right),$ (16) $\displaystyle
u^{(k+2/3)}$ $\displaystyle=u^{(k+1/3)}+S_{g}\left(b-Au^{(k+1/3)}\right),$
$\displaystyle u^{(k+1)}$
$\displaystyle=u^{(k+2/3)}+S_{c}\left(b-Au^{(k+2/3)}\right).$

Rearranging,

$u^{(k+1)}=u^{(k)}+S_{gc}\left(b-Au^{(k)}\right),$ (17)

where

$S_{gc}=2S_{c}-S_{c}AS_{c}+(I-S_{c}A)S_{g}(I-AS_{c}).$ (18)

demonstrating that the smoother is symmetric.

###  4.3 Relationship to a domain decomposition method

The global–local combined smoother can be viewed as a domain decomposition by
considering two domains: the whole domain and several small subdomains, as
shown in fig. 5. Taking the whole domain as the first ‘subdomain’, then the
domain decomposition is based on

$\Omega=\Omega\,\cup\,\\{\Omega_{B}^{1}\,\cup\,\cdots\,\cup\,\Omega_{B}^{D}\\}=\Omega\,\cup\Omega_{B}.$
(19)

Figure 5 illustrates the domain decomposition view of the global–local
smoother.

Figure 5: Domain decomposition view of the combined global–local smoother,
which consists of the global smoother $S_{g}$ on the whole domain and the
local residual correction $S_{c}$ on the locally poor quality subdomains.

This combined smoother can be viewed as a full overlapped Schwarz-type domain
decomposition method [34, 9]. It is also noted in [9] that for the Schwarz-
type domain decomposition method the larger the subdomain overlap the faster
the convergence. From the algebraic point of view, the combined smoother is
also in the form of a block Gauss–Seidel. The combined smoother is similar to
a plane smoother [37, 24, 25], which is a variant of block Gauss–Seidel. The
plane smoother applies a standard smoother on multi-block structured grids,
e.g. on $x$\- or $y$-planes and is known for its effectiveness for strongly
anisotropic problems [29]. It is also noted that in [11, 12] a multigrid
preconditioner is proposed that solves systems on local patches for the
augmented momentum block of a finite element discretisation of the
incompressible Navier–Stokes equations.

###  4.4 Identifying low quality regions

There is no universal or binary measure of cell quality, so we heuristically
select a cell quality measure and set a threshold value for what constitutes a
low quality cell. For the examples in section 6, which use simplex cells, we
define low quality regions $\Omega_{B}$ by:

$\displaystyle\Omega_{b}$ $\displaystyle=\\{K:\ \gamma(K)<0.1\\},$ (20)
$\displaystyle\Omega_{B}$ $\displaystyle=\Omega_{b}\cup\\{K:K\text{ shares
vertex with }\Omega_{b}\\},$

where the region $\Omega_{b}$ contains all cells $K$ with a normalised radius
ratio $\gamma$ less than a specified threshold (see eq. 6 for the definition
of the radius ratio). In the case of Lagrange elements, the entire local
correction region $\Omega_{B}$ includes the region $\Omega_{b}$ and its one
layer (by vertex) extension, which involves all cells sharing degrees-of-
freedom with $\Omega_{b}$.

Figure 6 illustrates the local correction region for the mesh in fig. 1, with
low quality cells being those with normalised radius ratio $\gamma<0.1$ in
this case. The low quality cells in $\Omega_{b}$ are coloured blue. The cells
in the one layer extension are coloured orange. The local correction region
$\Omega_{B}$ is the entire coloured region.

Figure 6: Local correction region $\Omega_{B}$ in the low quality unit square
mesh (fig. 1). The area coloured blue is $\Omega_{b}$, the one-layer extension
of $\Omega_{b}$ is coloured orange. The local correction region $\Omega_{B}$
is the union of coloured cells.

##  5 Chebyshev smoothers: adjusted eigenvalues

The Chebyshev semi-iterative method [16] is commonly used as a multigrid
smoother [1], and requires knowledge of the largest and smallest eigenvalues
of the operator. For effective application as a smoother, only the largest
eigenvalue needs to be approximated [1]. In practice, when used as smoother
the estimated smallest eigenvalue is set to be specified fraction of the
estimated largest eigenvalue.

In the presence of localised low quality cells, the maximum eigenvalue
increases relative to a comparable resolution high quality mesh. For a mesh
with a small region of low quality cells, experiments indicate that the
increased maximum eigenvalue may render the smoother less effective in the
bulk of the domain. Consider a block-wise decomposition of the discretised
system

$A=\begin{bmatrix}A_{gg}&A_{gb}\\\ A_{bg}&A_{bb}\end{bmatrix},$ (21)

where $A_{bb}$ corresponds to the DOFs associated with the closure of
$\Omega_{B}$. The _adjusted eigenvalue_ is taken as the largest eigenvalue of
$A_{gg}$, i.e., $\lambda_{\max}\leftarrow\lambda_{\max}(A_{gg})$. We explore
numerically in section 6 differences in performance using estimates of the
largest eigenvalue of $A$ compared to the adjusted largest eigenvalue.

##  6 Numerical examples

We examine numerically the performance of non-nested geometric Galerkin
multigrid with the proposed smoother for model Poisson and elasticity
problems. Problems are solved using Lagrange elements on simplices using
linear ($P_{1}$) and quadratic ($P_{2}$) bases. The solver is terminated once
the relative residual reaches $10^{-10}$, measured in the $\ell^{2}$-norm. For
Chebyshev smoothers, the Krylov–Schur method [36, 18] is used to compute the
largest eigenvalue, using a tolerance of $10^{-8}$.

All meshes (grids) are unstructured and the levels are non-nested. We restrict
ourselves to problems in which the geometry can be exactly represented by the
coarsest grid. To study the influence of cell quality, we displace some
vertices in the generated meshes to change, controllably, mesh quality
characteristics. Low quality regions $\Omega_{B}$ are identified as regions
with normalised cell radius ratios of less than $0.1$ and their one layer
extension (given in eq. 20). We use $\gamma_{\min}$ to denote the minimum
normalised radius ratio in a mesh (grid). When distinguishing between high-
and low-quality meshes we use the annotation $\Omega_{\text{high}}$ and
$\Omega_{\text{low}}$, respectively.

We consider two cases:

  * Case A: All grids (meshes) have low quality regions.

  * Case B: All grids below the finest grid have low quality regions. The finest grid is high quality.

Reference performance is taken as the performance for problems with high
quality meshes for all levels.

The examples are computed using libraries from the FEniCS Project [2, 26, 27]
and PETSc [5, 4, 3]. The Galerkin finite element systems and the
restriction/prolongation operators are constructed with FEniCS, and the
remaining generic multigrid functionality is provided by PETSc. The largest
eigenvalue, used in the Chebyshev smoother, is estimated using the SLEPc
library [17, 32]. Meshes for each level are generated using Gmsh [14]. The
source code is freely available in the supporting material [7].

###  6.1 Poisson problem

For a domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ with boundary
$\Gamma:=\partial\Omega$ that is partitioned such that
$\Gamma=\Gamma_{D}\cup\Gamma_{N}$ and $\Gamma_{D}\cap\Gamma_{N}=\emptyset$,
the Poisson problem reads:

$\displaystyle-\nabla^{2}u$ $\displaystyle=f\quad\text{in}\ \Omega,$ (22)
$\displaystyle u$ $\displaystyle=g\quad\text{on}\ \Gamma_{D},$
$\displaystyle\nabla u\cdot n$ $\displaystyle=s\quad\text{on}\ \Gamma_{N}.$

where $f$, $g$ and $s$ are prescribed. We examine performance for two- and
three-dimensional problems. Multigrid is used as a solver for the Poisson
examples (as opposed to as a preconditioner).

####  6.1.1 Unit square domain

We test the homogeneous Poisson problem in section 3.1 on $\Omega=(0,1)^{2}$
with $f=0$, $\Gamma_{D}=\Gamma$ and $g=0$ using linear elements ($P_{1}$) and
a two-level V-cycle. An initial guess $u^{(0)}=\sin(10\pi x)\sin(10\pi y)$ is
interpolated onto the fine grid (see fig. 3(a)). The fine grid has $158$
vertices with $272$ cells and is shown in fig. 1 after perturbation of some
vertices. The coarse grid has $68$ cells. One iteration of the smoother
(symmetric Gauss–Seidel as the global smoother) is applied in the pre- and
post-smoothing steps. The norm of the relative residual after each multigrid
cycle is shown in fig. 7 for (i) a low quality mesh with the standard
smoother, (ii) a low quality mesh with the local correction smoother and (iii)
a high quality mesh (reference case).

Figure 7: Relative residual with and without the local correction for the
Poisson problem on a unit square.

Slow convergence in the residual with the low quality mesh is clear, whereas
the local correction smoother recovers the rate of the reference case. Figure
8 shows the absolute value of the solution error on the low quality fine grid
(Case A) after five multigrid V-cycles, with and without the local correction.

(a) without local correction

(b) with local correction

Figure 8: Absolute value of the error after five multigrid cycles on the low
quality unit square fine grid (a) without local correction and (b) with local
correction.

It is clear that the error persists in a localised region around the low
quality cells without the local correction, and this error is removed by the
local correction smoother.

We next consider a four-level V-cycle, with symmetric Gauss–Seidel as the
standard global smoother. The smoother is applied twice in the pre- and post-
smoothing steps. We perturb the position of some vertices on each level to
generate the low quality meshes (levels). The low quality finest grid is shown
in fig. 9, with the regions of low cell quality indicated.

Figure 9: A unit square mesh with three regions of low cell quality.

The number of cells for each level and the number of DOFs for $P_{1}$ and
$P_{2}$ elements are summarised in table 1, as well as the cell quality
measures and sizes of the degraded regions.

level | $\gamma_{\min}(\Omega_{\text{high}})$ | $\gamma_{\min}(\Omega_{\text{low}})$ |  | number of
---
cells in $\Omega$
| number of
---
cells in $\Omega_{B}$
| number of
---
DOFs in $\Omega$
$P_{1}/P_{2}$
| number of
---
DOFs in $\Omega_{B}$
$P_{1}/P_{2}$
$1$ (fine) | $0.661$ | $1.53\times 10^{-4}$ | $4236$ | $49$ |  $2199$/$8633$ |  $44$/$134$
$2$ | $0.673$ | $3.10\times 10^{-4}$ | $1016$ | $41$ |  $549$/$2133$ |  $37$/$112$
$3$ | $0.772$ | $2.09\times 10^{-4}$ | $254$ | $36$ |  $148$/$549$ |  $34$/$100$
$4$ (coarse) | $0.773$ | $9.08\times 10^{-4}$ | $68$ | $11$ |  $45$/$157$ |  $10$/$30$
Table 1: Cell quality of minimum normalized radius ratio $\gamma$ and the
problem size on each level of the four–level unit square hierarchy meshes.

We solve the Poisson problem on $\Omega=(0,1)^{2}$ with
$\Gamma_{D}=\\{(x,y)\in\partial\Omega:y=0,1\\}$, $f=2\pi^{2}\cos(\pi
x)\sin(\pi y)$, $g=0$ and $s=0$. Figure 10 presents the computed relative
residual for Case A and Case B, with and without local correction, for $P_{1}$
and $P_{2}$ elements. It is clear that the convergence rate is slow for low
quality meshes, and particularly so for quadratic elements. In all cases the
local correction smoother recovers the convergence rate of the reference case
with high quality meshes.

(a) $P_{1}$ elements

(b) $P_{2}$ elements

Figure 10: Relative residual after each multigrid cycle (symmetric
Gauss–Seidel as global smoother) for the Poisson problem on a unit square
domain with low quality regions on all levels (Case A) and with low quality
regions on all levels except the finest level (Case B).

We consider the error in the solution vector for Case A, Case B and the
reference case after 12 multigrid cycles, with the error given by $\lVert
u_{\rm MG}-u_{\rm LU}\rVert_{\ell^{2}}$, where $u_{\rm MG}$ is the multigrid
solution vector and $u_{\rm LU}$ is the solution vector computed using LU
factorisation. The computed discrete error is presented in table 2.

solver error | $\lVert u_{\rm MG}-u_{\rm LU}\rVert_{\ell^{2}}$
---|---
element type | $P_{1}$ | $P_{2}$
reference case | $2.86\times 10^{-8}$ | $7.04\times 10^{-9}$
case | Case A | Case B | Case A | Case B
without local correction | $2.09\times 10^{-6}$ | $1.57\times 10^{-7}$ | $3.91\times 10^{-5}$ | $3.76\times 10^{-5}$
with local correction | $8.04\times 10^{-8}$ | $6.35\times 10^{-8}$ | $1.10\times 10^{-8}$ | $1.01\times 10^{-8}$
Table 2: Solution vector error in the $\ell^{2}$ norm after $12$ multigrid
cycles with and without local correction for the Poisson problem on a unit
square domain.

With the proposed smoother, the solver error is reduced, particularly for
$P_{2}$ elements. The finite element solution error in the $\mathcal{L}^{2}$
norm for $P_{2}$ elements is shown in table 3 after different numbers of
cycles.

FEM solution error $\lVert u_{h}-u\rVert_{\mathcal{L}^{2}(\Omega)}$
---
| number of
---
cycles ($k$)
reference case |  | Case A
---
without local correction
| Case A
---
with local correction
2 | $1.28\times 10^{-3}$ | $2.18\times 10^{-3}$ | $1.46\times 10^{-3}$
5 | $6.05\times 10^{-6}$ | $3.67\times 10^{-5}$ | $1.04\times 10^{-5}$
10 | $2.56\times 10^{-6}$ | $3.83\times 10^{-6}$ | $2.70\times 10^{-6}$
| number of
---
cycles ($k$)
reference case |  | Case B
---
without local correction
| Case B
---
with local correction
2 | $1.28\times 10^{-3}$ | $2.11\times 10^{-3}$ | $1.37\times 10^{-3}$
5 | $6.05\times 10^{-6}$ | $3.64\times 10^{-5}$ | $9.74\times 10^{-6}$
10 | $2.56\times 10^{-6}$ | $3.61\times 10^{-6}$ | $2.56\times 10^{-6}$
Table 3: Finite element solution error in the $\mathcal{L}^{2}$-norm for
solving the Poisson equation on a unit square domain with and without local
correction for $P_{2}$ elements, where $u_{h}$ is the computed solution (after
a specified number of multigrid cycles) and $u$ is the exact solution to the
Poisson problem.

The error is greatest for low quality meshes without the local correction.
With the local correction, the finite element solution error is reduced. For
Case B, which has a high quality fine grid, the accuracy with the local
correction is the same as the reference case after ten cycles. Low quality
intermediate grids have not degraded accuracy and the local correction
overcomes slow solver convergence. Figure 11 shows the absolute value of
discrete residual on the finest grid after 10 V-cycles for the case where low
quality cells appear on all levels and local correction is only applied on
coarse grids. Unsurprisingly, the large value residual persists in the regions
with low quality cells.

(a) $P_{1}$ element

(b) $P_{2}$ element

Figure 11: Absolute value of residual at each vertex of the low quality finest
unit square grid after ten multigrid cycles, without local correction, for the
four-level Poisson problem using $P_{1}$ and $P_{2}$ elements.

####  6.1.2 Unit cube domain

We consider $\Omega=(0,1)^{3}$ with
$\Gamma_{D}=\\{(x,y,z)\in\partial\Omega:x=0,1\\}$,
$f=10\exp[-((x-0.5)^{2}+(y-0.5)^{2}+(z-0.5)^{2})/0.02]$, $g=0$ and $s=0$. We
use a four-level multigrid V-cycle for this problem. The coarsening rate in
terms of number of degrees-of-freedom is in the range of 7–8 at each level.
The cell quality measures and cell sizes for all grid levels
($\Omega_{\text{high}}$ and $\Omega_{\text{low}}$) are summarised in table 4.
A histogram of the normalised radius ratio for the Case A (low quality finest
grid) is shown in fig. 12.

level | $\gamma_{\min}(\Omega_{\text{high}})$ | $\gamma_{\min}(\Omega_{\text{low}})$ |  | number of
---
cells in $\Omega$
| number of
---
cells in $\Omega_{B}$
| number of
---
DOFs in $\Omega$
$P_{1}/P_{2}$
| number of
---
DOFs in $\Omega_{B}$
$P_{1}/P_{2}$
1 (fine) | $0.275$ | $7.20\times 10^{-6}$ | $582\,730$ | $625$ |  $104\,976$/$814\,775$ |  $231$/$1262$
2 | $0.288$ | $2.77\times 10^{-6}$ | $65\,259$ | $604$ |  $13\,361$/$97\,422$ |  $244$/$1294$
3 | $0.278$ | $1.74\times 10^{-6}$ | $7165$ | $568$ |  $1776$/$11\,845$ |  $233$/$1226$
4 (coarse) | $0.336$ | $4.64\times 10^{-6}$ | $792$ | $116$ |  $251$/$1501$ |  $50$/$258$
Table 4: Cell quality of minimum normalized radius ratio $\gamma$ and the
problem size on each level of the unit cube hierarchy of meshes.

(a) Full range of normalised radius ratio.

(b) Number of cells with a normalised radius ratio of less than $0.1$.

Figure 12: Histogram of the cell normalised radius ratio $\gamma$ (eq. 6) for
a unit cube fine grid with low quality regions.

Figure 13 presents the relative residual after each multigrid cycle using
symmetric Gauss–Seidel smoother as the global smoother, with and without the
local correction. Two iterations of the smoother are applied in pre- and in
post-smoothing. The reduction of the residual is again slow with the standard
smoother for cases with a low quality grid level, particularly for $P_{2}$
elements. The local correction restores the convergence rate to that of the
high quality mesh reference case.

(a) $P_{1}$ element

(b) $P_{2}$ element

Figure 13: Relative residual after each multigrid cycle (symmetric
Gauss–Seidel as global smoother) for for the Poisson problem on the unit cube
domain with low quality regions on all levels (Case A) and with low quality
regions on all levels except the finest level (Case B).

We now replace the symmetric Gauss–Seidel smoother with a Jacobi
preconditioned Chebyshev smoother and test performance for Case A (all levels
contain low quality cells). The smallest eigenvalue is approximated as one
tenth of the largest eigenvalue. We use the (i) largest eigenvalue
(unadjusted) and the (ii) adjusted largest eigenvalue in the Chebyshev
smoother. The largest eigenvalues for high and low quality meshes, as well as
the adjusted largest eigenvalues (section 5), are presented in table 5.

largest eigenvalue
---
element type | $P_{1}$ | $P_{2}$
level | level 1 | level 2 | level 3 | level 1 | level 2 | level 3
high quality mesh | 2.0787 | 2.0227 | 2.3897 | 2.4544 | 2.5602 | 2.2711
low quality mesh | 3.7103 | 2.6002 | 3.8458 | 4.4079 | 9.5812 | 10.1333
| adjusted largest eigenvalue
---
for low quality mesh
2.0788 | 2.1126 | 2.3756 | 2.4543 | 3.4003 | 2.8149
Table 5: Largest eigenvalues on high and low quality unit cube meshes, and the
adjusted largest eigenvalues for low quality unit cube meshes (Poisson
problem).

The residual after each cycle is presented in fig. 14. The convergence rate
improves using the adjusted eigenvalues in the Chebyshev smoother, and matches
the reference case.

(a) $P_{1}$ element

(b) $P_{2}$ element

Figure 14: Relative residual after each multigrid cycle (Chebyshev global
smoother, using maximum eigenvalues (unadjusted) and using adjusted
eigenvalues) for Poisson problem on the unit cube mesh with low quality
regions on all levels (Case A).

###  6.2 Linear elasticity

We now consider the linearised elasticity problem,

$-\nabla\cdot\sigma(u)=f\quad\text{in}\quad\Omega,$ (23)

where $u$ is the displacement field and $\sigma(u)$ is the stress tensor
satisfying the isotropic elastic law

$\sigma(u):=2\mu\epsilon(u)+\lambda\text{tr}(\epsilon(u))I,$ (24)

where $\epsilon(u)$ is the strain,

$\epsilon(u):=\frac{1}{2}\left(\nabla u+(\nabla u)^{T}\right),$ (25)

and $\mu:=E/2(1+\nu)$ and $\lambda:=E\nu/(1+\nu)(1-2\nu)$ are the Lamé
parameters, $E$ is the Young’s modulus and $\nu$ is Poisson’s ratio. For each
example we use $E=6.9\times 10^{10}$ and $\nu=0.33$.

For elasticity examples, multigrid is used as a preconditioner for the
conjugate gradient (CG) method, with one multigrid V-cycle applied at each CG
iteration. In each multigrid cycle, the smoother is applied twice in pre- and
twice in post-smoothing. All examples use four multigrid levels.

####  6.2.1 Lattice domain

The domain tested here is a truss-like lattice structure
$\Omega=a^{3}\setminus(a\times b^{2}\cup b\times a\times b\cup b^{2}\times
a)$, with $a=[0,6]$ and $b=[1,5]$, which is illustrated in fig. 15. The
locations of low quality regions for the low quality fine grid are indicated
in fig. 15.

(a) Low quality cells

(b) unstructured mesh

Figure 15: Lattice geometry and the positions of low quality cells on the
finest grid.

We consider $f=(0,0,0)$ and boundary conditions

$\displaystyle u=(0,0,0)$
$\displaystyle\text{on}\quad\Gamma_{1}=\\{(x,y,z)\in\partial\Omega:x=0\\},$
(26) $\displaystyle\sigma\cdot n=(10^{3},0,0)$
$\displaystyle\text{on}\quad\Gamma_{2}=\\{(x,y,z)\in\partial\Omega:x=6\\},$
$\displaystyle\sigma\cdot n=0$
$\displaystyle\text{on}\quad\partial\Omega\setminus\\{\Gamma_{1}\cup\Gamma_{2}\\}.$

Data on the size and cell quality of each grid level for both high and low
quality meshes is summarised in table 6.

level | $\gamma_{\min}(\Omega_{\text{high}})$ | $\gamma_{\min}(\Omega_{\text{low}})$ |  | number of
---
cells in $\Omega$
| number of
---
cells in $\Omega_{B}$
| number of
---
DOFs in $\Omega$
$P_{1}/P_{2}$
| number of
---
DOFs in $\Omega_{B}$
$P_{1}/P_{2}$
1 (fine) | $0.209$ | $2.64\times 10^{-8}$ | $711.683$ | $732$ |  $437\,946$/$3\,192\,258$ |  $825$/$4506$
2 | $0.228$ | $2.70\times 10^{-8}$ | $78\,100$ | $585$ |  $58\,992$/$391\,767$ |  $798$/$4101$
3 | $0.250$ | $6.43\times 10^{-8}$ | $8341$ | $444$ |  $8292$/$49\,131$ |  $612$/$3027$
4 (coarse) | $0.344$ | $1.07\times 10^{-6}$ | $1080$ | $96$ |  $1284$/$7116$ |  $183$/$816$
Table 6: Cell quality of minimum normalized radius ratio $\gamma$ and the
problem size on each level of the lattice hierarchy meshes.

With symmetric Gauss–Seidel applied as the global smoother, the computed
residual after each CG iteration, for $P_{1}$ and $P_{2}$ elements, is shown
in fig. 16 for Case A and Case B.

(a) $P_{1}$ element

(b) $P_{2}$ element

Figure 16: Relative residual after each multigrid (symmetric Gauss–Seidel as
global smoother) preconditioned CG iteration with and without local correction
for the linear elasticity problem on the lattice mesh with low quality regions
on all levels (Case A) and low quality regions on all levels except the finest
level (Case B).

The convergence rate with low quality meshes at all levels (Case A) is poor
compared to the reference case. The performance for Case B (high quality fine
grid, low quality coarse grids) appears reasonable for $P_{1}$ elements, but
is considerably worse than the reference case for $P_{2}$ elements. Using the
local correction smoother, the reference rate of convergence is recovered.
Figure 17 shows the locations on the finest grid where, after 10 iterations,
the residual is large for Case A (absolute value of the residual is greater
than $10^{-6}$). Comparing to fig. 15, the large residuals coincide with
regions of low cell quality.

(a) $P_{1}$ element

(b) $P_{2}$ element

Figure 17: Locations on the finest lattice grid where the residual is large
after ten iterations of multigrid preconditioned CG.

####  6.2.2 Dumbbell-like structure

Finally we consider a dumbbell-like structure with hexagonal ends connected by
three slender bars, shown in fig. 18. We take $f=(0,0,-10^{3})$ and boundary
conditions: $u=(0,0,0)$ on the left-hand most boundary, $\sigma\cdot
n=(10^{3},0,0)$ on the right-hand most boundary, and $\sigma\cdot n=(0,0,0)$
on other parts of boundary.

(a) level 1 for $P_{1}$ element

(b) level 1 for $P_{2}$ element

Figure 18: A dumbbell-like structure and the positions of low quality cells on
the finest grid of the $P_{1}$ and $P_{2}$ elements.

Summary data for the meshes is given in table 7. Different meshes are used for
the $P_{1}$ and $P_{2}$ element cases. Note that the coarsening rates vary
between levels, with more aggressive coarsening from the simulation (finest)
level to the second level.

$P_{1}$ element
---
level | cell size | $\gamma_{\min}(\Omega_{\text{high}})$ | $\gamma_{\min}(\Omega_{\text{low}})$ |  | number of
---
cells in $\Omega$
| number of
---
cells in $\Omega_{B}$
| number of
---
DOFs in $\Omega$
| number of
---
DOFs in $\Omega_{B}$
1 (fine) | $0.025$ | $0.150$ | $3.18\times 10^{-8}$ | $3\,982\,354$ | $788$ | $2\,130\,582$ | $876$
2 | $0.08$ | $0.204$ | $6.94\times 10^{-8}$ | $139\,469$ | $538$ | $90\,027$ | $699$
3 | $0.2$ | $0.273$ | $2.93\times 10^{-8}$ | $12\,627$ | $454$ | $10\,455$ | $615$
4 (coarse) | $0.5$ | $0.179$ | $8.92\times 10^{-8}$ | $1345$ | $164$ | $1515$ | $216$
$P_{2}$ element
level | cell size | $\gamma_{\min}(\Omega_{\text{high}})$ | $\gamma_{\min}(\Omega_{\text{low}})$ |  | number of
---
cells in $\Omega$
| number of
---
cells in $\Omega_{B}$
| number of
---
DOFs in $\Omega$
| number of
---
DOFs in $\Omega_{B}$
1 (fine) | $0.04$ | $0.220$ | $3.26\times 10^{-8}$ | $1\,035\,068$ | $600$ | $4\,453\,713$ | $3771$
2 | $0.12$ | $0.211$ | $6.27\times 10^{-8}$ | $48\,880$ | $381$ | $237\,942$ | $2664$
3 | $0.25$ | $0.267$ | $8.34\times 10^{-8}$ | $5059$ | $253$ | $28\,863$ | $1926$
4 (coarse) | $0.5$ | $0.179$ | $8.92\times 10^{-8}$ | $1345$ | $164$ | $8532$ | $1116$
Table 7: Cell quality data (minimum normalized radius ratio $\gamma$) and
problem size for each level of the dumbbell–like hierarchy meshes.

We use a Jacobi–preconditioned Chebyshev global smoother for this example. The
smallest eigenvalue in the Chebyshev smoother is set as one tenth of the
largest eigenvalue. The computed residuals at each CG iteration with the
adjusted largest eigenvalue and with and without the local correction are
shown in fig. 19. Convergence is slow with low quality meshes and a standard
smoother, and particularly so when the fine grid contains low quality cells
(Case A). As for the lattice example, for Case B the $P_{2}$ convergence rate
is more heavily impacted by the low quality cells than the $P_{1}$ case.
Applying the local correction restores the observed convergence close to the
reference case.

(a) $P_{1}$ element

(b) $P_{2}$ element

Figure 19: Relative residual after each Krylov iteration for the multigrid
(Chebyshev global smoother with adjusted largest eigenvalue) preconditioned CG
solver for solving linear elasticity on the dumbbell-like structure. Case A
(low quality regions on all levels) and Case B (low quality regions on levels
except the finest level) are presented.

The number of CG iterations required to reduce the relative residual to below
$10^{-10}$ are listed in table 8, including with and without the adjusted
largest eigenvalues in the Chebyshev smoother. The convergence rate is
improved by using the adjusted maximum eigenvalue compared to using the
largest eigenvalue of the whole system.

element type |  $P_{1}$ element |  $P_{2}$ element
---|---|---
reference case | $17$ |  | $24$ |
case | Case A | Case B | Case A | Case B
low quality mesh without correction | $>100$ | $25$ | $\gg 100$ | $49$
| low quality mesh with correction
---
using unadjusted $\lambda_{\max}$ in Chebyshev
$25$ | $22$ | $48$ | $40$
| low quality mesh with correction
---
using adjusted $\lambda_{\max}$ in Chebyshev
$18$ | $17$ | $25$ | $24$
Table 8: Number of multigrid (Chebyshev global smoother) preconditioned CG
iterations needed to reduce the relative residual to $10^{-10}$ for linear
elasticity on the the dumbbell–like structure with low quality regions on all
levels (Case A) and low quality regions on levels except the finest level
(Case B).

##  7 Conclusions

Geometric multigrid on non-nested, unstructured grids has been considered in
the presence of a small number of low quality cells, which can characterise
meshes of geometrically complex domains (at the very least in intermediate
grids). It was observed that the performance of the geometric multigrid method
degrades significantly when a mesh contains a small number of low quality
cells, with the poor convergence attributed to the local failure of smoothers
in regions close to the low quality cells. A global–local combined smoother
was developed to overcome this issue. The smoother involves application of a
global smoother on the entire grid combined with a local correction on
subdomains with low cell quality.

We have demonstrated the proposed smoother on several numerical examples. It
eliminates errors in low cell quality regions that are not removed by a
standard smoother. It was shown that the slow convergence rate for low quality
meshes can be restored to the high quality mesh reference rate by applying the
combined smoother. In particular, if the fine grid is high quality and only
coarse grids have low quality regions, the finite element discretisation error
is barely influenced and the local correction improves the convergence rate of
multigrid. This work improves the robustness for multigrid on complex
geometric domains and opens up the possibility of the high performance,
geometric, scalable multigrid solvers to solve complicated engineering
applications at a system level.

The non-nested geometric approach can offer an alternative to algebraic
multigrid in some cases. A question is the extension of the approach to
algebraic multigrid (AMG). Work is currently underway investigating coarsening
strategies, transfer operators and smoothers for AMG that are robust with
respect to cell quality.

## Acknowledgements

YC was supported by the Youth Program of the Natural Science Foundation of
Jiangsu Province (No. BK20230466), the Jiangsu Funding Program for Excellent
Postdoctoral Talent (No. 2022ZB584), and Jiangsu Shuangchuang Project
(JSSCTD202209). GNW acknowledges support from the Engineering and Physical
Sciences Research Council under grants EP/V001396/1, EP/S005072/1,
EP/W00755X/1 and EP/W026635/1.

## References

  * Adams et al. [2003] Mark Adams, Marian Brezina, Jonathan Hu, and Ray Tuminaro.  Parallel multigrid smoothing: polynomial versus Gauss–Seidel.  _Journal of Computational Physics_ , 188(2):593–610, 2003.  doi: 10.1016/S0021-9991(03)00194-3.
  * Alnæs et al. [2015] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells.  The FEniCS Project Version 1.5.  _Archive of Numerical Software_ , 3(100):9–23, 2015.
  * Balay et al. [1997] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith.  Efficient management of parallelism in object oriented numerical software libraries.  In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, _Modern Software Tools in Scientific Computing_ , pages 163–202. Birkhäuser Press, 1997.
  * Balay et al. [2023a] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc/TAO users manual.  Technical Report ANL-21/39 - Revision 3.20, Argonne National Laboratory, 2023a.
  * Balay et al. [2023b] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil M. Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc Web page, 2023b.  URL https://petsc.org/.
  * Briggs et al. [2000] William L. Briggs, Van E. Henson, and Steve F. McCormick.  _A Multigrid Tutorial_.  SIAM, second edition, 2000.
  * Chen [2023] Yuxuan Chen.  Supporting material.  URL https://github.com/yc397/multigrid_examples, 2023.
  * Dickopf [2010] Thomas Dickopf.  _On multilevel methods based on non-nested meshes_.  PhD thesis, University of Bonn, 2010.
  * Dolean et al. [2015] Victorita Dolean, Pierre Jolivet, and Frédéric Nataf.  _An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation_.  SIAM, 2015.
  * Duprez et al. [2019] Michel Duprez, Vanessa Lleras, and Alexei Lozinski.  Finite element method with local damage on the mesh.  _ESAIM: M2AN_ , 53(6):1871–1891, 2019.  doi: 10.1051/m2an/2019023.
  * Farrell et al. [2019] Patrick E. Farrell, Lawrence Mitchell, and Florian Wechsung.  An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier–Stokes equations at high Reynolds number.  _SIAM Journal on Scientific Computing_ , 41(5):A3073–A3096, 2019.  doi: 10.1137/18M1219370.
  * Farrell et al. [2021] Patrick E. Farrell, Lawrence Mitchell, L. Ridgway Scott, and Florian Wechsung.  A Reynolds-robust preconditioner for the Scott-Vogelius discretization of the stationary incompressible Navier-Stokes equations.  _SMAI Journal of Computational Mathematics_ , 7:75–96, 2021.  doi: 10.5802/smai-jcm.72.
  * Freitag and Ollivier-Gooch [2000] Lori A. Freitag and Carl Ollivier-Gooch.  A cost/benefit analysis of simplicial mesh improvement techniques as measured by solution efficiency.  _International Journal of Computational Geometry & Applications_, 10(4):361–382, 2000.  doi: 10.1142/S0218195900000218.
  * Geuzaine and Remacle [2009] Christophe Geuzaine and Jean-François Remacle.  Gmsh: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities.  _International Journal for Numerical Methods in Engineering_ , 79(11):1309–1331, 2009.  doi: 10.1002/nme.2579.
  * Hackbusch [2013] Wolfgang Hackbusch.  _Multi-Grid Methods and Applications_ , volume 4.  Springer Science & Business Media, 2013.
  * Hageman and Young [2012] Louis A. Hageman and David M. Young.  _Applied Iterative Methods_.  Courier Corporation, 2012.
  * Hernandez et al. [2005] Vicente Hernandez, Jose E. Roman, and Vicente Vidal.  SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems.  _ACM Transactions on Mathematical Software_ , 31(3):351–362, 2005.  doi: 10.1145/1089014.1089019.
  * Hernández et al. [2007] Vicente Hernández, Jose E Román, Andrés Tomás, and Vicente Vidal.  Krylov-Schur methods in SLEPc.  _Universitat Politecnica de Valencia, Tech. Rep. STR-7_ , 2007.
  * Katz and Sankaran [2011] Aaron Katz and Venkateswaran Sankaran.  Mesh quality effects on the accuracy of CFD solutions on unstructured meshes.  _Journal of Computational Physics_ , 230(20):7670–7686, 2011.  doi: 10.1016/j.jcp.2011.06.023.
  * Klingner and Shewchuk [2008] Bryan Matthew Klingner and Jonathan Richard Shewchuk.  Aggressive tetrahedral mesh improvement.  In _Proceedings of the 16th International Meshing Roundtable_ , pages 3–23. Springer, Berlin, Heidelberg, 2008.
  * Knupp [2000a] Patrick M. Knupp.  Achieving finite element mesh quality via optimization of the Jacobian matrix norm and associated quantities. Part I – A framework for surface mesh optimization.  _International Journal for Numerical Methods in Engineering_ , 48(3):401–420, 2000a.  doi: 10.1002/(SICI)1097-0207(20000530)48:3¡401::AID-NME880¿3.0.CO;2-D.
  * Knupp [2000b] Patrick M. Knupp.  Achieving finite element mesh quality via optimization of the Jacobian matrix norm and associated quantities. Part II – A framework for volume mesh optimization and the condition number of the Jacobian matrix.  _International Journal for Numerical Methods in Engineering_ , 48(8):1165–1185, 2000b.  doi: 10.1002/(SICI)1097-0207(20000720)48:8¡1165::AID-NME940¿3.0.CO;2-Y.
  * Knupp [2001] Patrick M. Knupp.  Algebraic mesh quality metrics.  _SIAM Journal on Scientific Computing_ , 23(1):193–218, 2001.  doi: 10.1137/S1064827500371499.
  * Llorente and Melson [2000] Ignacio Martín Llorente and Norman Duane Melson.  Behavior of plane relaxation methods as multigrid smoothers.  _Electronic Transactions on Numerical Analysis_ , 10:92–114, 2000.
  * Llorente et al. [2000] Ignacio Martín Llorente, Boris Diskin, and Norman Duane Melson.  Alternating plane smoothers for multiblock grids.  _SIAM Journal on Scientific Computing_ , 22(1):218–242, 2000.  doi: 10.1137/S106482759935736X.
  * Logg and Wells [2010] Anders Logg and Garth N. Wells.  DOLFIN: Automated finite element computing.  _ACM Transactions on Mathematical Software_ , 37(2):20, 2010.  doi: 10.1145/1731022.1731030.
  * Logg et al. [2012] Anders Logg, Kent-Andre Mardal, and Garth N. Wells.  _Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book_ , volume 84.  Springer Science & Business Media, 2012.
  * Mavriplis [1992] D. J. Mavriplis.  Three-dimensional unstructured multigrid for the Euler equations.  _AIAA Journal_ , 30(7):1753–1761, 1992.  doi: 10.2514/3.11133.
  * Oosterlee [1997] Cornelis Willebrordus Oosterlee.  A GMRES–based plane smoother in multigrid to solve 3D anisotropic fluid flow problems.  _Journal of Computational Physics_ , 130(1):41–53, 1997.  doi: 10.1006/jcph.1996.5442.
  * Parthasarathy et al. [1994] V. N. Parthasarathy, C. M. Graichen, and A. F. Hathaway.  A comparison of tetrahedron quality measures.  _Finite Elements in Analysis and Design_ , 15(3):255–261, 1994.  doi: 10.1016/0168-874X(94)90033-7.
  * Richardson et al. [2019] Chris N. Richardson, Nathan Sime, and Garth N. Wells.  Scalable computation of thermomechanical turbomachinery problems.  _Finite Elements in Analysis and Design_ , 155:32–42, 2019.  doi: 10.1016/j.finel.2018.11.002.
  * Roman et al. [2023] J. E. Roman, C. Campos, L. Dalcin, E. Romero, and A. Tomas.  SLEPc users manual.  Technical Report DSIC-II/24/02 - Revision 3.20, D. Sistemes Informàtics i Computació, Universitat Politècnica de València, 2023.
  * Shewchuk [2002] Jonathan Shewchuk.  What is a good linear finite element? interpolation, conditioning, anisotropy, and quality measures.  In _Proceedings of the 11th International Meshing Roundtable_ , pages 115–126. IMR, 2002.
  * Smith et al. [2004] Barry Smith, Petter Bjorstad, and William D. Gropp.  _Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations_.  Cambridge University Press, 2004.
  * Smith et al. [1993] K. M. Smith, W. K. Cope, and Surya Pratap Vanka.  A multigrid procedure for three-dimensional flows on non-orthogonal collocated grids.  _International Journal for Numerical Methods in Fluids_ , 17(10):887–904, 1993.  doi: 10.1002/fld.1650171005.
  * Stewart [2002] G. W. Stewart.  A Krylov–Schur algorithm for large eigenproblems.  _SIAM Journal on Matrix Analysis and Applications_ , 23(3):601–614, 2002.  doi: 10.1137/S0895479800371529.
  * Trottenberg et al. [2000] Ulrich Trottenberg, Cornelius W. Oosterlee, and Anton Schuller.  _Multigrid_.  Academic Press, London, 2000.
  * Xu [1992] Jinchao Xu.  Iterative methods by space decomposition and subspace correction.  _SIAM Review_ , 34(4):581–613, 1992.  doi: 10.1137/1034116.
  * Xu [2001] Jinchao Xu.  The method of subspace corrections.  _Journal of Computational and Applied Mathematics_ , 128(1):335–362, 2001.  doi: 10.1016/S0377-0427(00)00518-5.