*Article* **Spectral Analysis of the Finite Element Matrices Approximating 3D Linearly Elastic Structures and Multigrid Proposals**

**Quoc Khanh Nguyen 1,\*, Stefano Serra-Capizzano 2,3, Cristina Tablino-Possio <sup>4</sup> and Eddie Wadbro 1,5**

	- <sup>5</sup> Department of Mathematics and Computer Science, Karlstad University, SE-651 88 Karlstad, Sweden
	- **\*** Correspondence: qukh0001@student.umu.se

**Abstract:** The so-called material distribution methods for topology optimization cast the governing equation as an extended or fictitious domain problem, in which a coefficient field represents the design. In practice, the finite element method is typically used to approximate that kind of governing equations by using a large number of elements to discretize the design domain, and an elementwise constant function approximates the coefficient field in that domain. This paper presents a spectral analysis of the coefficient matrices associated with the linear systems stemming from the finite element discretization of a linearly elastic problem for an arbitrary coefficient field in three spatial dimensions. The given theoretical analysis is used for designing and studying an optimal multigrid method in the sense that the (arithmetic) cost for solving the problem, up to a fixed desired accuracy, is linear in the corresponding matrix size. Few selected numerical examples are presented and discussed in connection with the theoretical findings.

**Keywords:** matrix sequences; spectral analysis; finite element approximations

#### **1. Introduction**

In our previous paper [1], we applied the theory of generalized locally Toeplitz (GLT) sequences to compute and analyze the asymptotic spectral distribution of the sequence of stiffness matrices {*Kn*}*n*, with *K<sup>n</sup>* being the finite element (FE) approximation of the considered one spatial dimension topology optimization problem, for a given fineness parameter associated to *n*. In a later contribution [2], we extended the analysis to the two-dimensional setting using so-called multilevel block GLT sequences. In this paper, we further expand the theory to cover the three-dimensional case.

Since the first material distribution method for topology design was introduced in the late 1980s [3], topology optimization [4,5], a well-known computational tool for finding the optimal distribution of material within a given design domain, has been studied extensively. The material distribution topology optimization has contributed to the development of several areas, such as electromagnetic [6–8], fluid–structure interaction [9,10], acoustics [11,12], additive manufacturing [13], and especially (non-)linear elasticity [14–16]. For problems in linear elasticity, which motivates the study in this paper, the most common method to solve this type of problem is the so-called density-based or material distribution approach. In this approach, a so-called material indicator function *α*(*x*)—typically referred to as the density or the physical design—models the presence/absence of material; *α* = 1 where material is present, else *α* = 0. However, the binary design problem is computationally intractable. A standard approach to make the problem computationally feasible is employing a combination of relaxation, penalization, and filtering, in which the physical density is defined as *ρ*(*x*) = *α* + (1 − *α*)*g* F(*α*)(*x*) by *ρ* ≥ 0, which is a constant, the penalization operator

**Citation:** Nguyen, K.N.; Serra-Capizzano, S.; Tablino-Possio, C.; Wadbro, E. Spectral Analysis of the Finite Element Matrices Approximating 3D Linearly Elastic Structures and Multigrid Proposals. *Math. Comput. Appl.* **2022**, *27*, 78. https://doi.org/ 10.3390/mca27050078

Academic Editors: Sascha Eisenträger and Gianluigi Rozza

Received: 13 April 2022 Accepted: 10 September 2022 Published: 14 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

is *g*, and the filtering procedure is F. In this approach, the finite element method (FEM) is the standard choice for generating numerical solutions of a linearly elastic boundary value problem. In practice, the physical design is typically represented as an element-wise constant function. Moreover, a combination of allowing intermediate values of the design, filtering, and penalization has been crucial for the success of these methods. In this paper, we limit our attention to studying the linear system that stems from the FE discretization of the governing equation and its solution. From a high performance computing perspective, this limitation is natural since the computational effort in solving topology optimization problems is dominated by the cost of solving the discretized governing equations—a fact that is particularly prominent when studying three-dimensional problems. Over the last decades, there has been significant work in improving the efficiency of topology optimization, aiming at solving large-scale design problems [17–21]. Although our study is directly motivated by the material-distribution or density-based approach to topology optimization, the work is relevant for other design descriptions, such as level sets [22,23] or moving morphable components [24–27], provided that the physics description uses the so-called ersatz material approach, which in some cases can be justified rigorously [4].

The theory of multilevel block GLT sequences [28,29] is a generalization of the GLT sequences theory [30,31] that is typically used for computing/analyzing the spectral distribution of matrix sequences arising from, for example, the numerical discretization, such as the FE approximation of partial differential equations (PDEs) with proper boundary conditions. In the considered matrix sequences, the size of the given linear systems *dn* increases with *n*, in which *dn* tends to infinity as *n* → ∞. In this paper, the entire sequence of linear systems with increasing size arising from the three spatial dimension problems is the primary consideration. Under mild conditions, essentially relying on regularity assumptions on the meshes, we show that the sequence of discretization matrices has an asymptotic spectral distribution. By leveraging the theory of multilevel block GLT sequences, we now extend our previous works [1,2] further to perform a detailed spectral analysis of the linear systems associated with the FE discretization of the governing equation in the three-dimensional setting. Similar to our previous work [2], we also make use of the information obtained from the spectral symbol *f* to design a fast, multigrid solver in three dimensions for optimizing the (arithmetic) cost to solve the related linear systems up to a fixed desired accuracy, which is proportional to the matrix–vector cost, which is linear in the corresponding matrix size. This solver is also verified and comes up with very satisfying numerical results, in terms of the linear cost and number of iterations, which are bounded independently of the matrix size and mildly depending on the desired accuracy.

In the following sections, we will go into detail on the problem description, spectral analysis, multigrid method, and a brief conclusion. The description of the continuous problem and the resulting coefficient matrices arising from our FE approximation is delivered in Sections 2 and 3 is devoted to the spectral analysis of the FE matrices from the perspective of the GLT theory. In Section 4, a brief account of multigrid methods with special attention to the block case encountered in the present context is given, and the spectral information is a core component in the development of the multigrid proposal for our specific 3D setting presented in Section 5. Eventually, the conclusions are reported in Section 6, and some relevant model information is explained in Appendice A and B.

#### **2. Problem Description**

We consider a linearly elastic structure that occupies (part of) the hyper-rectangular domain <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup>3. In particular, we are interested in the setting used in material distribution based topology optimization, where a function, typically denoted as the physical density, *ρ* : Ω → [0, 1], describes the layout of the unloaded structure. We assume that the structure is clamped along the boundary portion <sup>Γ</sup>*<sup>D</sup>* <sup>⊂</sup> *<sup>∂</sup>*Ω. Moreover, we let *<sup>b</sup>* <sup>∈</sup> *<sup>L</sup>*2(Ω)<sup>3</sup> be a given body load (a volume force) in <sup>Ω</sup>, *<sup>t</sup>* <sup>∈</sup> *<sup>L</sup>*2(Γ*F*)<sup>3</sup> be the surface traction acting on the

non-clamped boundary Γ*<sup>F</sup>* ⊂ *∂*Ω of the solid, and *u* denote the resulting equilibrium displacement, which solves the following problem.

$$\text{Find } u \in \mathcal{U} \text{ such that } a(u, v; \rho) = \ell(v) \,\,\forall v \in \mathcal{U}, \tag{1}$$

where the set of all kinematically admissible displacements of the structure is

$$\mathcal{U} = \left\{ \mu \in H^1(\Omega)^3 \mid \mu|\mathbf{r}\_D \equiv 0 \right\}.$$

The energy bilinear form *a* and the load linear form are defined as

$$a(\mu, \upsilon; \rho) = \int\_{\Omega} \rho \left( E\_{\mathsf{c}} \epsilon(\mathsf{u}) \right) : \epsilon(\upsilon),$$

$$\ell(\upsilon) = \int\_{\Gamma\_F} \mathsf{t} \cdot \upsilon + \int\_{\Omega} \mathsf{b} \cdot \upsilon$$

where (*u*) = <sup>∇</sup>*<sup>u</sup>* <sup>+</sup> <sup>∇</sup>*u<sup>T</sup>* /2 is the strain tensor of *u*. The colon ":" denotes the full contraction between the two tensors; when using the standard basis, the full contraction of the two matrices is their Frobenius scalar product. *Ec* is a constant fourth-order elasticity tensor. In this paper, we study an FE discretization of problem (1), in which the physical density is approximated by an element-wise constant function, which is typical in material distribution based topology optimization. The domain Ω is discretized into *n* trilinear hexahedral elements and then applying FE approximation, the variational problem (1) is reduced to the linear system

$$\mathbb{K}\_n(\rho)u = f\_\prime$$

where *u* and *f* are the nodal displacement and load vector, respectively, and *Kn*(*ρ*) is the stiffness matrix of element-wise constant physical density function *ρ*—the entries of the vector *ρ* = [*ρ*1, *ρ*2, ... , *ρn*] are the element values of *ρ*; that is, *ρ<sup>i</sup>* is the value of *ρ* in the *i*th element in the FE mesh. The stiffness matrix *K<sup>n</sup>* is typically assembled by looping over each element so that

$$\mathbb{K}\_n(\rho) = \sum\_{i=1}^n \rho\_i \mathbb{K}\_\varepsilon^{(i)}\prime$$

where *K*(*i*) *<sup>e</sup>* is the element stiffness matrix. We emphasize that the formal expression of the relevant matrices is a key ingredient for applying the multilevel block GLT theory to produce a global spectral description of the matrix sequences under consideration. More precisely, the non-zero blocks of three-dimensional element stiffness matrix can be expressed as

$$\mathbf{K}\_{\varepsilon} = \frac{h}{2} \frac{E\_0}{1+\nu} \begin{bmatrix} \mathbf{K}\_1 & \mathbf{K}\_2 & \mathbf{K}\_3 & \mathbf{K}\_4 & \mathbf{K}\_5 & \mathbf{K}\_6 & \mathbf{K}\_7 & \mathbf{K}\_8\\ & \mathbf{K}\_9 & \mathbf{K}\_{10} & \mathbf{K}\_{11} & \mathbf{K}\_{12} & \mathbf{K}\_{13} & \mathbf{K}\_{14} & \mathbf{K}\_7\\ & & \mathbf{K}\_{15} & \mathbf{K}\_{16} & \mathbf{K}\_{17} & \mathbf{K}\_{18} & \mathbf{K}\_{19} & \mathbf{K}\_6\\ & & & \mathbf{K}\_{19} & \mathbf{K}\_{20} & \mathbf{K}\_{17} & \mathbf{K}\_{12} & \mathbf{K}\_5\\ & & & & \mathbf{K}\_{19} & \mathbf{K}\_{16} & \mathbf{K}\_{11} & \mathbf{K}\_4\\ & & & & & \mathbf{K}\_{15} & \mathbf{K}\_{10} & \mathbf{K}\_3\\ & & & & & & & \mathbf{K}\_9 & \mathbf{K}\_2\\ & & & & & & & \mathbf{K}\_1 \end{bmatrix} \tag{2}$$

where each *Ki* block has size 3 × 3 and is associated with one node on a so-called reference element, *E*<sup>0</sup> is Young's modulus, and *ν* is Poisson's ratio. More details are provided in Appendix A for the explicit expressions for the element stiffness matrix and in Appendix B for the stress–strain relation and various bounds. With regard to Equation (2), the lower part of the matrix is not explicitly given because the matrix is globally real symmetric and so is every block *Kl*, *l* = 1, ... , 20, whose expressions are reported in Appendix A, Equation (A1).

#### **3. Spectral Analysis**

The current section is devoted to the spectral analysis of the FE coefficient matrices derived in the previous section and is complemented by a selection of numerical tests that confirm the theoretical analysis (more numerical experiments have been performed but here we report only a selection for the sake of brevity). We limit our focus to isotropic materials; for more information about various bounds, see the discussion in Appendix B. From a mathematics perspective, this means that we limit our attention to when the Poisson's ratio *ν* is in the range [0, 0.5). In particular, Section 3.1 contains the minimal set of preliminary concepts and tools, while Sections 3.2 and 3.3 are focused on the specific study in 3D in the constant and variable coefficient cases, respectively.

#### *3.1. Premises*

The premises include the formal definition of the eigenvalue and singular value distribution, the notion of multi-indexing, the concepts of multilevel block Toeplitz matrices, multilevel block sampling matrices, and multilevel block GLT matrix sequences.

#### 3.1.1. Singular Value/Eigenvalue Distributions

We first give the formal definitions, and then we briefly discuss their informal and practical meaning.

**Definition 1.** *Let r*, *t be two positive integers. Let* {*An*}*<sup>n</sup> be a sequence of matrices, with An of size dn with eigenvalues λ*1(*An*), ... , *λdn* (*An*) *and singular values σ*1(*An*), ... , *σdn* (*An*)*. Furthermore, let <sup>f</sup>* : *<sup>D</sup>* <sup>⊂</sup> <sup>R</sup>*<sup>t</sup>* <sup>→</sup> <sup>C</sup>*r*×*<sup>r</sup> be a measurable function defined on a set <sup>D</sup> with* <sup>0</sup> <sup>&</sup>lt; *<sup>μ</sup>t*(*D*) <sup>&</sup>lt; <sup>∞</sup>*, and with <sup>μ</sup>t*(·) *denoting the Lebesgue measure on* <sup>C</sup>*<sup>t</sup> .*

• *We say that* {*An*}*<sup>n</sup> has a (asymptotic) singular value distribution described by f , and we write* {*An*}*<sup>n</sup>* ∼*<sup>σ</sup> f , if*

$$\lim\_{n \to \infty} \frac{1}{d\_n} \sum\_{i=1}^{d\_n} F(\sigma\_i(A\_{\mathbb{R}})) = \frac{1}{\mu\_t(D)} \int\_D \frac{\sum\_{i=1}^r F(\sigma\_i(f(\mathbf{x})))}{r} d\mathbf{x}, \quad \forall F \in \mathbb{C}\_{\mathbb{C}}(\mathbb{R}), \tag{3}$$

*with σ*1(*f*), ... , *σr*(*f*) *being the singular values of f , each of them being a measurable function non-negative almost everywhere.*

• *We say that* {*An*}*<sup>n</sup> has a (asymptotic) spectral (or eigenvalue) distribution described by f , and we write* {*An*}*<sup>n</sup>* ∼*<sup>λ</sup> f , if*

$$\lim\_{n \to \infty} \frac{1}{d\_n} \sum\_{i=1}^{d\_n} F(\lambda\_i(A\_n)) = \frac{1}{\mu\_t(D)} \int\_D \frac{\sum\_{i=1}^r F(\lambda\_i(f(\mathbf{x})))}{r} d\mathbf{x}, \quad \forall F \in \mathbb{C}\_t(\mathbb{C}), \tag{4}$$

*with λ*1(*f*), ... , *λr*(*f*) *being the eigenvalues of f , each of them being a complex-valued measurable function.*

*If* {*An*}*<sup>n</sup> has both a singular value and an eigenvalue distribution described by f , we write* {*An*}*<sup>n</sup>* ∼*σ*,*<sup>λ</sup> f . (In practice, in the Toeplitz setting and often in the GLT setting, the parameter r can be read at a matrix level as the size of the elementary blocks which form the global matrix An, as it will be clear both in our stiffness matrices in Sections 3.2 and 3.3, and in the block Toeplitz/block diagonal sampling structures in Section 3.1.2.)*

The symbol *f* contains spectral/singular value information briefly described informally as follows. With reference to relation (4), the eigenvalues of *Kn* are partitioned into *r* subsets of the same cardinality, except possibly for a small number of outliers, such that the *i*th subset is approximately formed by the samples of *λi*(*f*) over a uniform grid in *D*, *i* = 1, ... ,*r*. Thus, provided that *n* is large enough, the symbol *f* provides a 'compact' and

quite accurate description of the spectrum of the matrices *Kn*. Similarly, relation (3) has the same meaning when talking of the singular values of *K<sup>n</sup>* and by replacing *λi*(*f*) with *σi*(*f*), *i* = 1, . . . ,*r*.

3.1.2. Multilevel Block Toeplitz Matrices, Multilevel Block Diagonal Sampling Matrices, and Multilevel Block GLT Sequences

The present section is specifically devoted to matrix-theoretic notations and definitions, which are essential when dealing with multi-dimensional problems and with the related sequences of matrices.

**Definition 2.** *A multi-index <sup>i</sup>* <sup>∈</sup> <sup>Z</sup>*d, also called a d-index, is a (row) vector in* <sup>Z</sup>*d, whose components are denoted by i*1,..., *id.*


$$\left[\begin{array}{c} \dots \dots \, \left[\begin{array}{c} \left[\begin{array}{c} (j\_1, \dots, \dots, j\_d) \end{array} \right]\_{j\_d = h\_d, \dots, k\_d} \end{array} \right]\_{j\_{d-1} = h\_{d-1}, \dots, k\_{d-1}} \dots \right]\_{j\_1 = h\_1, \dots, k\_1} \dots \right]$$

With regard to the previous definition, in the case *d* = 2, the lexicographic ordering is the following: (*h*1, *h*2), (*h*1, *h*<sup>2</sup> + 1), ... , (*h*1, *k*2), (*h*<sup>1</sup> + 1, *h*2),(*h*<sup>1</sup> + 1, *h*<sup>2</sup> + 1), ...(*h*<sup>1</sup> + 1, *k*2), ...(*k*1, *h*1), (*k*1, *h*<sup>1</sup> + 1), ... , (*k*1, *k*2). Notice that, in general, a multi-index range [*h*, ... , *k*], *h* ≤ *k* is used with *h* = **0** or with *h* = **1**.

**Definition 3** ((Multilevel) Block Toeplitz Matrices)**.** *Given* **<sup>n</sup>** <sup>∈</sup> <sup>N</sup>*d, a matrix of the form*

$$[A\_{\mathbf{i}-\mathbf{j}}]\_{\mathbf{i},\mathbf{j}=\mathbf{e}}^{\mathbf{n}} \in \mathbb{C}^{N(\mathfrak{n})r \times N(\mathfrak{n})r}$$

*with* **<sup>e</sup>** *vector of all ones, with blocks <sup>A</sup>***<sup>k</sup>** <sup>∈</sup> <sup>C</sup>*r*×*<sup>r</sup> ,* **k** = −(**n** − **e**), ... , **n** − **e***, is called a multilevel block Toeplitz matrix, or, more precisely, a d-level r-block Toeplitz matrix. Let φ* : [−*π*, *π*] *<sup>d</sup>* <sup>→</sup> <sup>C</sup>*r*×*<sup>r</sup> be a matrix-valued function in which each entry belongs to <sup>L</sup>*1([−*π*, *<sup>π</sup>*] *<sup>d</sup>*)*. We denote the Fourier coefficients of the generating function φ as*

$$\hat{f}\_{\mathbf{k}} = \frac{1}{(2\pi)^{d}} \int\_{[-\pi,\pi]^{d}} \phi(\theta) e^{-\mathbb{f}(\mathbf{k},\theta)} d\theta \in \mathbb{C}^{r \times r}, \quad \mathbf{k} \in \mathbb{Z}^{d}.$$

*where the integrals are computed componentwise, ı*ˆ <sup>2</sup> <sup>=</sup> <sup>−</sup>1*, and* (**k**, *<sup>θ</sup>*) = *<sup>k</sup>*1*θ*<sup>1</sup> <sup>+</sup> ... <sup>+</sup> *kdθd. For every* **<sup>n</sup>** <sup>∈</sup> <sup>N</sup>*d, the* **<sup>n</sup>***th Toeplitz matrix associated with <sup>φ</sup> is defined as*

$$T\_{\mathfrak{n}}(\phi) := [\hat{f}\_{\mathfrak{i}-\mathfrak{j}}]\_{\mathfrak{i},\mathfrak{j}-\mathfrak{e}}^{\mathfrak{n}}$$

*or, equivalently, as*

$$T\_{\mathbf{n}}(\boldsymbol{\phi}) = \sum\_{|j\_1| \le m\_1} \dots \sum\_{|j\_d| \le m\_d} [J\_{n\_1}^{(j\_1)} \otimes \dots \otimes J\_{n\_d}^{(j\_d)}] \otimes \hat{f}\_{(j\_1, \dots, j\_d) \times \dots}$$

*where* ⊗ *denotes the (Kronecker) tensor product of matrices, while J* (*l*) *<sup>m</sup> is the matrix of order m whose* (*i*, *<sup>j</sup>*) *entry equals* <sup>1</sup> *if <sup>i</sup>* <sup>−</sup> *<sup>j</sup>* <sup>=</sup> *<sup>l</sup> and zero otherwise. We call* {*T***n**(*φ*)}**n**∈N*<sup>d</sup> the family of (multilevel block) Toeplitz matrices associated with φ, which, in turn, is called the generating function of* {*T***n**(*φ*)}**n**∈N*<sup>d</sup> .*

**Definition 4** ((Multilevel) Block Diagonal Sampling Matrices)**.** *For n* ∈ N *and a* : [0, 1] → C*r*×*<sup>r</sup> , we define the block diagonal sampling matrix Dn*(*a*) *as the diagonal matrix*

$$D\_n(a) = \mathop{\text{diag}}\_{i=1,\dots,n} a\left(\frac{i}{n}\right) = \begin{bmatrix} a(\frac{1}{n}) & & & \\ & a(\frac{2}{n}) & & \\ & & \ddots & \\ & & & a(1) \end{bmatrix} \in \mathbb{C}^{rn \times rn} \text{ .} $$

*For <sup>n</sup>* <sup>∈</sup> <sup>N</sup>*<sup>d</sup> and <sup>a</sup>* : [0, 1] *<sup>d</sup>* <sup>→</sup> <sup>C</sup>*r*×*<sup>r</sup> , we define the multilevel block diagonal sampling matrix Dn*(*a*) *as the block diagonal matrix*

$$D\_{\mathfrak{u}}(a) = \underset{i=1,\ldots,n}{\text{diag}} \; a\left(\frac{\mathfrak{i}}{\mathfrak{u}}\right) \in \mathbb{C}^{rN(\mathfrak{u}) \times rN(\mathfrak{u})} \_{\text{'} }$$

*with the lexicographical ordering discussed at the beginning of Section 3.1.2.*

**Definition 5** (Zero-Distributed Sequences of Matrices)**.** *According to Definition 1, a sequence of matrices* {*Zn*}*<sup>n</sup> such that*

$$\{Z\_{\mathfrak{m}}\}\_{\mathfrak{m}} \sim\_{\sigma} 0$$

*is referred to as a zero-distributed sequence. Note that, for any r* ≥ 1*,* {*Zn*}*<sup>n</sup>* ∼*<sup>σ</sup>* 0 *is equivalent to* {*Zn*}*<sup>n</sup>* ∼*<sup>σ</sup> Or (notice that Om and Im denote the m* × *m zero matrix and the m* × *m identity matrix, respectively).*

Proposition 1 provides an important characterization of zero-distributed sequences together with a useful sufficient condition for detecting such sequences. Throughout this paper, we use the natural convention 1/∞ = 0.

**Proposition 1.** [30] *Let* {*Zn*}*<sup>n</sup> be a sequence of matrices, with Zn of size dn, and let* · *be the standard spectral matrix norm (the one induced by the Euclidean vector norm).*


**(Multilevel) Block GLT Matrix Sequences.** Now, we give a very concise and operational description of the multilevel block GLT sequences, from which it will be clear that the multilevel block Toeplitz structures, the zero-distributed matrix sequences, and the multilevel block diagonal sampling matrices represent the basic building components. All the material is taken from the books [30,31] and from the papers [28,32].

Let *d*,*r* ≥ 1 be fixed positive integers. A multilevel *r*-block GLT sequence (or simply a GLT sequence if *d*,*r* can be inferred from the context or we do not need/want to specify them) is a special *r*-block matrix sequence {*An*}*<sup>n</sup>* equipped with a measurable function *κ* : [0, 1] *<sup>d</sup>* <sup>×</sup> [−*π*, *<sup>π</sup>*] *<sup>d</sup>* <sup>→</sup> <sup>C</sup>*r*×*<sup>r</sup>* , the so-called symbol. We use the notation {*An*}*<sup>n</sup>* ∼GLT *κ* to indicate that {*An*}*<sup>n</sup>* is a GLT sequence with symbol *κ*. The symbol of a GLT sequence is unique in the sense that if {*An*}*<sup>n</sup>* ∼GLT *κ* and {*An*}*<sup>n</sup>* ∼GLT *ς* then *κ* = *ς* a.e. in [0, 1] *d* × [−*π*, *π*] *<sup>d</sup>*. The main properties of *r*-block GLT sequences proved in [28] are listed below. If *A* is a matrix, we denote by *A*† the Moore–Penrose pseudoinverse of *A* (recall that *A*† = *A*−<sup>1</sup> whenever *A* is invertible).

**GLT 1** If {*An*}*<sup>n</sup>* ∼GLT *κ* then {*An*}*<sup>n</sup>* ∼*<sup>σ</sup> κ*. If moreover each *An* is Hermitian, then {*An*}*<sup>n</sup>* ∼*<sup>λ</sup> κ*.

**GLT 2** We have the following:


• {*Zn*}*<sup>n</sup>* ∼GLT *κ*(**x**, *θ*) = *Or* if and only if {*Zn*}*<sup>n</sup>* ∼*<sup>σ</sup>* 0 (zero-distributed sequences coincide exactly with the GLT sequences having GLT symbol equal to *Or* a.e.).

**GLT 3** If {*An*}*<sup>n</sup>* ∼GLT *κ* and {*Bn*}*<sup>n</sup>* ∼GLT *ς*, then


### *3.2. Constant Coefficient Case ρ* ≡ 1

In the current section, we report the derivation and the formal expression of the symbol for the matrix sequences in the constant coefficient setting *ρ* ≡ 1, according to the standard notion of generating function in the Toeplitz theory (see Section 3.1.2 and the paper by Garoni et al. [29] for more details) and according to the notion of symbol reported in Definition 1.

#### 3.2.1. Symbol Definition

Recall that the computational domain Ω is a hyper-rectangular domain in 3D, and the boundary of Ω comprises six sides. Let *An*(1, DN5) be the stiffness matrix obtained with a Q<sup>1</sup> FE approximation with proper boundary conditions, Dirichlet 'D' on one of Ω's sides and Neumann 'N' in the other five sides of Ω (and hence the formal notation *An*(1, DN5)), where we have chosen a uniform meshing with *n* intervals in the *x*<sup>1</sup> direction, *n* intervals in the *x*<sup>2</sup> direction, and *n* intervals in the *x*<sup>3</sup> direction. According to the previously considered ordering of the nodes, the matrix *An*(1, DN5) is a three-level block tridiagonal structure of size *n* with two-level tridiagonal blocks of size *n* + 1, with tridiagonal blocks of size *n* + 1, whose elements are small matrices of size 3. We notice that the size is dictated by all the meshpoints, including those in the boundaries when considering Neumann boundary conditions, while only the internal meshpoints are involved when the Dirichlet boundary conditions are enforced.

In accordance with the 2D case [2], if all the boundary conditions are of the Dirichlet type, then we obtain a matrix *An*(1, D6), which is a three-level block Toeplitz structure with elementary blocks of size 3, having global dimension 3(*<sup>n</sup>* <sup>−</sup> <sup>1</sup>)3, since all the nodes on the boundaries are not unknowns. More precisely, we obtain

$$A\_{\boldsymbol{\eta}}(\mathbf{1}, \mathbf{D}^{\boldsymbol{\delta}}) = \begin{bmatrix} a\_0 & a\_{-1} \\ a\_1 & a\_0 & a\_{-1} \\ & \ddots & \ddots & \ddots \\ & & a\_1 & a\_0 & a\_{-1} \\ & & & a\_1 & a\_0 \end{bmatrix} = \text{triddiag} \ (a\_{\boldsymbol{\delta}})\_{\boldsymbol{i} = -1, 0, 1}$$

and, for *i* = −1, 0, 1, we have

$$a\_i = \begin{bmatrix} a\_{i0} & a\_{i-1} \\ a\_{i1} & a\_{i0} & a\_{i-1} \\ & \ddots & \ddots & \ddots \\ & & a\_{i1} & a\_{i0} & a\_{i-1} \\ & & & a\_{i1} & a\_{i0} \end{bmatrix} = \text{tridig}(a\_{ij})\_{j=-1,0,1,1}$$

while, for *i*, *j* = −1, 0, 1, the following block tridiagonal matrices are defined:

$$a\_{ij} = \begin{bmatrix} a\_{ij0} & a\_{ij-1} \\ a\_{ij1} & a\_{ij0} & a\_{ij-1} \\ & \ddots & \ddots & \ddots \\ & & a\_{ij1} & a\_{ij0} & a\_{ij-1} \\ & & & a\_{ij1} & a\_{ij0} \end{bmatrix} = \text{tridig}(a\_{ijk})\_{k=-1,0,1}.$$

By reading the entries of *An*(1, D6), analogously to the process in the 2D setting as in [2], we can compute explicitly the related generating function, following the rules given in Definition 3:

*<sup>f</sup>*Q<sup>1</sup> (*θ*1, *<sup>θ</sup>*2, *<sup>θ</sup>*3) = *fa*<sup>0</sup> (*θ*1, *<sup>θ</sup>*2) + *fa*−<sup>1</sup> (*θ*1, *<sup>θ</sup>*2)*<sup>e</sup>* <sup>−</sup>*ı*ˆ*θ*<sup>3</sup> + *fa*<sup>1</sup> (*θ*1, *<sup>θ</sup>*2)*<sup>e</sup> ı*ˆ*θ*<sup>3</sup> = - *fa*<sup>00</sup> (*θ*1) + *fa*<sup>0</sup>−<sup>1</sup> (*θ*1)*<sup>e</sup>* <sup>−</sup>*ı*ˆ*θ*<sup>2</sup> + *fa*<sup>01</sup> (*θ*1)*<sup>e</sup> ı*ˆ*θ*<sup>2</sup> + - *fa*−<sup>10</sup> (*θ*1) + *fa*−1−<sup>1</sup> (*θ*1)*<sup>e</sup>* <sup>−</sup>*ı*ˆ*θ*<sup>2</sup> <sup>+</sup> *fa*−<sup>11</sup> (*θ*1)*<sup>e</sup> ı*ˆ*θ*<sup>2</sup> *e* −*ı*ˆ*θ*<sup>3</sup> + - *fa*<sup>10</sup> (*θ*1) + *fa*<sup>1</sup>−<sup>1</sup> (*θ*1)*<sup>e</sup>* <sup>−</sup>*ı*ˆ*θ*<sup>2</sup> + *fa*<sup>11</sup> (*θ*1)*<sup>e</sup> ı*ˆ*θ*<sup>2</sup> *e ı*ˆ*θ*<sup>3</sup> = - *a*<sup>000</sup> + *a*00−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> + *a*001*e ı*ˆ*θ*1 + - *a*0−<sup>10</sup> + *a*0−1−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*0−11*<sup>e</sup> ı*ˆ*θ*1 *e* −*ı*ˆ*θ*<sup>2</sup> + - *a*<sup>010</sup> + *a*01−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> + *a*011*e ı*ˆ*θ*1 *e ı*ˆ*θ*<sup>2</sup> + - *a*−<sup>100</sup> + *a*−10−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*−101*<sup>e</sup> ı*ˆ*θ*1 + - *a*−1−<sup>10</sup> + *a*−1−1−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*−1−11*<sup>e</sup> ı*ˆ*θ*1 *e* −*ı*ˆ*θ*<sup>2</sup> + - *a*−<sup>110</sup> + *a*−11−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*−111*<sup>e</sup> ı*ˆ*θ*1 *e ı*ˆ*θ*<sup>2</sup> *e* −*ı*ˆ*θ*<sup>3</sup> + - *a*<sup>100</sup> + *a*10−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> + *a*101*e ı*ˆ*θ*1 + - *a*1−<sup>10</sup> + *a*1−1−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*1−11*<sup>e</sup> ı*ˆ*θ*1 *e* −*ı*ˆ*θ*<sup>2</sup> + - *a*<sup>110</sup> + *a*11−1*e* <sup>−</sup>*ı*ˆ*θ*<sup>1</sup> + *a*111*e ı*ˆ*θ*1 *e ı*ˆ*θ*<sup>2</sup> *e ı*ˆ*θ*3 ,

where every *arst*, *r*,*s*, *t* ∈ {−1, 0, 1} is a 3 × 3 matrix because three degrees of freedom are associated to each node of the mesh. Each *arst* is a sum of the building blocks of the elementary matrix (2). More precisely, in terms of the block matrices *Ki*, cf. definition (A1) in Appendix A, we can write

$$\begin{aligned} a\_{000} &= 2\left(K\_1 + K\_9 + K\_{15} + K\_{19}\right), & a\_{001} &= a\_{00-1} = 2\left(K\_2 + K\_{16}\right), \\ a\_{010} &= a\_{0-10} = 2\left(K\_3 + K\_{11}\right), & a\_{011} &= a\_{0-1-1} = 2K\_4, \\ a\_{01-1} &= a\_{0-11} = 2K\_{10}, & a\_{100} &= a\_{-100} = 2\left(K\_5 + 2K\_{13}\right), \\ a\_{101} &= a\_{-10-1} = 2K\_6, & a\_{10-1} &= a\_{-101} = 2K\_{12}, \\ a\_{110} &= a\_{-1-10} = 2K\_7, & a\_{1-10} &= a\_{-110} = 2K\_{17}, \\ a\_{111} &= a\_{-11-1} = K\_8, & a\_{11-1} &= a\_{-1-11} = K\_{14}, \\ a\_{1-11} &= a\_{-11-1} = K\_{18}, & a\_{-111} &= a\_{1-1-1} = K\_{20}. \end{aligned}$$

Thus, we have

$$f\_{\mathbb{Q}\_1}(\theta\_1, \theta\_2, \theta\_3) = \begin{bmatrix} f\_{11}(\theta\_1, \theta\_2, \theta\_3) & f\_{12}(\theta\_1, \theta\_2, \theta\_3) & f\_{13}(\theta\_1, \theta\_2, \theta\_3) \\ f\_{12}(\theta\_1, \theta\_2, \theta\_3) & f\_{22}(\theta\_1, \theta\_2, \theta\_3) & f\_{23}(\theta\_1, \theta\_2, \theta\_3) \\ f\_{13}(\theta\_1, \theta\_2, \theta\_3) & f\_{23}(\theta\_1, \theta\_2, \theta\_3) & f\_{33}(\theta\_1, \theta\_2, \theta\_3) \end{bmatrix},\tag{5}$$

where

*f*11(*θ*1, *θ*2, *θ*3) = 8*k*<sup>1</sup> + 8*k*<sup>3</sup> cos(*θ*1) + 2 cos(*θ*2) 4*k*<sup>3</sup> + 4*k*<sup>5</sup> cos(*x*) + 2 cos(*θ*3) - 4*k*<sup>6</sup> + 2 cos(*θ*2) 2*k*<sup>7</sup> + 2*k*<sup>9</sup> cos(*θ*1) + 4*k*<sup>7</sup> cos(*θ*1) , *f*12(*θ*1, *θ*2, *θ*3) = 4*k*<sup>2</sup> + 4*k*<sup>8</sup> + 4 cos(*θ*1) *k*<sup>4</sup> + *k*10) + 2*e* −*iθ*<sup>3</sup> - *e iθ*<sup>2</sup> *k*<sup>2</sup> + *k*<sup>4</sup> cos(*θ*1) + *e* −*iθ*<sup>2</sup> *k*<sup>8</sup> + *k*<sup>10</sup> cos(*θ*1) + 2*e iθ*<sup>3</sup> - *e iθ*<sup>2</sup> *k*<sup>8</sup> + *k*<sup>10</sup> cos(*θ*1) + *e* −*iθ*<sup>2</sup> *k*<sup>2</sup> + *k*<sup>4</sup> cos(*θ*1) , *f*13(*θ*1, *θ*2, *θ*3) = 4*k*<sup>2</sup> + 4*k*<sup>8</sup> + 8*k*<sup>0</sup> cos(*θ*1) + 4 cos(*θ*2) *k*<sup>4</sup> + *k*<sup>10</sup> + 2*k*<sup>0</sup> cos(*θ*1) + 2*e* −*iθ*<sup>3</sup> - *k*2*e <sup>i</sup>θ*<sup>1</sup> + *k*8*e* <sup>−</sup>*iθ*<sup>1</sup> + *k*4*e <sup>i</sup>θ*<sup>1</sup> + *k*10*e* −*iθ*<sup>1</sup> cos(*θ*2) + 2*e iθ*<sup>3</sup> - *k*2*e* <sup>−</sup>*iθ*<sup>1</sup> + *k*8*e <sup>i</sup>θ*<sup>1</sup> + *k*4*e* <sup>−</sup>*iθ*<sup>1</sup> + *k*10*e iθ*<sup>1</sup> cos(*θ*2) , *f*22(*θ*1, *θ*2, *θ*3) = 8*k*<sup>1</sup> + 8*k*<sup>3</sup> cos(*θ*1) + 8 cos(*θ*2) *k*<sup>6</sup> + *k*<sup>7</sup> cos(*θ*1) + 8 cos(*θ*3) - *k*<sup>3</sup> + *k*<sup>5</sup> cos(*θ*1) + cos(*θ*2) *k*<sup>7</sup> + *k*<sup>9</sup> cos(*θ*1) , *f*23(*θ*1, *θ*2, *θ*3) = 4*k*<sup>2</sup> + 4*k*<sup>8</sup> + 2*e* −*iθ*<sup>2</sup> *k*2*e <sup>i</sup>θ*<sup>1</sup> + *k*8*e* −*iθ*<sup>1</sup> + 2*e <sup>i</sup>θ*<sup>2</sup> (*k*2*e* <sup>−</sup>*iθ*<sup>1</sup> + *k*8*e <sup>i</sup>θ*<sup>1</sup> ) + 2 cos(*θ*3) - 2*k*<sup>4</sup> + 2*k*<sup>10</sup> + 4*k*<sup>0</sup> cos(*θ*1) + 4*k*<sup>0</sup> cos(*θ*2) + *e* −*iθ*<sup>2</sup> *k*4*e <sup>i</sup>θ*<sup>1</sup> + *k*10*e* −*iθ*<sup>1</sup> + *e iθ*<sup>2</sup> *k*4*e* <sup>−</sup>*iθ*<sup>1</sup> + *k*10*e iθ*<sup>1</sup> , *f*33(*θ*1, *θ*2, *θ*3) = 8*k*<sup>1</sup> + 8*k*<sup>6</sup> cos(*θ*1) + 8 cos(*θ*2) *k*<sup>3</sup> + *k*<sup>7</sup> cos(*θ*1) + 8 cos(*θ*3) - *k*<sup>3</sup> + *k*<sup>7</sup> cos(*θ*1) + cos(*θ*<sup>2</sup> (*k*<sup>5</sup> + *k*<sup>9</sup> cos(*θ*1) .

Finally, by expanding the *ki*s according to expressions (A2), we deduce the formal expression of the generating functions

$$\begin{split} f\_{11}(\theta\_{1},\theta\_{2},\theta\_{3}) &= 2\cos(\theta\_{3}) \left[ \cos(\theta\_{2}) \left( \frac{2\upsilon}{3(2\upsilon-1)} + 2\cos(\theta\_{1}) \left( \frac{\upsilon}{3(2\upsilon-1)} - \frac{1}{9} \right) - \frac{4}{9} \right) \right] \\ &+ \cos(\theta\_{1}) \left( \frac{2\upsilon}{3(2\upsilon-1)} - \frac{4}{9} \right) - \frac{8}{9} \\ &+ \frac{8}{9} \cos(\theta\_{1}) - \frac{16\upsilon}{3(2\upsilon-1)} \\ &+ \cos(\theta\_{2}) \left( 2\cos(\theta\_{1}) \left( \frac{2\upsilon}{3(2\upsilon-1)} + \frac{2}{9} \right) + \frac{8}{9} \right) + \frac{16}{9}, \\ f\_{22}(\theta\_{1},\theta\_{2},\theta\_{3}) &= 8\cos(\theta\_{3}) \left[ \cos(\theta\_{2}) \left( \frac{\upsilon}{6(2\upsilon-1)} + \cos(\theta\_{1}) \left( \frac{\upsilon}{6(2\upsilon-1)} - \frac{1}{18} \right) - \frac{1}{9} \right) \right] \\ &+ \cos(\theta\_{1}) \left( \frac{\upsilon}{6(2\upsilon-1)} + \frac{1}{18} \right) + \frac{1}{9} \Big] \\ &+ \frac{8}{9} \cos(\theta\_{1}) - \frac{16\upsilon}{3(2\upsilon-1)} \\ &+ 8 \cos(\theta\_{2}) \left( \cos(\theta\_{1}) \left( \frac{\upsilon}{6(2\upsilon-1)} - \frac{1}{9} \right) - \frac{2}{9} \right) + \frac{16}{9}, \end{split}$$

$$\begin{split} f\_{33}(\theta\_1, \theta\_2, \theta\_3) &= 8 \cos(\theta\_3) \left[ \cos(\theta\_2) \left( \frac{\upsilon}{6(2\upsilon - 1)} + \cos(\theta\_1) \left( \frac{\upsilon}{6(2\upsilon - 1)} - \frac{1}{18} \right) + \frac{1}{18} \right) \right] \\ &+ \cos(\theta\_1) \left( \frac{\upsilon}{6(2\upsilon - 1)} - \frac{1}{9} \right) + \frac{1}{9} \Big] \\ &- \frac{16\upsilon}{3(2\upsilon - 1)} - \frac{16}{9} \cos(\theta\_1) + \frac{16}{9} \\ &+ 8 \cos(\theta\_2) \left( \cos(\theta\_1) \left( \frac{\upsilon}{6(2\upsilon - 1)} - \frac{1}{9} \right) + \frac{1}{9} \right), \\ f\_{12}(\theta\_1, \theta\_2, \theta\_3) &= -\frac{4}{3(2\upsilon - 1)} \left( \upsilon \sin(\theta\_2) \sin(\theta\_3) \left( \cos(\theta\_1) + 2 \right) \right), \\ f\_{13}(\theta\_1, \theta\_2, \theta\_3) &= -\frac{4}{3(2\upsilon - 1)} \left( \upsilon \sin(\theta\_1) \sin(\theta\_3) \left( \cos(\theta\_2) + 2 \right) \right), \\ f\_{23}(\theta\_1, \theta\_2, \theta\_3) &= -\frac{4}{3(2\upsilon - 1)} \left( \upsilon \sin(\theta\_1) \sin(\theta\_2) \left( \cos(\theta\_3) + 2 \right) \right). \end{split}$$

The following proposition links in a precise way the involved Toeplitz structures and matrices *An*(1, D6) and *An*(1, DN5).

**Proposition 2.** *Let f*Q<sup>1</sup> (*θ*1, *θ*2, *θ*3) *be the symbol defined in the previous lines; see (5). Let us consider a uniform meshing in all the three directions with n subintervals. Then we have the following relationships:*

$$A\_{\mathfrak{n}}(1, \mathrm{D}^{6}) = T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_{1}}), \quad \mathfrak{n} = (n\_{1}, n\_{2}, n\_{3}), \; n\_{1} = n\_{2} = n\_{3} = n - 1,$$

$$A\_{\mathfrak{n}}(1, \mathrm{D}\mathrm{N}^{5}) = T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_{1}}) + R\_{\mathfrak{n}}, \quad \mathfrak{n} = (n\_{1}, n\_{2}), \; n\_{1} = n, n\_{2} = n\_{3} = n + 1,$$

*where R<sup>n</sup> has rank of order n*<sup>2</sup> = *o*(size(*Rn*))*, with* size(*Rn*) = size(*An*(1, DN5)) = *n*1*n*2*n*3*.*

*If the more general setting is considered with n*<sup>1</sup> *subintervals in the x*<sup>1</sup> *direction, n*<sup>2</sup> *subintervals in the x*<sup>2</sup> *direction, and n*<sup>3</sup> *intervals in the x*<sup>3</sup> *direction, then the analogous is true:*

$$A\_{\mathfrak{n}}(1, \operatorname{D}^{\mathfrak{f}}) = T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_{1}}), \quad \mathfrak{n} = (n\_{1} - 1, n\_{2} - 1, n\_{3} - 1),$$

$$A\_{\mathfrak{n}}(1, \operatorname{D}\operatorname{N}^{\mathfrak{f}}) = T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_{1}}) + R\_{\mathfrak{n}}, \quad \mathfrak{n} = (n\_{1}, n\_{2} + 1, n\_{3} + 1),$$

*where R<sup>n</sup> has a rank proportional to n*1*n*<sup>2</sup> + *n*1*n*<sup>3</sup> + *n*2*n*<sup>3</sup> = *o*(size(*Rn*))*, with* size(*Rn*) = size(*An*(1, DN5)) = *n*1*n*2*n*3*.*

**Proof.** Since the grid points in a cube are *n*1*n*2*n*3, if we consider Dirichlet boundary conditions in a facet orthogonal, for example, to the direction *x*1, then the number of points and equations which are affected by the Neumann boundary conditions on the other 5 facets are 2*l*3*n*1*n*<sup>2</sup> + 2*l*2*n*1*n*<sup>3</sup> + *l*1*n*2*n*3, where typically *l*<sup>3</sup> = *l*<sup>2</sup> = *l*<sup>1</sup> = 2 if the standard linear FEs are employed. Hence, the rank correction induced by the matrix *R<sup>n</sup>* is proportional, as claimed to *n*1*n*<sup>2</sup> + *n*1*n*<sup>3</sup> + *n*2*n*<sup>3</sup> = *o*(size(*Rn*)).

In the following section, we prove that the function *f*Q<sup>1</sup> is the eigenvalue symbol, in the sense of Definition 1, of the matrix sequences {*An*(1, <sup>D</sup>6)}*<sup>n</sup>* and {*An*(1, DN5)}*n*.

3.2.2. Symbol Spectral Analysis in 3D: Distribution, Extremal Eigenvalues, and Conditioning

We study the spectral distribution, extremal eigenvalues, and conditioning of our matrix sequences in 3D. All the results are based on the symbol and on its analytical features. **Proposition 3.** *Let f*Q<sup>1</sup> (*θ*1, *θ*2, *θ*3) *be the symbol defined in Section 3.2.1 According to the notation in Proposition 2, all the matrix sequences* {*Tn*(*f*Q<sup>1</sup> )}*n,* {*An*(1, D6)}*n,* {*An*(1, D6)}*n,* {*An*(1, DN5)}*n,* {*An*(1, DN5)}*<sup>n</sup> are spectrally distributed as f*Q<sup>1</sup> *in the sense of Definition 1.*

**Proof.** For the multilevel block Toeplitz sequences {*Tn*(*f*Q<sup>1</sup> )}*n*, {*An*(1, <sup>D</sup>6)}*n*, and {*An*(1, <sup>D</sup>6)}*<sup>n</sup>* refer to Item **GLT 2.**, part 1, and Item **GLT 1.**, part 2 in Section 3.1.2. For the sequences {*An*(1, DN5)}*n*, {*An*(1, DN5)}*<sup>n</sup>* first observe that the matrix sequence {*Rn*}*<sup>n</sup>* defined in Proposition 2 is zero distributed, thanks to Proposition 1 and Proposition 2 since rank(*Rn*)/size(*Rn*) → 0, as *n* → ∞. Then the claim follows thanks to the ∗-algebra structure of the GLT sequences and more specifically thanks to Item **GLT 3.**, part 2, Item **GLT 2.**, part 1, and Item **GLT 1.**, part 2.

Now, we identify a few key analytical features of the spectral symbol which are shared by all the matrix sequences mentioned in Proposition 3. Theorem 1 is crucial for the study of the extremal eigenvalues and of the conditioning of the same matrix sequences and, in Section 5, it is the main ingredient for designing ad hoc multigrid solvers when dealing with the associated large linear systems.

**Theorem 1.** *Let f*Q<sup>1</sup> (*θ*1, *θ*2, *θ*3) *be the symbol defined in (5). The following statements hold true:*


**Proof.** Claim 1. The function *f*Q<sup>1</sup> evaluated at (0, 0, 0) equals

$$f\_{\mathbb{Q}\_1}(0,0,0) = \begin{bmatrix} \alpha & \beta & \beta \\ \beta & \alpha & \beta \\ \beta & \beta & \alpha \end{bmatrix}$$

with *α* = 8(*k*<sup>1</sup> + 2*k*<sup>3</sup> + *k*<sup>5</sup> + *k*<sup>6</sup> + 2*k*<sup>7</sup> + *k*9) and *β* = 8(4*k*<sup>0</sup> + *k*<sup>2</sup> + *k*<sup>4</sup> + *k*<sup>8</sup> + *k*10), whose row-sum *α* + 2*β* = 0 according to (A2).

Claim 2. It is just a direct check. Indeed, it is enough to check that the quantities

$$\left(\text{tr}(f\_{\mathbb{Q}\_1}(\theta\_1, \theta\_2, \theta\_3))\right) \text{ and } \det(f\_{\mathbb{Q}\_1}(\theta\_1, \theta\_2, \theta\_3))$$

have a zero of orders two and six, respectively, by considering the Taylor expansion centered at (0, 0, 0).

In accordance to what was proven in the 2D setting, the minimal eigenvalue of the symbol *f*Q<sup>1</sup> (*θ*1, *θ*2, *θ*3) has a unique zero of order two at (*θ*1, *θ*2, *θ*3)=(0, 0, 0). In fact, the symbol *f*Q<sup>1</sup> is positive semi-definite, and all the eigenvalues of the symbol *f*Q<sup>1</sup> (*θ*1, *θ*2, *θ*3) have a unique zero of order two at (*θ*1, *θ*2, *θ*3)=(0, 0, 0), that is, there exist positive constants *c*(*j*), *C*(*j*), *j* = 1, 2, 3 such that

$$\mathcal{L}(j) \| (\theta\_1, \theta\_2, \theta\_3)^T \|^{\
\|} \le \lambda\_j (f\_{\mathbb{Q}\_1}(\theta\_1, \theta\_2, \theta\_3)) \le \mathcal{C}(j) \| (\theta\_1, \theta\_2, \theta\_3)^T \|^{\
\|}$$

uniformly in a proper neighborhood of (0, 0, 0).

Therefore, we can infer important information on the extremal eigenvalues and on the conditioning of the corresponding matrix sequences.

**Proposition 4.** *Let f*Q<sup>1</sup> (*θ*1, *θ*2, *θ*3) *be the symbol defined in* (5)*. Let us consider a uniform meshing in all the directions with n subintervals. Then we have the following relationships:*

$$\begin{aligned} A\_{\min}(A\_{\mathfrak{n}}(1, \mathbf{D}^{\delta})) &\sim n^{-2}, \\ \max f\_{\mathbb{Q}\_{1}} -\lambda\_{\max}(A\_{\mathfrak{n}}(1, \mathbf{D}^{\delta})) &\sim n^{-2}, \\ \mu(A\_{\mathfrak{n}}(1, \mathbf{D}^{\delta})) &\sim n^{2}, \\ A\_{\mathfrak{n}}(1, \mathbf{D}^{\delta}) &= T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_{1}}), \ \mathfrak{n} = (n\_{1}, n\_{2}, n\_{3}), \ \mathfrak{n}\_{1} = \mathfrak{n}\_{2} = n\_{3} = n - 1. \end{aligned}$$

*and*

$$\begin{aligned} &\lambda\_{\text{min}}(A\_{\text{il}}(1,\text{DN}^{5})) \sim n^{-2}, \\ &\max f\_{\mathbb{Q}\_{1}} - \lambda\_{\text{max}}(A\_{\text{il}}(1,\text{DN}^{5})) \sim n^{-2}, \\ &\mu(A\_{\text{il}}(1,\text{DN}^{5})) \sim n^{2}, \\ &A\_{\text{il}}(1,\text{DN}^{5}) = T\_{\text{n}}(f\_{\mathbb{Q}\_{1}}) + R\_{\text{n}}, \quad \mathfrak{n} = (n\_{1},n\_{2},n\_{3}), \; n\_{1} = n\_{\text{r}} \\ &\text{ }n\_{\text{s}} = n\_{\text{s}} = n\_{\text{s}} = n+1, \; \mu(A\_{\text{l}}(1,\text{DN}^{5})) \sim n^{-2}, \; \mu(A\_{\text{l}}(1,\text{DN}^{5})) \sim n^{-2}, \\ \end{aligned}$$

*with R<sup>n</sup> as in Proposition 2 and μ*(·) *denoting the condition number in spectral norm (the one induced by the Euclidean vector norm).*

*If the more general setting is considered with n*<sup>1</sup> *subintervals in the x*<sup>1</sup> *direction, n*<sup>2</sup> *subintervals in the x*<sup>2</sup> *direction, and n*<sup>3</sup> *subintervals in the x*<sup>3</sup> *direction, then analogous relations are true:*

$$\begin{aligned} \lambda\_{\min}(A\_{\mathfrak{n}}(1, \mathbf{D}^{\mathfrak{f}})) &\sim [\min n\_{j}]^{-2}, \\ \max f\_{\mathbb{Q}\_{1}} - \lambda\_{\max}(A\_{\mathfrak{n}}(1, \mathbf{D}^{\mathfrak{f}})) &\sim [\min n\_{j}]^{-2}, \\ \mu(A\_{\mathfrak{n}}(1, \mathbf{D}^{\mathfrak{f}})) &\sim n^{2}, \\ A\_{\mathfrak{n}}(1, \mathbf{D}^{\mathfrak{f}}) &= T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_{1}}), \ \mathfrak{n} = (n\_{1} - 1, n\_{2} - 1, n\_{3} - 1) \end{aligned}$$

*and*

$$\begin{aligned} \lambda\_{\text{min}}(A\_{\mathfrak{n}}(1, \text{DN}^5)) &\sim [\min n\_{j}]^{-2}, \\ \max f\_{\mathbb{Q}\_1} - \lambda\_{\text{max}}(A\_{\mathfrak{n}}(1, \text{DN}^5)) &\sim [\min n\_{j}]^{-2}, \\ \mu(A\_{\mathfrak{n}}(1, \text{DN}^5)) &\sim [\min n\_{j}]^2, \\ A\_{\mathfrak{n}}(1, \text{DN}^5) &= T\_{\mathfrak{n}}(f\_{\mathbb{Q}\_1}) + R\_{\mathfrak{n}}, \quad \mathfrak{n} = (n\_1, n\_2 + 1, n\_3 + 1), \end{aligned}$$

*with Rn as in Proposition 2.*

#### *3.3. Non-Constant Coefficient ρ Case—3D*

It should be observed that the natural extension of the previous analysis refers to the case of a non-constant coefficient *ρ*. According to a standard assembling procedure in FEs, given the triangulation T made by all the basic elements *τ*, we write the stiffness matrix as

$$A\_{\mathfrak{n}}(\rho) = \sum\_{\mathfrak{r} \in \mathcal{T}} \rho\_{\mathfrak{r}} A\_{\mathfrak{n}, \mathfrak{r}\prime}^{El} \tag{6}$$

where *AEl <sup>n</sup>*,*<sup>τ</sup>* is the elementary matrix *Kn* in (2) (possibly properly cut when nodes on the boundary are involved), but widened to size *N*(*n*) according to the chosen global ordering of nodes. Here, *n* = (*n* − 1, *n* − 1, *n* − 1) in the case of the Dirichlet boundary conditions and *n* subintervals in all the directions such that *An*(*ρ*) = *An*(*ρ*, D6), while *n* = (*n*, *n* + 1, *n* + 1) in the case of Dirichlet boundary conditions on one-side Neumann boundary conditions in the remaining ones with *n* subintervals in all the directions so that *An*(*ρ*) = *An*(*ρ*, DN5).

Clearly, the elementary matrix *Ke* is positive semi-definite for 0 ≤ *ν* < 1/2. More precisely, the not zero eigenvalues are

$$\begin{array}{llll} \text{1 (double)}, & \frac{(1-\nu)}{3(1-2\nu)} \text{ (triple)}, & \frac{2\nu}{(1-2\nu)} \text{ (triple)}, & \frac{\nu}{3(1-2\nu)} \text{ (double)},\\ \frac{4\nu}{3(1-2\nu)} \text{ (simple)}, & \frac{(\nu+1)}{(1-2\nu)} \text{ (simple)}, & \frac{(\nu+1)}{3(1-2\nu)} \text{ (triple)}, & \frac{(\nu+1)}{9(1-2\nu)} \text{ (triple)}. \end{array}$$

In addition, every elementary matrix, which has been cut due to nodes on the boundary, is positive semi-definite as well since it is a principal submatrix of *Ke* in (2).

Thus, on the basis of (6) and by applying the Courant–Fisher theorem, we can claim again

$$\begin{aligned} \rho\_{\min} \lambda\_{\min}(A\_n(1)) \le \lambda\_{\min}(A\_n(\rho)) \le \rho\_{\max} \lambda\_{\min}(A\_n(1)),\\ \rho\_{\min} \lambda\_{\max}(A\_n(1)) \le \lambda\_{\max}(A\_n(\rho)) \le \rho\_{\max} \lambda\_{\max}(A\_n(1)). \end{aligned}$$

with *ρ*min and *ρ*max minimum and maximum of *ρ*, respectively (see Figure 1 for an upperbound of the maximal eigenvalues of *An*(1) as a function of *ν*).

**Figure 1.** Maxima of eigenvalue surfaces of the symbol *f*Q<sup>1</sup> vs. *ν*.

By combining the above inequalities and Proposition 4, we infer that the extremal eigenvalues and the conditioning have the same asymptotical behavior as in the constant coefficient setting.

The whole eigenvalues distribution is sketched below by referring to the basics of the GLT theory [29] reported in Section 3.1.2. Let *Dn*(*ρ*) be a multilevel block diagonal sampling matrix according to the notions introduced in Section 3.1.2 and let *An*(1) be the multilevel block Toeplitz matrix *Tn*(*f*Q<sup>1</sup> ) if the Dirichlet boundary conditions are used or *Tn*(*f*Q<sup>1</sup> ) + *R<sup>n</sup>* in the other case. Then the following facts hold:

**Fact 1** {*Dn*(*ρ*)}*<sup>n</sup>* ∼GLT *ρ* according to Item **GLT 2.**

**Fact 2** {*Rn*}*<sup>n</sup>* ∼GLT 0 according to Proposition 2, Proposition 1, and Item **GLT 2.**


sequence {*An*(*ρ*)}*<sup>n</sup>* is made up by Hermitian matrices, by Item **GLT 1.** it follows that {*An*(*ρ*)}*<sup>n</sup>* ∼*σ*,*<sup>λ</sup> ρ f*Q<sup>1</sup> .

We stress that a unique notation has been employed for the sake of notational compactness. However, the main statement, Fact 6, holds for all the matrix sequences reported in Proposition 3 in the case of a variable *ρ*.

Finally, it is worth noticing that non-square domains can be treated as well, using the reduced GLT theory (see pages 398–399 in [33], Section 3.1.4 in [34], and the recent analysis [35]): here, we do not go in the details and this will be the subject of future investigations.

#### **4. Two-Grid and Multigrid Methods**

In this section, we concisely report few relevant results concerning the convergence theory of algebraic multigrid methods with special attention of the case of multilevel block Toeplitz structures generated by a matrix-valued symbol *f* .

We start by taking into consideration the generic linear system *Amxm* = *bm* with large dimension *<sup>m</sup>*, where *Am* <sup>∈</sup> <sup>C</sup>*m*×*<sup>m</sup>* is a Hermitian positive definite matrix and *xm*, *bm* <sup>∈</sup> <sup>C</sup>*m*. Let *<sup>m</sup>*<sup>0</sup> <sup>=</sup> *<sup>m</sup>* <sup>&</sup>gt; *<sup>m</sup>*<sup>1</sup> <sup>&</sup>gt; ... <sup>&</sup>gt; *ms* <sup>&</sup>gt; ... <sup>&</sup>gt; *ms*min and let *<sup>P</sup>s*+<sup>1</sup> *<sup>s</sup>* <sup>∈</sup> <sup>C</sup>*ms*+1×*ms* be a full-rank matrix for any *s*. At last, let us denote by V*<sup>s</sup>* a class of stationary iterative methods for given linear systems of dimension *ms*.

In accordance with [36], the algebraic two-grid method (TGM) can be easily seen as a stationary iterative method whose generic steps are reported below.

$$\mathbf{x}\_s^{\rm out} = \mathbf{T} \mathbf{G} \mathbf{M} (\mathbf{s}\_\prime \mathbf{x}\_s^{\rm in}, b\_s)$$

*x* pre *<sup>s</sup>* <sup>=</sup> <sup>V</sup>*ν*pre *s*,pre(*x*in *<sup>s</sup>* , *bs*) Pre-smoothing iterations *rs* = *Asx* pre *<sup>s</sup>* − *bs rs*+<sup>1</sup> <sup>=</sup> *<sup>P</sup>ms*+<sup>1</sup> *ms rs As*+<sup>1</sup> <sup>=</sup> *<sup>P</sup>ms*+<sup>1</sup> *ms As*(*Pms*+<sup>1</sup> *ms* )*<sup>H</sup>* Solve *As*+1*ys*+<sup>1</sup> = *rs*+<sup>1</sup> *x*ˆ*<sup>s</sup>* = *x* pre *<sup>s</sup>* <sup>−</sup> (*Pms*+<sup>1</sup> *ms* )*Hys*<sup>+</sup><sup>1</sup> Exact coarse grid correction (CGC) *x*out *<sup>s</sup>* <sup>=</sup> <sup>V</sup>*ν*post *<sup>s</sup>*,post(*x*ˆ*s*, *bs*) Post-smoothing iterations

Here, we refer to the dimension *ms* by means of its subscript *s*.

In the first and last steps, a *pre-smoothing iteration* and a *post-smoothing iteration* are applied *ν*pre times and *ν*post times, respectively, taking into account the considered stationary iterative method in the class V*s*. Furthermore, the intermediate steps define the *exact coarse grid correction operator*, which is dependent on the considered projector operator *P<sup>s</sup> <sup>s</sup>*+1. The resulting iteration matrix of the TGM is then defined as

$$\begin{aligned} TGM\_s &= V\_{s, \text{post}}^{\text{"post"}} \mathcal{C} \mathcal{C} \mathcal{C}\_s \mathcal{V}\_{s, \text{pre}}^{\text{"pre"}}\\ \mathcal{C} \mathcal{G} \mathcal{C}\_s &= I^{(s)} - (P\_{m\_s}^{m\_{s+1}})^H A\_{s+1}^{-1} P\_{m\_s}^{m\_{s+1}} A\_s\\ A\_{s+1} &= P\_{m\_s}^{m\_{s+1}} A\_s (P\_{m\_s}^{m\_{s+1}})^H \end{aligned}$$

where *Vs*,pre and *Vs*,post represent the pre-smoothing and post-smoothing iteration matrices, respectively, and *I*(*s*) is the identity matrix at the *s*th level.

By employing a recursive procedure, the TGM leads to a multigrid method (MGM). Indeed, the standard V-cycle can be expressed in the following way:

From a computational viewpoint, to reduce the related costs, it is more efficient that the matrices *As*+<sup>1</sup> = *Ps*+<sup>1</sup> *<sup>s</sup> As*(*Ps*+<sup>1</sup> *<sup>s</sup>* )*<sup>H</sup>* are computed in the so-called *setup phase*.

According to the previous setting, the global iteration matrix of the MGM is recursively defined as

$$\begin{aligned} MGM\_{s\_{\min}} = O &\in \mathbb{C}^{s\_{\min} \times s\_{\min}},\\ MGM\_{s} = V\_{s\_{\text{post}}}^{\text{Vpsot}} \Big[ I^{(s)} - (P\_{m\_s}^{m\_{s+1}})^H \Big( I^{(s+1)} - MGM\_{s+1} \Big) A\_{s+1}^{-1} P\_{m\_s}^{m\_{s+1}} A\_s \Big] V\_{s, \text{pre}}^{\text{Vpsor}} \\ s &= s\_{\min} - 1, \dots, 0. \end{aligned}$$

Lastly, the W-cycle is just a variation of the previous V-cycle considering two recursive calls in the coarse grid correction as follows:

$$\mathfrak{x}\_{\mathrm{s}}^{\mathrm{out}} = \mathsf{M}\mathbf{G}\mathsf{M}\mathsf{W}(\mathsf{s}, \mathfrak{x}\_{\mathrm{s}}^{\mathrm{in}}, b\_{\mathrm{s}}),$$

In the following remark, we emphasize the relevant computational properties of the considered methods.

**Remark 1.** *The first part of the current remark concerns the computational cost. In the V-cycle, there is just one recursive call, while in the W-cycle, there are two recursive calls, which in principle, is more expensive. Let us analyze the related costs. Since our matrices are sparse and they remain*

*sparse at the lower levels, the number of arithmetic operations in the V-cycle and W-cycle algorithms is of the type*

$$\mathbb{C}\_V(m) \le \mathbb{C}\_V(m/q) + \alpha m, \qquad \qquad \mathbb{C}\_W(m) \le 2\mathbb{C}\_W(m/q) + \alpha m, \dots$$

*where m* = *m*<sup>0</sup> *and ms*+<sup>1</sup> = *ms*/*q. Now in the one-dimensional setting q* = 2 *and CV*(*m*) = *O*(*m*)*, while CW*(*m*) = *O*(*m* log *m*) *and thus the W-cycle is asymptotically more expensive. Conversely, in a d-dimensional setting with <sup>d</sup>* <sup>≥</sup> <sup>2</sup>*, we have <sup>q</sup>* <sup>=</sup> <sup>2</sup>*<sup>d</sup> and, as a consequence of the recursive relations, both CV*(*m*) *and CW*(*m*) *are linear in the matrix-size m. To be precise, we have*

$$\mathbb{C}\_V(m) \le \frac{2^d}{2^d - 1} \alpha m\_\prime \qquad \qquad \mathbb{C}\_W(m) \le \frac{2^{d-1}}{2^{d-1} - 1} \alpha m\_\prime$$

*so that, in the current three-dimensional setting, the bounds of the V-cycle and W-cycle complexity are very close, i.e.,*

$$\mathbb{C}\_V(m) \le \frac{8}{7} \alpha m\_\prime \qquad\qquad \mathbb{C}\_W(m) \le \frac{4}{3} \alpha m.$$

*Another issue to be considered is the case where the discretization matrices appear in checkerboard fashion (see, for example, [37] and the references therein). We emphasize that our analysis, which is based on the Toeplitz generating function, is not directly applicable since we lose the Toeplitz character of the approximation matrices when a checkerboard (called also red–black) ordering is used. However, this is just a matter of a similarity transformation by a permutation matrix. Hence, with a careful work and without increasing the complexity analyzed in the previous lines, the algorithm discussed here and the related theoretical analysis can be adapted to the checkerboard context.*

**Remark 2.** *In the relevant literature (see, for instance, [38]), the convergence analysis of the TGM splits into checking two separate conditions: the smoothing property and the approximation property. Regarding the latter and regarding scalar structured matrices [38,39], the TGM optimality is given in terms of choosing the proper conditions that the symbol p of a family of projection operators has to fulfill. Indeed, let Tn*(*f*) *with <sup>n</sup>* = (2*<sup>t</sup>* <sup>−</sup> <sup>1</sup>)*, where <sup>f</sup> is a non-negative trigonometric polynomial. Let θ*<sup>0</sup> *be the unique zero of f . Then the TGM optimality applied to Tn*(*f*) *is guaranteed if we choose the symbol p of the family of projection operators such that*

$$\begin{aligned} \limsup\_{\theta \to \theta^0} \frac{|p(\eta)|^2}{f(\theta)} &< \infty, \ \eta \in \mathcal{M}(\theta), \\ \sum\_{\eta \in \Omega(\theta)} |p(\eta)|^2 &> 0, \end{aligned} \tag{7}$$

*where, for d* = 1*, the sets* Ω(*θ*) *and* M(*θ*) *are the following corner and mirror points*

$$\Omega(\theta) = \{ \eta \in \{\theta, \theta + \pi\} \}, \qquad \mathcal{M}(\theta) = \Omega(\theta) \backslash \{\theta\}.$$

*respectively. In the general case of d* > 1*, we have*

$$\Omega(\theta) = \{ \eta \in \{ \theta + \tau \mathbf{s} \}, \mathbf{s} = (\mathbf{s}\_1, \dots, \mathbf{s}\_d), \mathbf{s}\_j \in \{ 0, 1 \}, j = 1, \dots, d \}$$

*with* <sup>M</sup>(*θ*) = <sup>Ω</sup>(*θ*) \ {*θ*}*, so that the cardinality of* <sup>Ω</sup>(*θ*) *and* <sup>M</sup>(*θ*) *is* <sup>2</sup>*<sup>d</sup> and* <sup>2</sup>*<sup>d</sup>* <sup>−</sup> <sup>1</sup>*, respectively.*

Informally, for *d* = 1, it means that the optimality of the two-grid method is obtained by choosing the family of projection operators associated to a symbol *p* such that |*p*| <sup>2</sup>(*ϑ*) + <sup>|</sup>*p*<sup>|</sup> <sup>2</sup>(*<sup>ϑ</sup>* <sup>+</sup> *<sup>π</sup>*) does not have zeros and <sup>|</sup>*p*<sup>|</sup> <sup>2</sup>(*ϑ* + *π*)/ *f*(*ϑ*) is bounded (if we require the optimality of the V-cycle, then the second condition is a bit stronger); see [38]. For approximation of an elliptic differential operator of order 2*α* by local methods (e.g., finite differences, FEs, and isogeometric analysis), the previous conditions mean that *p* has a unique zero of order at least *α* at *ϑ* = *π* whenever *f* has a unique zero at *θ*<sup>0</sup> = 0 of order 2*α*.

In our specific block setting, by interpreting the analysis given in [40], all the involved symbols are matrix valued, and the conditions which generalize (7) and are sufficient for the TGM convergence and optimality are the following:


Even if the theoretical extension to the V-cycle and W-cycle convergence and optimality is not given, in the subsequent section, we propose specific choices of the projection operators numerically showing how this leads to two-grid, V-cycle, and W-cycle procedures converging optimally or quasi-optimally with respect to all the relevant parameters (size, dimensionality, and polynomial degree *k*).

Our choices are in agreement with the mathematical conditions set in items **(A)**, **(B)**, and **(C)**. We remark that the violation of condition **C)** is discussed by Donatelli et al. [40], in their conclusion section.

#### **5. Multigrid Proposals**

Based on the theory in Section 4, we define the prolongation operator as

$$P\_{2h}^h = P \otimes P \otimes P \otimes I\_3 / \sqrt{2}$$

in the case of matrices *An*(1, D6), and

$$P\_{2h}^{\text{lt}} = P\_{\text{c}} \otimes P\_{\text{t}} \otimes P\_{\text{t}} \otimes I\_3 / \sqrt{2}$$

in the case of matrices *An*(1, DN5). In the expressions above,

$$P = \begin{bmatrix} \frac{1}{2} \\ 1 \\ \frac{1}{2} & \frac{1}{2} \\ & 1 \\ & 1 \\ & \frac{1}{2} & \ddots & \frac{1}{2} \\ & & \frac{1}{2} \\ & & & \frac{1}{2} \end{bmatrix}, P\_t = \begin{bmatrix} 1 \\ \frac{1}{2} & \frac{1}{2} \\ & 1 \\ & \frac{1}{2} & \frac{1}{2} \\ & & 1 \\ & & \frac{1}{2} & \ddots & \frac{1}{2} \\ & & & 1 \\ & & & \frac{1}{2} & \frac{1}{2} \\ & & & & \frac{1}{2} & \frac{1}{2} \\ & & & & & 1 \end{bmatrix}, \text{and } P\_t = \begin{bmatrix} \frac{1}{2} \\ 1 \\ 1 \\ \frac{1}{2} \\ 1 \\ \frac{1}{2} \\ \vdots \\ \frac{1}{2} \end{bmatrix}.$$

Note that *A*2*<sup>h</sup>* = *PTAhP* equals the same FEM approximation on the coarse mesh in both cases, independently of the boundary conditions.

Independently of the boundary conditions, we observe that the symbol of the prolongation operator is a 3 × 3 matrix-valued function as the symbol of the coefficient matrix sequence of the linear systems to be solved: more precisely we define

$$p(\theta\_1, \theta\_2, \theta\_3) = \sqrt{2^{-1}} (1 + \cos(\theta\_1)(1 + \cos(\theta\_2)(1 + \cos(\theta\_3))),$$

Since the matrix-valued function *p*(*θ*1, *θ*2, *θ*3) is a scalar function times the identity, it automatically follows that the difficult condition **(C)** is satisfied. Moreover, taking into account Theorem 1, since the scalar function 1 + cos(*θ*) has a zero of order 2 at *θ* = *π*, also condition **(A)** holds, while the positive definiteness of <sup>∑</sup>*η*∈Ω(*θ*) *ppH*(*η*) is verified as well, so that also condition **(B)** is met.

Finally, we choose one iteration of Gauss–Seidel as the pre and post smoother. Given the analysis reported in the previous section, we expect that our multigrid method is convergent in an optimal way, that is, with a convergence speed independent of the matrix size. In Table 1, we find a plain confirmation of our theoretical expectations, which holds not only for a constant *ρ*, but also in the variable coefficient setting as long as *ρ* remains positive.


**Table 1.** Number of iteration for matrices *An*(1, D6) and *An*(1, DN5) of increasing dimension *N*(*n*), with the numeric precision *ε* = 10−4, 10−6, 10−8, respectively.

#### **6. Conclusions**

We provided a quite complete spectral analysis of the (large) coefficient matrices associated with the linear systems stemming from the FE discretization of a linearly elastic problem for an arbitrary element-wise constant coefficient field. Our interest in this problem stems from the fact that in the material distribution method for topology optimization, such a problem is solved at each iteration. The solution for these linear systems typically dominates the computational effort required to solve the topology optimization problems. Based on the spectral information, we proposed a specialized multigrid method, which turned out to be optimal in the sense that the (arithmetic) cost for solving the related linear systems, up to a fixed desired accuracy, is proportional to the computational cost of matrix–vector products, which is linear in the corresponding matrix size and mildly depending on the given accuracy. The method was tested, and the preliminary numerical results are very satisfactory and promising, in terms of having a linear cost and a number of iterations that is bounded by a constant independent of the matrix size and only lightly influenced by the desired accuracy.

Finally, we mention future lines of research:


of high order and intermediate regularity (see [30,31] and references therein). Hence, we are convinced that such extensions can be formally obtained, also in the case of the current topology optimization problem.

**Author Contributions:** Conceptualization, Q.K.N., S.S.-C., C.T.-P. and E.W.; Investigation, Q.K.N., S.S.-C., C.T.-P. and E.W.; Methodology, Q.K.N., S.S.-C., C.T.-P. and E.W.; Writing—original draft, Q.K.N., S.S.-C., C.T.-P. and E.W.; Writing—review and editing, Q.K.N., S.S.-C., C.T.-P. and E.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is financially supported by the Swedish strategic research program eSSENCE and the Italian Institution for High Mathematics (INdAM).

**Data Availability Statement:** The data in this paper are results of numerical computations and not experimental measurements. The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **Appendix A. Explicit Expressions for the Element Stiffness Matrix**

Recall that for the studied problem, the three-dimensional element stiffness matrix is

$$
\mathbb{K}\_{n} = \frac{h}{2} \frac{E\_0}{1+\nu} \begin{bmatrix}
\begin{matrix}
K\_1 & K\_2 & K\_3 & K\_4 & K\_5 & K\_6 & K\_7 & K\_8 \\
& K\_9 & K\_{10} & K\_{11} & K\_{12} & K\_{13} & K\_{14} & K\_7 \\
& & & K\_{15} & K\_{16} & K\_{17} & K\_{18} & K\_{13} & K\_6 \\
& & & & K\_{19} & K\_{20} & K\_{17} & K\_{12} & K\_5 \\
& & & & & K\_{19} & K\_{16} & K\_{11} & K\_4 \\
& & & & & & K\_{15} & K\_{10} & K\_3 \\
& & & & & & & K\_9 & K\_2 \\
& & & & & & & K\_1
\end{matrix}
\end{bmatrix}
$$

where

$$\begin{aligned} \begin{array}{c} K\_{1} = \begin{bmatrix} k\_{1} \ k\_{2} \ k\_{2} \\ k\_{1} \ k\_{2} \end{bmatrix}, \quad K\_{2} = \begin{bmatrix} k\_{3} \ k\_{4} \ 0 \\ k\_{3} \ 0 \\ k\_{6} \end{bmatrix}, \quad K\_{3} = \begin{bmatrix} k\_{3} \ 0 \ 0 \\ k\_{6} \ 0 \\ k\_{7} \end{bmatrix}, \quad K\_{4} = \begin{bmatrix} k\_{5} \ 0 \ 0 \\ k\_{7} \ 0 \\ k\_{8} \end{bmatrix}, \\\ K\_{5} = \begin{bmatrix} k\_{6} \ 0 \ 0 \\ k\_{5} \ 0 \\ k\_{3} \end{bmatrix}, \quad K\_{6} = \begin{bmatrix} k\_{7} \ 0 \ 0 \\ k\_{5} \ 0 \\ k\_{7} \end{bmatrix}, \quad K\_{7} = \begin{bmatrix} k\_{7} \ 0 \ 0 \\ k\_{7} \ 0 \\ k\_{7} \end{bmatrix}, \quad K\_{8} = \begin{bmatrix} k\_{9} \ k\_{10} \ 0 \\ k\_{9} \ 0 \\ k\_{9} \end{bmatrix}, \\\ K\_{9} = \begin{bmatrix} k\_{1} \ \ k\_{2} \ \kappa \\ k\_{1} \ \kappa \\ k\_{1} \end{bmatrix}, \quad K\_{10} = \begin{bmatrix} k\_{5} \ 0 \ 0 \\ k\_{7} \ \kappa \\ k\_{7} \end{bmatrix}, \quad K\_{11} = \begin{bmatrix} k\_{9} \ 0 \ 0 \\ k\_{6} \ 0 \\ k\_{9} \end{bmatrix}, \quad K\_{12} = \begin{bmatrix} k\_{9} \ 0 \ k\_{10} \\ k\_{5} \ 0 \\ k\_{7} \end{bmatrix}, \quad K\_{13} = \begin{bmatrix} k\_{7} \ 0 \ 0 \\ k\_{5} \ 0 \\ k\_{7} \end{bmatrix}, \quad K\_{12} = \begin{bmatrix} k\_{1} \$$

We stress that all the 3 × 3 blocks mentioned before are real symmetric and hence the lower part is defined accordingly; see Section 3.2.1 for the use of these blocks in the definition of the GLT symbol. We observe that the quantity *hE*0/2(1 + *ν*) is a multiplicative term that will be neglected in the spectral analysis since it will be simplified and included in the right-hand vector when solving the related large linear systems). An analytical evaluation of the entries of the stiffness matrix is reported in the following lines:

$$\begin{aligned} k\_1 &= \frac{2\dot{\lambda}}{3} + \frac{4\dot{\mu}}{9}, & k\_2 &= \frac{\dot{\lambda}}{3}, & k\_3 &= \frac{2\dot{\mu}}{9}, & k\_4 &= \frac{\dot{\lambda}}{6}, & k\_5 &= \frac{\dot{\mu}}{9} - \frac{\dot{\lambda}}{6}, \\ k\_6 &= \frac{-4\dot{\mu}}{9}, & k\gamma &= \frac{-\dot{\lambda}}{6} - \frac{2\dot{\mu}}{9}, & k\mathfrak{s} &= -\frac{\dot{\lambda}}{3}, & k\mathfrak{s} &= -\frac{\dot{\lambda}}{6} - \frac{\dot{\mu}}{9}, & k\_{10} &= -\frac{\dot{\lambda}}{6}, \end{aligned} \tag{A2}$$

where

$$
\lambda = \frac{\nu}{(1 - 2\nu)}, \quad \hat{\mu} = \frac{1}{2}.
$$

Note that *λ*ˆ and *μ*ˆ in (A2) do not represent the Lamé parameters because they have been modified to facilitate the spectral analysis conducted in Section 3. In fact, in this article, *λ* = *λ*ˆ *<sup>E</sup>*<sup>0</sup> <sup>1</sup>+*<sup>ν</sup>* and *<sup>μ</sup>* <sup>=</sup> *<sup>μ</sup>*<sup>ˆ</sup> *<sup>E</sup>*<sup>0</sup> <sup>1</sup>+*<sup>ν</sup>* are the actual Lamé parameters.

#### **Appendix B. Stress–Strain Relation and Various Bounds**

In the three-dimensional setting, there are six independent strain components in total at a point in an element, and they are written as a vector

$$\boldsymbol{\epsilon} = [\boldsymbol{\epsilon}\_{11} \ \boldsymbol{\epsilon}\_{22} \ \boldsymbol{\epsilon}\_{33} \ 2\boldsymbol{\epsilon}\_{12} \ 2\boldsymbol{\epsilon}\_{23} \ 2\boldsymbol{\epsilon}\_{31}]^T.$$

Similarly, corresponding to the six strain components above, there are also six independent stress components written in vector form as

$$
\sigma = [\sigma\_{11} \; \sigma\_{22} \; \; \sigma\_{33} \; \; \; \sigma\_{12} \; \; \; \sigma\_{23} \; \; \; \sigma\_{31}]^T \; .
$$

By the generalized Hooke's law, the most general linear relation among components of the stress and strain tensor can then be written as

$$
\sigma = \mathsf{E}\mathfrak{e},
\tag{\mathsf{A3}}
$$

where *E* is a matrix that corresponds to the constant fourth-order elasticity tensor *Ec*. The relationship between stresses and strains is

$$\begin{aligned} \varepsilon\_{11} &= \frac{1}{E\_0} (\sigma\_{11} - \nu (\sigma\_{22} + \sigma\_{33})), & \varepsilon\_{12} &= \frac{\sigma\_{12}}{2G}, & \varepsilon\_{13} &= \frac{\sigma\_{13}}{2G}, \\ \varepsilon\_{22} &= \frac{1}{E\_0} (\sigma\_{22} - \nu (\sigma\_{11} + \sigma\_{33})), & \varepsilon\_{23} &= \frac{\sigma\_{23}}{2G}, & \varepsilon\_{33} &= \frac{1}{E\_0} (\sigma\_{33} - \nu (\sigma\_{11} + \sigma\_{22})), \end{aligned} \tag{A4}$$

where *ν* is Poisson's ratio, *E*<sup>0</sup> is Young's modulus, and the shear modulus *G* is

$$G = \frac{E\_0}{2 + 2\nu}.\tag{A5}$$

The relationships above can be rewritten in matrix form as

$$
\mathfrak{e} = \frac{1}{E\_0} \begin{bmatrix}
1 & -\nu & -\nu & 0 & 0 & 0 \\
0 & 0 & 0 & 2(1+\nu) & 0 & 0 \\
0 & 0 & 0 & 0 & 2(1+\nu) & 0 \\
0 & 0 & 0 & 0 & 0 & 2(1+\nu)
\end{bmatrix} \sigma. \tag{A6}
$$

Furthermore, the equation system (A6) can be inverted to obtain Hooke's law (A3) with

$$E = \begin{bmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0\\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0\\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0\\ 0 & 0 & 0 & \lambda & 0 & 0\\ 0 & 0 & 0 & 0 & \lambda & 0\\ 0 & 0 & 0 & 0 & 0 & \lambda \end{bmatrix} \tag{A7}$$

where

$$\begin{aligned} \lambda &= \frac{E\_0 \nu}{(1+\nu)(1-2\nu)},\\ \mu &= \frac{E\_0}{2(1+\nu)}. \end{aligned} \tag{A8}$$

In this case, the bulk modulus *K* can be expressed as

$$K = \frac{(\sigma\_{11} + \sigma\_{22} + \sigma\_{33})/3}{\epsilon\_{11} + \epsilon\_{22} + \epsilon\_{33}} = \frac{E\_0}{3(1 - 2\nu)}. \tag{A9}$$

The Young's (*E*0), shear (*G*), and bulk (*K*) moduli need to be positive. Thus, Equations (A5) and (A9) imply that the Poisson's ratio in three dimensions must satisfy −1 < *ν* < 0.5.

**Remark A1.** *Physically, there are no known isotropic materials with a negative Poisson's ratio. Therefore, from a practical perspective [43], we can limit our study to ν, satisfying the inequalities* 0 ≤ *ν* < 0.5*. Furthermore, in some more delicate circumstances, the Poisson ratio lies in the interval* 0.2 ≤ *ν* < 0.5*, as proved experimentally in view of the elastic properties of real materials [44].*

#### **References**

