Spectral analysis of the finite element matrices approximating 2D linearly elastic structures and multigrid proposals

Topology optimization aims to find the best material layout subject to given constraints. The so‐called material distribution methods cast the governing equation as an extended or fictitious domain problem, in which a coefficient field represents the design. When solving the governing equation using the finite element method, a large number of elements are used to discretize the design domain, and an element‐wise constant function approximates the coefficient field in the considered design domain. This article presents a spectral analysis of the (large) coefficient matrices associated with the linear systems stemming from the finite element discretization of a linearly elastic problem for an arbitrary coefficient field. Based on the spectral information, we design a multigrid method which turns 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 matrix‐vector cost, which is linear in the corresponding matrix size. The method is tested, and the numerical results are very satisfactory in terms of linear cost and number of iterations, which is bounded by a constant independent of the matrix size.

conditions. For example, the corresponding technique is used for designing advanced lightweight components in the car and aeronautical industries. 13 In the material distribution setting, the decision variable in the optimization problem is a design function ∶ Ω → {0, 1}, that defines the physical density .
Typically, topology optimization problems are solved numerically by discretizing Ω into N elements and by an element-wise constant function with element values = [ 1 , 2 , … , N ]. After discretization, the use of a binary-valued material indicator function is not feasible for many reasons, one is that the problem becomes computationally intractable. A standard solution for this issue is to allow the element values i to take values in the range [0, 1] and to define the physical density by using penalization and filtering. In particular, the physical density is defined as (x) = + ( solid − )g ( ( )(x)), where  is a filter operator, g is a penalty function, and ≥ 0 is a constant. The relaxation of the material indicator function enables the use of gradient-based optimization methods, which, together with the efficient adjoint-based computation of design sensitivities, empowers the efficient solution of large-scale problems. The standard approach to attain a unique solution for the state problem is to set to a small but strictly positive value. Berggren and Kasolis 14 studied a linearly elastic boundary value problem (BVP) defined on Ω and proved that the approximation error is bounded by the sum of three terms: a standard finite element (FE) approximation error term and two additional terms, both tending to zero as → 0 + . However, there are some drawbacks when setting to a small positive value, which have been indicated and discussed by Buhl et al. 15 and Bruns and Tortorelli. 16 To avoid those issues, there are few studies for vanishing lower bound. Bruns 17 let = 0 and treated the problem by using spectral decomposition and singular value decomposition (SD/SVD), to construct a generalized inverse (pseudoinverse) stiffness matrix K N for solving the FE system of equations. Nguyen et al. 18 have introduced a preconditioning approach that allows the element densities to take values from a continuum and use a specific ad hoc preconditioner when solving a standard test problem. However, this approach has been only proved analytically in one spatial dimension (d = 1), by using a suitable spectral analysis. In fact, there is a very limited number of studies that use spectral analysis for treating the considered classes of stiffness matrices {K N } N , arising from FE approximation in topology optimization. Nguyen et al. 18 have used a new spectral analysis method that is the so-called theory of Generalized Locally Toeplitz (GLT) sequences to compute and analyze the asymptotic spectral distribution of the stiffness matrices {K N } N for one spatial dimension problems.
The GLT sequences theory 19,20 is applied here 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 d N increases with N, and it tends to infinity as N → ∞. Hence, what needs to be considered is not just a single linear system, but an entire sequence of linear systems with increasing size. Under suitable conditions, the sequence of discretization matrices, such as {K N } N , has an asymptotic spectral distribution. More precisely, for a large set of test functions F (usually, for all continuous functions F with bounded support), it often happens that the following limit relation holds: where j (K N ), j = 1, … , d N are the eigenvalues of K N , t (D) ∈ (0, ∞) is the measure (t-dimensional volume) of D, and f ∶ D → C is the spectral symbol of the sequence {K N } N . The spectral symbol f contains spectral information briefly described informally as follows: assuming that N is large enough, the eigenvalues of K N , except possibly for a small number of outliers, are approximately equal to the samples of f over a uniform grid in D. It is then clear that the symbol f provides a "compact" and quite accurate description of the spectrum of the matrices K N (for N large enough), where, of course, D is related to the design domain Ω and t is related to the dimensionality d of Ω.
Here, we extend the results of our previous work 18 to the two-dimensional case (d = 2). A key component is the theory of multilevel block GLT sequences, 21,22 which provides the tools for computing the spectral distribution of block-structured matrices arising from the FE approximation in two-dimensional topology optimization, where the symbol f appearing in relation (1) is matrix-valued instead of scalar-valued, according to the concepts provided in Definition 1. By using this theory and more results, 23,24 we perform a detailed spectral analysis of the linear systems associated with the FE discretization of the governing equation. Moreover, the information obtained from the spectral symbol f is exploited to design a fast multigrid solver. More precisely, the proposed multigrid technique turns 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 the matrix-vector products, which is linear with respect to the corresponding matrix size. The method is tested and the numerical results are promising, in terms of linear cost and number of iterations, which is bounded by a constant independent of the matrix size and only mildly depending on the desired accuracy.
It is well-known that the spectral condition numbers of FE matrices grows as O(h −2 ), where h is the element size. This type of result is quite old (see Reference 25 and references therein): for instance, in some settings, it is enough to work with elementary tools such as the Gerschgorin theorems, at least for proving the positive definiteness. However, in our setting, the matrices have an expression, which is less standard, due to the tensor structure of the operator (u), and, in fact, the involved matrices do not possess the form described by Dorostkar et al. 26 In reality, if one applies the Gerschgorin theorems, the conclusions are less general and parameter ( ) dependent, as it is clear looking at the intricate expression of the matrix in Equation (6); refer also to (7). Hence we need more sophisticated mathematical instruments along the lines followed in Reference 26. Indeed, we observe that our results and those in Reference 26 are similar in the sense that the same block multilevel GLT machinery is employed, but different just because the matrices have different mathematical expressions and the derivation of the symbol has to be performed from scratch. We finally stress that determining the symbol is important for two additional reasons, which go beyond the conditioning analysis: • the position and the order of the zeros, not only for understanding the behavior of the spectral condition numbers but especially for designing the projection/restriction operators, with the target of obtaining optimal solvers, as done in Section 5; • the global distribution results for the eigenvalues, which of course are not classical and cannot be recovered using classical tools and, moreover, they have an impact in the modern analysis of the convergence speed of (preconditioned) Krylov methods; see References 27,28.
The article is organized as follows. In Section 2, we describe the continuous problem and the resulting coefficient matrices arising from our FE approximation. Section 3 is devoted to the spectral analysis of the FE matrices in the two-dimensional setting, from the perspective of the GLT theory. In Section 4, we give a brief account of multigrid methods, with special attention to the block case encountered in the present context. Section 5 contains a multigrid proposal based on the spectral information given in Section 3 and on the conditions reported in Section 4. Finally, conclusions are reported in Section 6. In addition, Appendix A is devoted to some relevant model information.

PRELIMINARIES AND DISCRETIZATION
In this section, we report the description of the continuous problem (Section 2.1), its approximation by a basic FE procedure (Section 2.2), and finally, the formal expression of the resulting FE matrices (Section 2.3). We emphasize that the formal expression of the relevant matrices is a key ingredient for applying the GLT theory, in order to produce a global spectral description of the matrix sequences under consideration.

A problem in linear elastostatics
A material that deforms under a load and resumes its undeformed shape, when removing the load, is called elastic. Provided that the load-induced deformations are small, many solids are linearly elastic, which means that the relation between the deformation and the applied load is linear. In this article, we consider a linearly elastic structure that (when unloaded) occupies the hyper-rectangular domain Ω ⊂ R d . Henceforth, let b ∈ L 2 (Ω) d be a given body load (a volume force) in Ω, t ∈ L 2 (Γ F ) d be the surface traction acting on the non-clamped boundary Γ F ⊂ Ω of the solid, and u denote the resulting equilibrium displacement. In linear elasticity, deformations are characterized by the so-called strain tensor We remark that the skew-symmetric part of the displacement gradient is related to rigid rotations (and not to the deformations) of the body. Hooke's generalized law states the relationship between the strain tensor and the so-called stress tensor that is with E being the fourth-order plane stress elasticity tensor having in total d 4 components. In this article, we limit our attention to the analysis of FE matrices arising in the two-dimensional plane stress setting, which is the standard in material distribution topology optimization problems. 1 By invoking symmetries, the number of independent components in E can be reduced 29(p. 194-197) . The equilibrium assumption, which states that the forces acting on any sub-body must be in balance, together with Cauchy's theorem, which states that the surface force density depends linearly on the normal derivative n, and the divergence theorem implies In this article, we assume that the structure is clamped along the boundary portion Γ D ⊂ Ω and subject to a surface traction load t on Γ F = Ω ⧵ Γ D . The boundary conditions above together with Equations (2) and (3) yield the BVP Here, the application in focus is material distribution based topology optimization. More precisely, we assume that the elasticity tensor is E c , where E c is a constant fourth-order elasticity tensor and ∈  ⊂ L ∞ (Ω). A typical choice is to where is a constant that satisfies 0 < ≪ 1. In the context of material distribution topology optimization, the finite element method (FEM) is the standard choice for generating numerical solutions of the BVP (4). This solution process uses the variational form of the BVP. Since the structure is clamped along the boundary portion Γ D ⊂ Ω, the set of all kinematically admissible displacements of the structure is Under the above assumptions, the steady-state displacement u of the structure is the solution to Here, the energy bilinear form a and the load linear form are defined as where 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.

Discretization by Q 1 finite elements
Let  h be a structured grid consisting of tensor product elements. More precisely, we define  h to be the space of all continuous functions that are zero on the boundary Γ D and 2-vectors with each component being linear in each coordinate direction on each element in  h . More precisely, we set where P 1 (S) is the space of all functions that are linear in each spatial component on element S. We define N = n 1 n 2 to be the number of elements in  h and M = m 1 m 2 to be the number of basis functions in  h where subscripts 1 and 2 denote the number of elements and the number of basis functions in x 1 and x 2 direction, respectively.
In order to obtain a discrete equation system corresponding to the variational problem (5), we approximate functions u ∈  and v ∈  by functions u h ∈  h and v h ∈  h , and approximate the physical density by a function h that is constant on each element in  h . We let u ∈ R M and v ∈ R M be the coefficient vectors of u h and v h , respectively. We define ∈ R N to be the vector whose entries coincide with the values of h in each element. In other words, we define  h ⊂ R N to be the set of vectors that correspond to element-wise constant functions h that are in . By applying the above approximations, we deduce that the variational problem (5) reduces to the linear system where the stiffness matrix K( ) and the right hand vector f have entries and

respectively.
A standard procedure for assembling the stiffness matrix K is to loop over each element so that where K n is the elementary stiffness matrix for the element S n . In the studied setup, we remark that all elements have the same shape and size, and hence all elements in the stiffness matrices possess the same non-zero values. The positions holding these values in the element matrices correspond to the indices of the basis function having support in the element.

Explicit expressions for the element stiffness matrix
In this article, we limit our attention to two spatial dimensions. Hence, if not otherwise stated, by default d = 2. In the subsequent lines, we give an explicit expression of the element stiffness matrix K n with respect to the reference element illustrated in Figure 1A, which will be assembled into stiffness matrix K using the global node numbering illustrated in Figure 1B. For example, for the first element, the nodes with local node numbers 5, 6, 7, and 8 have the global node numbers 2m 2 + 1, 2m 2 + 2, 2m 2 + 3, and 2m 2 + 4, respectively.
Here, we consider the so-called plane stress case. In this setting, the element stiffness matrix where is the Poisson ratio, E 0 ∕(1 − 2 ) is a constant that will be neglected without affecting spectral analysis in the following sections, and the entries are , F I G U R E 1 Local node numbering on the reference element (left) and the global node numbering for the full discretization (right) Appendix A provides more information regarding the stress-strain relation and the assumptions used to arrive at the element stiffness matrix above. Moreover, Appendix A includes a motivation why the Poisson's ratio must be in the range ∈ (−1, 1) in the two-dimensional setting. This is in contrast to the typical bound ∈ (−1, 0.5), which holds for the three-dimensional setting. From an application viewpoint, we are mainly interested in the case where the material is isotropic. Therefore, we primarily consider Poisson's ratios in the range ∈ [0, 0.5).

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 numerical tests, that confirm the theoretical analysis. In particular, Section 3.1 contains the necessary preliminary concepts and tools, while Sections 3.2 and 3.3 focus on the specific study in 2D, where Section 3.2 covers the constant coefficient case and Section 3.3 covers the variable coefficient case.

Premises
The premises include the formal definition of 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.

Singular value/eigenvalue distributions
We first give the formal definitions, and then we briefly discuss the informal and practical meaning.

Definition 1.
Let {A n } n be a sequence of matrices, with A n of size d n , and let f ∶ D ⊂ R t → C r×r be a measurable function defined on a set D with 0 < t (D) < ∞.
• We say that {A n } n has a (asymptotic) singular value distribution described by f , and we write • We say that {A n } n has a (asymptotic) spectral (or eigenvalue) distribution described by f , and we write If {A n } n has both a singular value and an eigenvalue distribution described by f , then we write {A n } n ∼ , f . Notice that (9) is a generalization of (1). More specifically (1) reduces to (9) when the size r of the matrix-valued symbol f is equal to 1. (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 A n . This also holds for the problem studied in this article, as will be further discussed in terms of the stiffness matrices in Sections 3.2 and 3.3 as well as in the block Toeplitz/block diagonal sampling structures in Section 3.1.2.) As already mentioned in the introduction, the symbol f contains spectral/singular value information briefly described informally as follows. With reference to (9), assuming that N is large enough, the eigenvalues of K N are partitioned into r subsets of the same cardinality, except possibly for a small number of outliers, and the ith subset is approximately formed by the samples of i (f ) over a uniform grid in D, i = 1, … , r. It is then clear that the symbol f provides a "compact" and a quite accurate description of the spectrum of the matrices K N for N large enough. Relation (8) has the same meaning when referring to the singular values of K N and by replacing

3.1.2
Multilevel block Toeplitz matrices, multilevel block diagonal sampling matrices, and multilevel block GLT sequences We start by introducing the multi-index notation, which is quite useful for treating sequences of matrices arising from the discretization of PDEs. A multi-index i ∈ Z d , also called a d-index, is simply a (row) vector in Z d ; its components are denoted by i 1 , … , i d .
• 0, 1, 2, … are the vectors of all zeros, all ones, all twos, … (their size will be clear from the context).
lexicographic ordering is assumed uniformly For instance, in the case d = 2 the ordering is the following:

Block Toeplitz matrices.
We now briefly summarize the definition and few relevant properties of multilevel block Toeplitz matrices, that we will employ in the analysis of the stiffness matrices. Given n ∈ N d , a matrix of the form with e vector of all ones, with blocks A k ∈ C r×r , k = −(n − e), … , n − e, is called a multilevel block Toeplitz matrix, or, more precisely, a d-level r-block Toeplitz matrix. Let ∶ [− , ] d → C r×r a matrix-valued function in which each entry belongs to L 1 ([− , ] d ). We denote the Fourier coefficients of the generating function aŝ where the integrals are computed component-wise and (k, ) = k 1 1 + · · · + k d d . For every n ∈ N d , the nth Toeplitz matrix associated with is defined as or, equivalently, as where ⊗ denotes the (Kronecker) tensor product of matrices, while J (l) m is the matrix of order m whose (i, j) entry equals 1 if i − j = l and zero otherwise. We call {T n ( )} n∈N d the family of (multilevel block) Toeplitz matrices associated with , which, in turn, is called the generating function of {T n ( )} n∈N d .
(Multilevel) block diagonal sampling matrices. For n ∈ N and a ∶ [0, 1] → C r×r , we define the block diagonal sampling matrix D n (a) as the diagonal matrix For n ∈ N d and a ∶ [0, 1] d → C r×r , we define the multilevel block diagonal sampling matrix D n (a) as the block diagonal matrix with the lexicographical ordering (10) as discussed at the beginning of Section 3.1.2. Zero-distributed sequences. According to Definition 1, a sequence of matrices {Z n } n such that is referred to as a zero-distributed sequence and we emphasize that this notion applies in the sense of the singular values. Note that, for any r ≥ 1, {Z n } n ∼ 0 is equivalent to {Z n } n ∼ O r (notice that O m and I m 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 article, we use the natural convention 1∕∞ = 0.

Proposition 1.
Let {Z n } n be a sequence of matrices, with Z n of size d n , 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.
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 {A n } n equipped with a measurable function ∶ [0, 1] d × [− , ] d → C r×r , the so-called symbol. We use the notation {A n } n ∼ GLT to indicate that {A n } n is a GLT sequence with symbol . The symbol of a GLT sequence is unique in the sense that if {A n } n ∼ GLT and The main properties of r-block GLT sequences proved in Reference 21 are listed below. If A is a matrix, we denote by A † the Moore-Penrose pseudoinverse of A (recall that A † = A −1 whenever A is invertible).
and only if {Z n } n ∼ 0 (zero-distributed sequences coincide exactly with the GLT sequences having GLT symbol equal to O r a.e. and hence equal to O s a.e. for any positive integer s).

The constant coefficient case ( ≡ 1)
Hereafter we are looking for identifying the symbol underlying the constant coefficient case ≡ 1, according to the standard notion of generating function in the Toeplitz theory (see Section 3.1.2 and Reference 22 for more details) and according to the notion of symbol reported in Definition 1.

Symbol definition
Let A n (1, DN 3 ) be the stiffness matrix obtained with a Q 1 FEs approximation as described in Section 2.2 with proper boundary conditions, Dirichlet in one side "D" and Neumann "N" in the other three (and hence the formal notation A n (1, DN 3 )), where we have chosen a uniform meshing with n intervals in the x 1 direction and n intervals in the x 2 direction. According to the previously considered ordering of the nodes, the matrix A n (1, DN 3 ) is a two-level block tridiagonal structure of size n with tridiagonal blocks of size n + 1, whose elements are small matrices of size 2. We notice that the size is dictated by all the mesh points (including those lying in the boundaries) when considering Neumann boundary conditions, and by all the internal mesh points (excluding those lying in the boundaries), when considering Dirichlet boundary conditions, so that the precise structure is the following: with Â 1 = Â T −1 andÃ 0 slightly differing from Â 0 due to the Neumann boundary condition on the right vertical border, and where Again with the superscript "∽" we are denoting slightly different entries due to Neumann boundary conditions. We notice that the dimension of the blocks A i equals double the number of nodes on the vertical border, while their number equals the number of nodes on the horizontal border minus one due to the Dirichlet boundary condition on the left vertical border. Finally, every a s,t , s, t ∈ {−1, 0, 1} is a 2 × 2 matrix, because two degrees of freedom are associated to each node of the mesh. More precisely, it holds . Now setting n = (n 1 , n 2 ), n 1 = n, n 2 = n + 1 we define Therefore, in accordance with the multi-index notation and with the definition of multilevel block Toeplitz matrices in Section 3.1.2, we can easily read the generating function f Q 1 , just by having in mind that we are dealing with a matrix-valued one of size 2 × 2. More precisely, we have = (a 0,0 + a 0,−1 e −̂1 + a 0,−1 ê1) + (a −1,0 + a −1,−1 e −̂1 + a −1,1 ê1)e −̂2 + (a 1,0 + a 1,−1 e −̂1 + a 1,1 ê1)ê2 where f 11 ( 1 , 2 ) = 4k 1 + 4k 3 cos 1 + 4k 5 cos 2 + 4k 7 cos 1 cos 2 , f 22 ( 1 , 2 ) = 4k 1 + 4k 5 cos 1 + 4k 3 cos 2 + 4k 7 cos 1 cos 2 , Finally, by recalling the expression of k i in (7), we deduce that cos 1 cos 2 , Of course, when the boundary conditions are uniformly of Dirichlet type, then A n (1, D 4 ) = T n (f Q 1 ) with n = (n 1 , n 2 ), n 1 = n 2 = n − 1, since only the internal meshpoints are involved in the matrix. By collecting all the previous results, the following proposition links in a precise way the Toeplitz structures with matrices A n (1, D 4 ) and A n (1, DN 3 ).
Proposition 2. Let f Q 1 ( 1 , 2 ) be the symbol defined in (11). Let us consider a uniform meshing in both directions with n sub-intervals. Then we have the following relationships where R n has rank bounded by 2(n + 1) due to the termÃ 0 − A 0 plus 4(n − 1) due to the terms Â j − A j , j = −1, 0, 1.
If the more general setting is considered with n 1 sub-intervals in the x 1 direction and n 2 sub-intervals in the x 2 direction, then analogous relations are true: where R n has rank bounded by 2(n 2 + 1) due to the termÃ 0 − A 0 plus 4(n 1 − 1) due to the terms Â j − A j , j = −1, 0, 1. In the following Section 3.2.

Numerics: The eigenvalue distribution
We start our analysis by performing a few numerical tests. First of all, in Figure 2 we draw the two eigenvalues surfaces of the symbol f Q 1 with a sampling in (− , ) × (− , ) with respect to different values together with corresponding contour lines. The figure clearly shows that there is a minimum at (0, 0) for both surfaces, regardless of the value of . Moreover, the eigenvalue surfaces for different values of share the same general features. However, there are some noticeable trends that one might find interesting. In particular, for the second surface, both the maximum value and the minimum decrease as increases in (0, 0.5) with the minimum decreasing faster than the maximum. Concerning the maximum values of the two surfaces, Figure 3 highlights that the maximum of the first surface is 4 independently of , while the maximum of the second surface is decreasing from 2 to 1.6 as long as increases in (0, 0.5). The picture is similar if we consider the full interval (−1, 0.5).
We have computed the minimal eigenvalue and the spectral condition number of matrices of increasing dimension with = 0.4, both in the case of Dirichlet boundary conditions and Dirichlet+Neumann boundary conditions (see Table 1). In both cases it is evident that the ratio min (h)∕ min (h∕2) tends to 4 and the ratio 2 (h)∕ 2 (h∕2) to 1∕4 in accordance with the above conjecture. No significant dependency on is present.
Finally, we would like to stress as the match between the samplings of the symbol eigenvalues and the matrix eigenvalues is really sharp even for the moderate size of the involved matrices. To this end in Figure 4, we report the ordered union of equispaced samplings of the two surfaces 1 ( ) (red line) side by side to the ordered eigenvalues of the matrix A n (1, D 4 ) for N(n) = 7938 and different values of . No significant dependency on is observed for both types of boundary conditions (see Figure 5).

3.2.3
Symbol spectral analysis: Distribution, extremal eigenvalues, and conditioning This section comprises three propositions that focus on spectral distribution, extremal eigenvalues, and conditioning of our matrix sequences in 2D. All the results are based on the symbol and its analytical features. We start with a result devoted to the spectral distribution, whose numerical evidences have been discussed already in Section 3.2.2.  DN 3 )} n , first observe that the matrix sequence {R n } n defined in Proposition 2 is zero-distributed thanks to Propositions 1 and 2, since rank(R n )∕size(R n ) → 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 analyze few key analytical features of the spectral symbol, which is shared by all the matrix sequences mentioned in Proposition 3. Such a study is important for the analysis of the extremal eigenvalues and of the conditioning of the same matrix sequences and, in Section 5, it will be the main ingredient for designing ad hoc multigrid solvers, when dealing with the related large linear systems.
Theorem 1. Let f Q 1 ( 1 , 2 ) be the symbol defined in (11). The following statements hold true: both eigenvalues of f Q 1 have a zero of order 2 at (0, 0).

F I G U R E 5 Ordered equispaced samplings of
, j = 1, 2 (red line) and ordered eigenvalues l (A n (1, DN 3 )), l = 1, … , N(n), N(n) = 8320 (blue line) and the row-sum vanishes as ∑ 8 i=1 k i = 0 according to (7). Claim 2. It is just a direct check. Indeed we can evaluate the eigenvalues of the symbol as where tr(f Q 1 ( 1 , 2 )) and det(f Q 1 ( 1 , 2 )) denote the trace and the determinant of the matrix-valued symbol respectively. By expanding the previous expression, we obtain that (cos 1 + cos 2 + 2 cos 1 2 − 4) ± ( + 1)| cos 1 cos 2 − 1| ) with > −1. Finally, by considering the Taylor expansion centered at (0, 0) we have whereas the second order derivatives form a positive definite Hessian matrix at (0, 0). ▪ Since the minimal eigenvalue of the symbol f Q 1 ( 1 , 2 ) has a unique zero of order two at ( 1 , 2 ) = (0, 0) (in fact both eigenvalues of the symbol f Q 1 ( 1 , 2 ) have a unique zero of order two at ( 1 , 2 ) = (0, 0)) and since the symbol is positive semi-definite, we can deduce and discuss important information on the extremal eigenvalues and on the conditioning of the corresponding matrix sequences.
Proposition 4. Let f Q 1 ( 1 , 2 ) be the symbol defined in (11). Let us consider a uniform meshing in both directions with n subintervals. Then we have the following relationships with 2 (⋅) denoting the condition number in spectral norm (the one induced by the Euclidean vector norm || ⋅ || 2 ), and a n ∼ b n for a n , b n ≥ 0 indicating that there exist positive constants c, C independent of n such that c ⋅ a n ≤ b n ≤ C ⋅ a n for any index n (and analogously in the case of multi-indices).
If the more general setting is considered with n 1 subintervals in the x 1 direction and n 2 subintervals in the x 2 direction, then analogous relations are true Proof. The proof relies essentially on the fact that the involved operators are linear and positive (see Reference 30 for a general treatment of the subject in a matrix theoretic context).
By Theorem 1, second item, it follows that f Q 1 ( 1 , 2 ) can be factored as where l is the generating function of the standard 2D discrete Laplacian by centered finite differences of order two and minimal bandwidth, and where g Q 1 has eigenvalues given by which are bounded away from zero and infinity since j (f Q 1 ), j = 1, 2, l( 1 , 2 ) are all nonnegative and with a unique zero of order 2 at (0, 0). More precisely there exist Therefore we infer in the sense of the partial ordering in the (real) vector space of the Hermitian matrices. Now, since T n (⋅) is a linear positive operator, it is also monotone (see e.g., References 23,24), and hence where min (T n (l( 1 , 2 )I 2 )) = 4 − 2 cos By combining (12) and (13), we deduce min (A n (1, D 4 )) ∼ [min n j ] −2 , which reduces to if n 1 = n 2 = n − 1. The other claim that is follows in the same way by duality, simply because Therefore the statements on the conditioning represent a direct consequence of the estimates on the extreme of the spectrum, given the positive definite character of the involved matrices. ▪ Remark 1. For the case of A n (1, DN 3 ) the numerical experiments reported in Table 1 confirms that the extremal eigenvalue behavior and the conditioning should be the same as in the case of the pure Dirichlet boundary conditions. Hence the following relations should hold: The proof could be conducted using the Sherman-Morrison-Woodbury formula since A n (1, DN 3 ) = T n (f Q 1 ) + R n , n = (n 1 , n 2 ), n 1 = n, n 2 = n + 1, with R n of relatively small rank with respect to the matrix size, as in Proposition 2. The idea relies on a clever writing of R n as XY , with X, Y T rectangular matrices of the same sizes, the smaller dimension coinciding with the rank of R n . The same splitting should be computed also in the general case of A n (1, DN 3 ) = T n (f Q 1 ) + R n , n = (n 1 , n 2 + 1). We believe that this direction is worth to be explored in a future investigation.

The non-constant coefficient case
The natural extension of the previous analysis refers to the case of a non-constant coefficient . In the physical model, we are assuming the function piecewise constant with respect to the mesh elements ∈  . Thus, according to a standard assembling procedure, we can write the stiffness matrix as where A El n, is the elementary matrix K n in (6) (possibly properly cut when nodes on the boundary are involved), but widened to size N(n) following the chosen ordering of nodes. Here n = (n − 1, n − 1) in the case of Dirichlet boundary conditions and n subintervals in all the directions so that A n ( ) = A n ( , D 4 ), while n = (n, n + 1) in the case of Dirichlet boundary conditions on one side, Neumann boundary conditions in the other three with n subintervals in all the directions so that A n ( ) = A n ( , DN 3 ).
Clearly the elementary matrix K n is positive semi-definite. More precisely according to values k i in (7), the non-zero eigenvalues are In addition, every elementary matrix, which has been cut due to nodes on the boundary, is positive semi-definite as well being a principal submatrix of K n in (6). Thus, based on (14), we infer min x T A n (1)x ≤ x T A n ( )x ≤ max x T A n (1)x, for all x ≠ 0, with min and max minimum and maximum of , respectively, so that by applying the Courant-Fischer theorem we claim that min min (A n (1)) ≤ min (A n ( )) ≤ max min (A n (1)), min max (A n (1)) ≤ max (A n ( )) ≤ max max (A n (1)).
Now putting together the inequalities in (15) and Proposition 4, we deduce that the extremal eigenvalues and the conditioning have the same asymptotical behavior, as in the constant coefficient setting with Dirichlet boundary conditions. When Neumann boundary conditions are used, the result would be true if the conjecture reported in Remark 1 will be proven.
A deeper analysis of the whole eigenvalues distribution can be conveniently performed by referring to the quite general and powerful GLT theory, 22 whose basics have been introduced in Section 3.1.2. Let D n ( ) be a multilevel block diagonal sampling matrix according to the notions introduced in Section 3.1.2 and let A n (1) be the multilevel block Toeplitz matrix T n (f Q 1 ) if Dirichlet boundary conditions are used or T n (f Q 1 ) + R n in the other case. Then the following facts hold: Fact 1 {D n ( )} n ∼ GLT according to Item GLT 2. Fact 2 {R n } n ∼ GLT 0 according to Proposition 2, Proposition 1, and Item GLT 2.
Fact 4 {A n (1)} n ∼ GLT f Q 1 according to Fact 2, Fact 3, and to the * -algebra structure of GLT sequences that is Item GLT 3. Fact 5 given Δ n = A n ( ) − D n ( )A n (1) a simple check shows that {Δ n } n ∼ GLT 0 since it is obviously true that {Δ n } n ∼ 0, by invoking Proposition 1. Fact 6 {A n ( )} n ∼ GLT f Q 1 as a consequence of Fact 5, Fact 1, Fact 4, and of the * -algebra structure of GLT sequences that is Item GLT 3.; moreover since the matrix sequence {A n ( )} n is made up by Hermitian matrices, by Item GLT 1. it follows that {A n ( )} n ∼ , f Q 1 . Fact 7 {A −1 n (1)A n ( )} n ∼ GLT as a direct consequence of Fact 4 and Fact 6, taking into account * -algebra structure of GLT sequences that is Item GLT 3. Here, for the present preconditioned matrix sequences, we can conclude that the eigenvalues are distributed as , even if the involved matrices are not Hermitian. The reasoning is as follows: (1)} n ∼ GLT since the square root of a positive definite GLT matrix sequence is still a positive definite GLT matrix sequence (see Reference 19). Therefore and by similarity, we infer that In the preceding series of statements we have used the notation {A n ( )} n for several matrix sequences. More in detail, we stress that the main facts that is Fact 6 and Fact 7 apply to all the matrix sequences In fact, the previous matrix sequences represent the variable coefficient versions of the different matrix sequences reported in Proposition 3, for the constant coefficient case when ≡ 1.

TWO-GRID AND MULTIGRID METHODS
In this section, we concisely report few relevant results concerning the convergence theory of algebraic multigrid methods 31 , with special attention to the case of multilevel block Toeplitz structures generated by a matrix-valued symbol f . We start by taking into consideration the generic linear system A m x m = b m with large dimension m, where A m ∈ C m×m is a Hermitian positive definite matrix and x m , b m ∈ C m . Let m 0 = m > m 1 > · · · > m s > · · · > m s min and let P s+1 s ∈ C m s+1 ×m s be a full-rank matrix for any s. At last, let us denote by  s a class of stationary iterative methods for given linear systems of dimension m s .
In accordance with Reference 32, the algebraic Two-Grid Method (TGM) can be easily seen as a stationary iterative method whose generic steps are reported below, where we refer to the dimension m s by means of its subscript s. More V pre s,pre , s = s min − 1, … , 0. Remark 2. In the relevant literature (see, for instance, Reference 33), the convergence analysis of the TGM splits into the validation of two separate conditions: the smoothing property and the approximation property. Regarding the latter, with reference to scalar structured matrices, 33,34 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 T n (f ) with n = (2 t − 1), f a nonnegative trigonometric polynomial. Let 0 be the unique zero of f . Then the TGM optimality applied to T n (f ) is guaranteed if we choose the symbol p of the family of projection operators such that lim sup where, for d = 1, the sets Ω( ) and ( ) are the following corner and mirror points 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| 2 ( ) + |p| 2 ( + ) does not have zeros and |p| 2 ( + )∕f ( ) is bounded (if we require the optimality of the V-cycle then the second condition is a bit stronger); see Reference 33. When discretizing differential operators, the previous conditions mean that p has a zero of order at least at = , whenever f has a zero at 0 = 0 of order 2 , since the involved Toeplitz-like structures have a spectral symbol with a unique zero at 0 = 0 of order 2 .
In our specific block setting, by interpreting the analysis given in Reference 35, all the involved symbols are matrix-valued and the conditions which generalize (16) 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, W-cycle procedures converging optimally or quasi-optimally with respect to all the relevant parameters (size, dimensionality, 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 in the conclusion section in Reference 35.

MULTIGRID PROPOSALS
According to the theory in Section 4, we define the prolongation operator as follows in the case of matrices A n (1, D 4 ) and in the case of matrices A n (1, DN 3 ). In the expressions above, Note that A 2h = P T A h P equals the same FEM approximation on the coarse mesh in both cases, independently of the boundary conditions. Again, independently of the boundary conditions, the symbol of the prolongation operator is p( 1 , 2 ) = (1 + cos( 1 ))(1 + cos( 2 ))I 2 .
Since the matrix-valued function p( 1 , 2 ) is a scalar function times the identity, it is evident that the quite delicate condition C) is automatically satisfied. Furthermore, taking into account Theorem 1, since the scalar function 1 + cos( ) has a zero of order 2 at = , we deduce that condition A) holds, while the positive definiteness of ∑ ∈Ω( ) pp H ( ) is also trivially verified so that also condition B) is met.
In addition, one iteration of Gauss-Seidel as pre and post smoother is considered. In view of 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. Table 2 confirms such a theoretical forecast.
We have performed numerical experiments also with respect to other solvers. The incomplete Cholesky factorization preconditioning and the block circulant preconditioning with proper Strang-type correction (see Table 3), as expected, both improve the convergence with respect to the classical conjugate gradient. However, the iteration count grows (according to the theory) as the square root of the matrix size of the involved matrices. On the other hand, our multigrid strategy is not only optimal (number of iterations independent of the matrix size), but also the concrete iteration count is very small (less than 10 iterations) and the computational cost is optimal that is linear and proportional to the cost of the matrix-vector product.
In the last numerical tests reported in Tables 4 and 5, we compare our approach both as a multigrid solver and a one-step iteration in a preconditioned flexible conjugate gradient (FCG) to professional software (see References 36-40 and references therein). As we can see from the subsequent tables, our procedures employ fewer iterations, with a comparable cost per iteration, and show clear optimality. The optimality is not observed using the Notay codes. For the sake of completeness, we again remark that the CPU timings are not reported since our codes are really basic and indeed they would need further optimization for a fair comparison with professional and well-established codes. At any rate, even if the CPU times are not reported, we are sure that also the cost of a single iteration is favorable in our approach both in terms of theoretical complexity and in terms of timing. The reason is that we do not manipulate matrices when using our prolongation/restriction operators but only polynomial symbols: whence only O(1) coefficients are used, and the related timings are very compressed.

CONCLUSIONS
In this work, we have considered a problem that stems from topology optimization, which aims to find the best material layout subject to assigned constraints. When solving the related governing equation using the FEM, a large number of elements is employed to discretize the design domain, and an element-wise constant function approximates the coefficient field in the considered 2D design domain. First, we have provided a spectral analysis of the coefficient matrices associated with the linear systems stemming from the FE discretization of a linearly elastic problem, for an arbitrary coefficient field. Based on the spectral information, we have 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 matrix-vector cost. The method has been tested, and the numerical results are very satisfactory, in terms of linear cost number of iterations, which is bounded by a constant independent of the matrix size and lightly influenced by the desired accuracy. We finally remark that the used tools are very flexible, and in particular, they do not depend on the number of spatial dimensions of the underlying problem. Therefore, we are confident that a generalization of our spectral analysis and the related algorithmic proposals can be obtained in a three-dimensional setting as well, even if to the price of quite heavy notations. Finally, what is conjectured in Remark 1 and it is numerically observed in Table 1 will be the subject of future investigations, in order to complete the spectral analysis of the considered matrices, with Neumann boundary conditions on at least one side.

ACKNOWLEDGMENT
This work is financially supported by the Swedish strategic research programme eSSENCE, a strategic collaborative eScience program funded by the Swedish Research Council, and the Italian Institution for High Mathematics (INdAM).

CONFLICT OF INTEREST
This study does not have any conflicts to disclose.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. Similarly, corresponding to the six strain components above, there are also six independent stress components written in vector form as By the generalized Hooke's law, the most general linear relation among components of the stress and strain tensor can then be written as

ORCID
where E is a matrix that corresponds to the constant fourth-order elasticity tensor. The relationship between stresses and strains is 11 = 1 E 0 ( 11 − ( 22 + 33 )) , 12 where is Poisson's ratio, E 0 is Young's modulus, and is the shear modulus. However, in this article, we consider the two-dimensions plane stress condition that implies Moreover, by combining relations for 11 and 22 , we find that 33 = − 1− ( 11 + 22 ). Hence, the stress-strain relation can be reduced to the stress and strain components in the x 1 x 2 -plane.
Analogous the three-dimensional case, equation system (A4) can be inverted to obtain Hooke's law (A1) with In this case, the bulk modulus K can be expressed as Also in this case, the Young (E 0 ), shear (G), and bulk (K) moduli need to be positive. Thus, Equations (A3) and (A5) imply that the Poisson ratio in the two dimensional plane stress setting must satisfy −1 < < 1. We remark that the upper-bound in this case is larger than the upper-bound in the three dimensional setting.
Remark 3. Doing a similar analysis as above for the two-dimensional plane strain as well as in the three-dimensional case yields that the Poisson's ratio belongs to an open interval of the form −1 < < 0.5. As a consequence, the plane stress setting allows for exotic (unphysical) materials with Poisson's ratios larger than 0.5, which is the incompressibility limit for the underlying physical problem (the three-dimensional case). To the best of our knowledge, there are no known natural occurring isotropic materials for < 0. On the other hand, one can design isotropic material with negative Poisson's ratio using theoretical analysis and simulations. 41 Nevertheless as a matter of practice, 42 we limit our attention to Poisson's ratio satisfying the following inequalities 0 ≤ < 0.5. However, in some more delicate circumstances, the Poisson ratio should lie in the interval 0.2 ≤ < 0.5, as proved experimentally in view of elastic properties of real isotropic materials. 43