Fast multiscale contrast independent preconditioners for linear elastic topology optimization problems

06/23/2020 ∙ by Miguel Zambrano, et al. ∙ Lawrence Livermore National Laboratory Universidad Nacional de Colombia 0

The goal of this work is to present a fast and viable approach for the numerical solution of the high-contrast state problems arising in topology optimization. The optimization process is iterative, and the gradients are obtained by an adjoint analysis, which requires the numerical solution of large high-contrast linear elastic problems with features spanning several length scales. The size of the discretized problems forces the utilization of iterative linear solvers with solution time dependant on the quality of the preconditioner. The lack of clear separation between the scales, as well as the high-contrast, imposes severe challenges on the standard preconditioning techniques. Thus, here we propose new methods for the high-contrast elasticity equation with performance independent of the high-contrast and the multi-scale structure of the elasticity problem. The solvers are based on two-levels domain decomposition techniques with a carefully constructed coarse level to deal with the high-contrast and multi-scale nature of the problem. The construction utilizes spectral equivalence between scalar diffusion and each displacement block of the elasticity problems and, in contrast to previous solutions proposed in the literature, is able to select the appropriate dimension of the coarse space automatically. The new methods inherit the advantages of domain decomposition techniques, such as easy parallelization and scalability. The presented numerical experiments demonstrate the excellent performance of the proposed methods.



There are no comments yet.


page 6

page 9

page 14

page 18

page 19

page 20

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Topology optimization Bendsøe and Sigmund (2003) is an iterative design process aiming to find close to optimal material distribution by minimizing an objective function and fulfilling a set of constraints. More precisely, the discrete optimization problem, considered here, can be written as


where is the objective function equal to the compliance of the mechanical system, represents the associated physical problem or the state equations written in residual form, and

is an additional constraint on the volume of the solid material distributed in the design domain. The vectors

and represents the state solution and the material distribution respectively. As the focus here is on linear elastic problems discretized using the finite element method, the state equations can be written as


where is the vector with the external forces applied on the system and is the so-called stiffness matrix obtained using standard finite element assembly. The computation domain is discretized using finite elements and a density value, bounded between zero and one, is assigned to each of them. Void elements are modeled by assigning the density to zero, and parts occupied with solid material are modeled with density values equal to one. All density variables, or also called design variables, are collected in the density vector . The density values are allowed to vary continuously between zero and one in order to utilize gradient optimization techniques for finding a material distribution fulfilling the constraints and minimizing the objective. The vector collects all areas/volumes of the discrete finite elements.

The actual physical material distribution is calculated by a set of transformations applied on the original density field , e.g., Lazarov et al. (2016). The first transformation is usually obtained by convoluting the density distribution with a filter functionsBendsøe and Sigmund (2003) resulting in the so-called filtered density field and providing a mesh independent solution of the optimization problem. The filtered density can be utilized directly for modeling the stiffness by using the SIMP Bendsøe (1989)

(Solid Isotropic Material with Penalisation) material interpolation scheme where the modulus of elasticity for every element is computed as


where is the modulus of elasticity of the solid material and is a small regularization parameter ensuring the existence of a solution to the associated linear elastic problem, and is a penalization parameter often taken to be equal to . Alternatively, as the filtered field consists of many elements with densities between zero and one, additional projection step Lazarov et al. (2016) can provide a sharper transition between void and solid. Here, the physical density is modeled directly by the filtered field, however, the presented preconditioning techniques can be applied to formulations with projections, penalization techniques different than SIMP Bendsøe and Sigmund (2003), to the so-called robust formulation in topology optimization Wang et al. (2011), and other problems with high-contrast and multi-scale coefficients like level-set type of topology optimization formulations, simple parametric studies, and simulations van Dijk et al. (2013).

As stated earlier, the solution of the topology optimization problem is obtained iteratively. The optimization process starts with some admissible initial design . The state solution is computed by solving Equation 5 and the gradients of the objective are evaluated by solving an adjoint equation Bendsøe and Sigmund (2003) in the general case. For the minimum compliance problem, the gradients of the objective are given as


where refers to the element index. The densities are usually updated using the Method of Moving Asymptotes (MMA) Svanberg (1987) or the optimality criteria method (OC) Bendsøe and Sigmund (2003). For a more detailed introduction to topology optimization, the interested readers are referred to Bendsøe and Sigmund (2003) and Andreassen et al. (2011).

1.1 Physical problem and its discretization

The systems considered here (represented above as and simplified to

) are linear elastic and their response is obtained by solving the Navier-Cauchy partial differential equation



is the stress tensor,

is the strain tensor given by


and is an elastic material properties tensor, denotes the displacement field and is the input supplied to the system, i.e., distributed and concentrated forces. The mechanical system occupies the bounded physical domain . The boundary , is decomposed into two disjoint subsets for each component , with prescribed displacements , and with prescribed traction .

The elastic material properties tensor is isotropic, and for a point in the computational domain is computed as


In the above equation, is a constant tensor obtained for predefined Poisson ratio and modulus of elasticity one. For a point located in an element the elastic modulus is obtained using Equation 6. We refer to Buck et al. (2013, 2014) for design for domain decomposition preconditioners for the elasticity equation in the general case.

The variational formulation Braess (2007) of Equation 8 reads


with bilinear form and linear form


where . The weak formulation is discretized using finite element space with vector valued shape function defined on a uniform mesh . Each basis function is a scalar bilinear Lagrange function in one of the components and zero in the other. Substituting the shape functions in the integrals given by Equation 12 for all finite elements in the mesh leads to the linear system of equations previously introduced as Equation 5.

1.2 Iterative solvers in topology optimization

Topology optimization is a computationally heavy process. The resolution of the obtained designs depends both on the transformation of the design field and the discretization of the design domain. Fine discretization is capable of representing small design features which leads to better utilization of the design freedom and at the same time to a larger system on linear equations. The total computational time is usually proportional to the number of design iterations. Every update of the design requires a negligible amount of time compared to the time necessary for solving the state problem Aage and Lazarov (2013). The solution of the linear system Equation 5 dominates the computational cost and requires careful selection of a solution algorithm and scalable implementation for large 3D problems Aage et al. (2017).

Factorization techniques are serial by nature and hard to parallelize. On the other hand, iterative linear solvers Saad (2003)

are relatively easy to parallelize and provide a scalable alternative to direct factorization techniques. The convergence rate depends on the condition number of the stiffness matrix and the clustering of the eigenvalues. The introduction of weak background material, with a modulus of elasticity orders of magnitude softer than the distributed solid, leads to an ill-conditioned system of linear equations. More precisely the condition number is of order

where is a characteristic measure of the mesh size, and measures the contrast in the coefficients . Furthermore, the optimal design often consists of multiscale segments, see Figure 1, which together with the bad condition number makes the solution of the linear system extremely challenging Galvis and Efendiev (2010b, a); Efendiev et al. (2012).

Preconditioning techniques alleviate the slow convergence speed. Therefore, the art of solving large sparse linear systems in parallel lies in the construction of computationally cheap, parallelizable and effective preconditioners. A preconditioned system is obtained from Equation 5 by premultiplication with a preconditioner . The most effective preconditioner is the exact inverse of the stiffness matrix . However, direct construction requires matrix factorization which as already discussed is not scalable and is extremely expensive for large problems. Thus, a preconditioner in the form of a multigrid procedure Vassilevski (2008) or a domain decomposition Toselli and Widlund (2004) provides the most efficient solution procedure.

1.3 High-contrast coefficients in topology optimization

The focus here is on domain decomposition preconditioners, in particular, on two-levels domain decomposition techniques. The classical variants of these preconditioners Toselli and Widlund (2004) do not perform well for high-contrast problems Efendiev et al. (2012)

. The condition number estimates for the traditional domain decomposition case depends on the contrast

if the high-conductivity regions are not aligned with the coarse mesh of the decomposition. We refer to Efendiev et al. (2012, 2013a); Buck et al. (2013, 2014) where it is demonstrated that if the material properties with local variations are enclosed in a coarse block the performance of the standard preconditioner is not affected by the contrast. However, for cases with several extended high-conductivity areas (high stiffness regions for linear elasticity) crossing the coarse block boundaries, the performance deteriorates significantly. That is precisely the case for topology optimized linear elastic structures. These design features can be observed in Figure 1 as well as in previous articles on the topic Lazarov (2014); Alexandersen and Lazarov (2015b, a).

We apply the Generalized Multiscale Finite Element methods (GMsFEM) framework introduced in Efendiev and Galvis (2011); Galvis and Efendiev (2010a); Efendiev et al. (2013a) to construct robust and fast solution algorithms adapted to linear elasticity. Similar to the constructions of diffusion type multilevel domain decomposition preconditioners, the proposed design depends on the behavior of the coefficient inside local coarse node neighborhoods. We demonstrate that the preconditioned solvers perform, in terms of iterations and condition number estimates, independently of the contrast in the media properties.

The most important part of the construction of a multilevel domain decomposition preconditioner is well know to be the coarse level. The coarse level of the preconditioner has to provide a good local approximation of the kernel of the elasticity operator Efendiev and Galvis (2010). Also, it should contain all eigenmodes with corresponding eigenvalues dependant on the contrast of the coefficient. Thus, an adequately chosen eigenvalue problem is constructed and solved locally to ensure the desired behavior.

Two major approaches can be identified in the current literature. For coarse spaces with standard dimension (for scalar problems that is, one basis function per coarse node and diffusion coefficient ) it is imperative to have coarse basis functions such that is bounded independently of the contrast. For high-conductivity regions (high stiffness regions for linear elasticity) restricted within the coarse element, i.e., there are not any long channels crossing the edges of the coarse element, the above requirement is fulfilled by classical multiscale finite element basis functions obtained by energy minimization Buck et al. (2013). For high-contrast coefficients with extended channels, even for isotropic problems, the above condition cannot be fulfilled. Therefore, to achieve robust behavior, an enrichment procedure is implemented by adding basis functions that approximate the contrast dependent eigenmodes of the operator locally.

Apart from the fact that the material coefficients in topology optimization show multiscale variations combined with high-contrast, an additional complication comes from the fact that throughout the optimization iteration the density field and correspondingly are changing as the optimization converges to the optimized design. A topology represented by a particular density distribution evolves with the iterations leading to iteration dependent multiscale structure. High-contrast channels may break apart or joint together during the optimization iterations within a globally connected subdomain restricted to a coarse neighborhood. Dealing with such an additional complication requires re-computation of the preconditioner as the optimization iteration advances towards the final solution. Similar to Lazarov (2014); Alexandersen and Lazarov (2015a), we pay extra attention to the building cost of the preconditioner, especially to the construction of the coarse level. Therefore, two new alternatives are proposed to reduce the cost of recomputing the basis functions for the coarse space:

  • Computation of eigenvalue problems using a randomized algorithm. The randomized approximation of the local eigenvalue problems for GMsFEM is proposed in Calo et al. (2016); Efendiev et al. (2013a) and realizes a significant reduction of the cost of computing the coarse basis functions that generate the coarse space. A detailed introduction can be found in the overview article Halko et al. (2011).

  • The utilization of a preconditioner constructed for the heat equation in a similar high-contrast multiscale media to precondition the elasticity equation. Here, we follow the ideas presented in Gustafsson and Lindskog (1998, 2002) for homogeneous media. More precisely, we utilize the coarse spaces generated for the heat equation with a diffusion coefficient defined using the solid material region . The local eigenvalue problem is related to the diffusion operator and not to the elasticity equation operator. Therefore, the local eigenvalue problem, for the same resolution, is twice smaller in two-dimensions and three times smaller in three dimensions reducing the computational cost significantly. The combination with randomized algorithm results in a solver for the topology optimization an order of magnitude faster than the one presented earlier in Alexandersen and Lazarov (2015a).

Remark 1.

An unexpected advantage of using the heat equation to precondition the elasticity equation in the context of GMsFEM framework is that in the case of the heat operator it is easier to define threshold strategies for the decay of the local eigenvalues, and therefore adaptivity to coefficients is easier to obtain.

1.4 Topology optimization example

To demonstrate the performance of the preconditioner we consider a square 2D plane stress example with homogeneous Dirichlet boundary condition and distributed force over the domain, as shown in Figure 1. The setup even though difficult to find in the engineering practice is an excellent example to test preconditioners as the optimization results in a highly sophisticated multiscale design. The design domain is partitioned on square mesh elements with

coarse elements. The maximum number of eigenvectors per coarse neighborhood is taken to be

. The topology optimization algorithm was set to stop at iterations.

As the focus here is on the effective preconditioning techniques, we limit the topology optimization problems to the one presented in Figure 1. The presented topology is obtained using the standard and the preconditioned solvers presented here. Depending on the preconditioner update strategy, the number of the iterations reduced an order of magnitude. The fastest total solution time is obtained by employing the strategy presented in Alexandersen and Lazarov (2015a), i.e., either updating the preconditioner after a fixed number of optimization iterations or updating it after the number of the iterations of the local iterative solver exceeds a predefined threshold. More extensive 3D topology optimization studies will be presented in the following articles, and here we will focus on demonstrating the numerical properties of the developed preconditioners and the associated computational complexity.

The proposed replacement of a standard elasticity with a diffusion solver for the local eigenvalue problem results in a time reduction factor of . The above result emphasizes the contribution of the local solves to the total solution time and is in line with the theoretical predictions. The diffusion problem is three times smaller compared to the elasticity, and the cost of the local eigenvalue solves is proportional to , where

is the number of local degrees of freedom. The randomized solver reduces further the computational cost resulting in an additional factor of two to three. It should be pointed out that the problems considered here are 2D, and the randomized modification is expected to have a more significant effect on the preconditioner time in larger 3D problems.

Figure 1: Design domain and topology optimization of a 2D plane stress problem with homogeneous Dirichlet boundary conditions and distributed force over the domain.

2 GMsFEM two-levels domain decomposition for the elasticity equation

The focus In this section is on the utilization of GMsFEM coarse spaces in the construction of two-levels domain decomposition preconditioners. In particular, we show that the proposed preconditioners yield a contrast-independent condition number, and thus they are optimal in terms of physical parameters. Additional theoretical details can be found in Galvis and Efendiev (2010b); Efendiev et al. (2011b) and further extensions of the results to multilevel methods are presented in Efendiev et al. (2011a).

The computational domain, shown in Figure 2, is partitioned into finite elements which are agglomerated into larger non-overlapping blocks . Based on the above decompotions an overlapping decomposition, shown in Figure 3, denoted as is obtained from an original non-overlapping decomposition by enlarging each subdomain to


where is some distance function. The overlapping subdomains and the coarse triangulation are not related in general settings. Two partitions of unity covering the whole computational domain, one for and one for , can be chosen independently of each other. However, in the numerical experiments presented later in the paper, we assume that the overlapping subdomains coincide with the coarse vertex neighborhoods of , and in this case , where is the overlapping parameter.

Based on the above decomposition, the preconditioned operator is , and the preconditioner matrix is defined as


The part corresponding to the first level is


where , , and the part corresponding to the second (or coarse) level is

where . The matrix consists of vectors defining the so-called coarse space which is obtained by interpolating the coarse functions onto the fine mesh. The matrices are rectangular and consist of zeros and ones and utilized to extract the degrees of freedom that lie inside the subdomains .

Figure 2: Decomposition of the computational domain into non-overlapping coarse blocks (subdomains) .

The fine-scale linear system Equation 5 is solved iteratively with the preconditioned conjugate gradient (PCG) method. The application of the preconditioner involves solving a coarse-scale system and solving a set of local problems in each iteration. The main goal is to reduce the number of iterations in the solution process. Without the coarse space , the preconditioner usually acts as a smoother and the convergence depends on the number of the coarse blocks. Thus, the appropriate construction of the coarse space plays a crucial role in obtaining robust iterative domain decomposition methods and ensures iteration number independent on the contrast and the characteristic mesh size. See Toselli and Widlund (2004); Efendiev et al. (2012, 2013a)

Figure 3: Definition of overlapping domain decomposition and a neighborhood of a coarse node .

2.1 Generalized multiscale coarse finite element spaces for the elasticity equation

The construction of the coarse space for linear elasticity follows the strategy outlined in Efendiev et al. (2011b, 2013b) for the diffusion equation. The process starts with the selection of an initial set of basis functions that form a partition of unity and are associated with the coarse domains . Additional sets of coarse basis functions are defined with respect to the fine mesh . These are computed by solving a local eigenvalue problem


with homogeneous Neumann boundary condition on and Dirichlet boundary condition on if is not empty. In matrix form we have


with eigenvalues and the eigenvectors denoted as and respectively. The eigenvalues are ordered as


The multiscale coarse basis functions are constructed by multiplying selected eigenfunctions

with the partition of unity function associated with subdomain . The coarse GMsFEM space for the entire problem is defined by


where and denote the number of eigenvectors and the number of coarse blocks respectively.

For the diffusion case, the contrast dependent eigenvalues are well separated, i.e., a clear jump can be observed between the contrast dependent and the rest of the eigenmodes computed on a given patch . Thus, selecting the low order contrast dependent modes results in a preconditioner which provides a condition number for the preconditioned operator independent of the contrast. However, for the linear elastic case, such behavior for the low order modes is difficult to be observed. The optimal number of low order modes is related to the disconnected high-stiffness regions and the RBM (rigid body motion) of the region Efendiev et al. (2012). We illustrate this fact in Figure 4 where we picture a coarse node neighborhood with two disconnected high-stiffness regions and its contrast dependent modes. There are three modes for each high-stiffness region, and these correspond to the rigid body modes. The next mode, the one corresponding to the eigenvalue number seven in increasing order, is contrast independent as observed in our numerical tests. Note that, restricted to high-contrast regions, the first 6 modes correspond to three linearly independent rigid body modes per inclusion. The space RBM of rigid body modes on a set is defined for , by

Remark 2.

In most of our numerical experiments using the eigenvalue problems in Equation 17 we did not find any automatic way to implement a threshold to select the adequate number of modes in each coarse neighborhood. The numerical experiments were performed by specifying the number of basis functions based on the number of disconnected high-stiffness regions present in the specific coarse node neighborhood. Later we will introduce the use of GMsFEM basis functions constructed for the heat equations where an automatic threshold can be implemented to select the adequate number of basis functions in each coarse neighborhood.

Figure 4: Two disconnected high-stiffness regions and their contrast dependent modes. Three for each region corresponding to the RBMs. The next mode is contrast independent as observed in our numerical tests. Note that, restricted to high-contrast regions, these contrast dependent modes correspond to linearly independent rigid body motions.
Figure 5: Basis function construction for an element of the domain. To illustrate the process only one displacement component, along the x-direction, of the eigenvector is shown in the coarse region. The values for the other components are obtained following the same procedure.

For the elasticity problem the GMsFEM approximate the solution on the coarse grid as,

where are the unknown constants. The coarse matrix is constructed,


and the multiscale finite element solution is the finite element projection of the fine-scale solution on , that is, where

and .

3 A randomized algorithm for eigenvalues computation

Inspired by Halko et al. (2011), randomized algorithms have been introduced in GMsFEM in Efendiev et al. (2013b); Calo et al. (2016). See also Galvis et al. (2017). The basic idea is to utilize random sampling to generate low-rank approximations to the set of matrices utilized for finding the coarse shape functions Equation 17. The main difference between deterministic factorization techniques is the convergence criteria. For randomized algorithms, the converge is probabilistic. However, as stated in Halko et al. (2011)

, the probability of failure is often negligible of order

. Therefore, for practical applications involving singular values or eigenvalues decompositions, a randomized algorithm offer a computationally cheap alternative to direct factorization methods.

The algorithm for finding approximations of the eigenvalues and the eigenvectors starts with dimension reduction. For each subdomain the following sequence of steps is executed in parallel

  • Generate forcing terms using (for instance) an uniform random distribution and .

  • Compute local solutions

  • Generate

  • Solve a reduced eigenvalue problem

The reduced size matrices are generated as:


where every column of is a vector from . The approximations of the eigenvalues and the eigenvectors are computed as


Usually, the matrix

holds several vectors obtained by singular value decomposition (SVD) of the snapshot matrix

enriched with the three rigid body modes in 2D and six in 3D. The actual number of vectors depends on the desired number of shape functions. As discussed in Halko et al. (2011), see also Calo et al. (2016), for a target number of shape functions , the rank of can be selected to be as low as . The computational time of a standard generalized eigenvalue algorithm is proportional to , where is the size of the problem. Thus, the reduced problem offers significant speed up. However, the cost associated with the solution of the linear system cannot be reduced further. A possible speed up based on splitting of the displacement field is discussed in the following section.

4 Displacement field splitting preconditioner

The basic idea behind the proposed preconditioner for linear elastic isotropic problems is presented in Axelsson and Gustafsson (1978); Blaheta (1994); Gustafsson and Lindskog (1998, 2002). The displacement vector is split in two blocks (three blocks in 3D) using the displacement components aligned with x, and y-directions, i.e., . Following the above decomposition, the stiffness matrix and the load vector , from Equation 5, can be written in a block form as


The displacement split preconditioner matrix is constructed by keeping only the diagonal block matrices.


As demonstrated in Gustafsson and Lindskog (1998, 2002), the condition number of the preconditioned operator depends only on the Poisson’s ratio and is given as


where . Furthermore, the condition number does not depend on the contrast of the materials distributed in the computational domain. The latter makes the block diagonal preconditioner perfect for topology optimization problems. Close inspection of the diagonal block matrices reveals that they are equal for isotropic elastic problems , and are equivalent to a stiffness matrix obtained by the discretization of a scalar Laplace problem. The last property is utilized in the proposed preconditioners to reduce further the computational cost. Instead of constructing the coarse bases by solving the linear elastic eigenvalue problem given by Equation 17, the idea here is to solve a reduced eigenvalue problem associated with the Laplace problem.

The MsFEM preconditioner and the procedure for finding the coarse bases can be applied directly to the diagonal blocks, can be combined with randomized algorithms, or can be utilized to form a coarse space and subsequently project the full elastic matrix. The above options lead to different preconditioners discussed in detail in the following and numerically tested in section 5.

4.1 Local problems and coarse basis from a diffusion operator

The first preconditioner derived from the field splitting preconditioner consists of a fine level block-diagonal preconditioner constructed for the diffusion equation and a coarse-level part obtained from a coarse set of basis functions obtained again using an eigenvalue problem for the diffusion case. The preconditioner can be written as




and is the stiffness matrix obtained by a finite element discretization of the Navier-Stokes equations for linear elasticity. The index denotes quantities obtained from the diffusion case and index quantities obtained from the elasticity equations. For instance, the image of the operator (or column space of the matrix) is the coarse space constructed using the GMsFEM procedure but starting with the diffusion operator . In this case the construction of the coarse basis function is done component-wise (x-and y-direction displacements) and each of those uses the local eigenvalue problem


where or some other quantity that captures the high-contrast and multiscale structure of . This eigenvalue problems is posed with homogeneous Neumann boundary condition on and Dirichlet boundary condition on if is not empty. In this case we have


See Section 2.1 for comparison with the construction of that uses local elasticity eigenproblems. The cost of constructing is less than that of constructing in (19) since we solve smaller local eigenvalue problems. Additionally, we also observed numerically that the selection of contrast-dependent modes can be performed automatically in the case of the local heat operator eigenvalue problem which is harder to do in the case of the local elasticity operator eigenproblem. For the numerical test, we consider the approximation of the eigenvalue problem (29) using a randomized method similar to the one described in Subsection 3.

The first level of the preconditioner (27) contains the operator that, in a similar manner to that of the construction of the coarse basis, it denotes the additive first level of the precoditioner constructed for the diffusion operator . Applying requires the solution of Dirichlet local diffusion equations and then add their extensions by zero outside the corresponding subdomain.

4.2 Local problems and coarse basis from a diffusion operator with rigid body motions enrichment

The previous preconditioner includes the null space for the Laplace operator, which is capable of representing only the rigid body translations. The null space of the elasticity operator is larger. The coarse space must be able to represent all rigid body modes to ensure convergence independent on the problem size Griebel et al. (2003). Thus, the coarse space is enriched with additional vector in 2D for the rigid body rotations. The construction of the vector elements is based on the vector which is the interpolation on the fine-grid of the vector function representing the rotation. Thus, to the basis functions previously constructed we add the vector obtained for each coarse neighborhood. For 3D problems, the coarse basis should include three additional rigid body rotational modes. The modified preconditioner, in this case, is given as


where is the enriched coarse space and


4.3 Elasticity operator local problems and coarse basis from a diffusion operator

The diagonal blocks associated with every coarse partition are relatively small. Even though they are three times larger compared to the diagonal blocks obtained for the Laplace problem, the computational time does not increase significantly. Thus the diagonal blocks of the first part of the preconditioner can be obtained directly from the elastic operator. The second part can be constructed similar to the previous two cases using a coarse space from the diffusion operator and enriched coarse space with rotations. That is we have,




where we use the definitions in (28) and (32), respectively. The first level of the preconditioner, , is defined in (15). These two cases complete the definitions of the preconditioners utilized in the numerical experiments. See Table 1 for a summary of all the implemented iterations.

Notation Definition Level 1 Eigenproblem Enrichment
Eq. (14) Elasticity Full elasticity (16) None
Eq. (27) Heat (blocks) Full heat (29) None
Eq. (31) Heat (blocks) Full heat (29) Rotations
Eq. (33) Elasticity Full heat (29) None
Eq. (34) Elasticity Full heat (29) Rotations
Similar to Eq. (34) Elasticity Randomized heat (29) Rotations
Similar to Eq. (14) Elasticity Randomized elasticity (16) None
Table 1: Information about the elasticity preconditioners used for numerical tests. All the methods solve an elasticity coarse problem, the differences can be found in the local solvers (Column 3) and in the construction of the coarse space where the changes correspond to the neighborhood eigenvalue problem used and its numerical approximation (Column 4). Some methods need enrichment of the coarse space (Column 5).

5 Numerical experiments and results

In addition to the topology optimized case discussed earlier, in order to demonstrate the contrast independence of the above preconditioners, we utilize high-contrast material distribution as shown in Figure 6. The computational domain is partitioned applying coarse mesh and each coarse-element is further partitioned utilizing fine-mesh. The fine discretization is performed with bilinear polynomials. Zero Dirichlet boundary conditions are applied to all boundaries of the computational domain. The tests are performed with forcing terms shown in Figure 7. Both forces are applied to the solid-material subdomain which is the expected case for topology optimization problems. For the solution of the linear system, we run PCG until the relative norm of the initial residual is reduced by a factor of . We use the Lanczos connection method to estimate the condition number of the preconditioned operator; see (Saad, 2003, Chapter 6).

In Table 1 we present the results with the implemented preconditioners for the elasticity equation. After the construction of the coarse basis functions, the second level of all the implemented methods consists in solving an elasticity coarse problem. Some of the preconditioners differ in the level-one local solvers: they either use elasticity equation local solvers or block diagonal heat equation local solvers. We then have several proposed coarse spaces constructions. The different constructions correspond to the local eigenvalue problem used and its numerical approximation (Column 4). We pose either a local elasticity eigenvalue problem or a local diffusion eigenvalue problem. For the numerical approximation of the eigenvalues and eigenvectors we use either a full eigenvalue solver of the fine scale local eigenvalue problem or the randomized method of Subsection 3. Some coarse spaces need to be enriched (in order to obtain the RBM) by adding rotations multiplied by partition of unity functions as additional coarse basis functions (Column 5).

In Table 2-5 we present iteration count and condition number estimates for the different iterations that have been introduced as they behave with respect to the contrast. We summarized this results in Table 6 and 7. We observe that the methods that are robust with respect to the contrast are , , and ; see Table 1. Therefore, they are computationally efficient alternatives to solve the elasticity equation and can be used in topology optimization problems such as the one considered in this paper. As it was mentioned before, an advantage of is that, according to our numerical experiments and for general multiscale configurations, it allows us to identify the contrast-asymptotically-vanishing eigenvalues and the corresponding eigenvectors.

Figure 6: Material distribution and coarse mesh for the numerical experiments.
Figure 7: Forcing term for the numerical experiments.
Preconditioner Iterations Condition Coarse dim.
None 293
14 4.6 243
using twice many eigenvectors 15 5.0 486
29 21.0 243
29 21.0 243
15 5.2 243
14 4.6 243
with 10 snapshots 24 13.8 243
with 15 snapshots 25 15.6 243
with 10 snapshots 20 9.4 243
with 15 snapshots 20 9.4 243
Table 2: Results for the elasticity problem with contrast . See Table 1 for the description of the preconditioners.
Preconditioner Iterations Condition Coarse dim.
None 1583
26 16.0 387
using twice many eigenvectors 19 7.7 774
86 387
81 387
35 72.0 387
28 18.3 387
with 10 snapshots 32 26.8 387
with 15 snapshots 33 27.0 387
with 10 snapshots 30 22.9 387
with 15 snapshots 30 22.8 387
Table 3: Results for the elasticity problem with contrast . See Table 1 for the description of the preconditioners.
Preconditioner Iterations Condition Coarse dim.
None >2000
53 113.8 387
using twice many eigenvectors 28 23.6 774
200 387
117 387
69 387
56 108.9 387
with 10 snapshots 62 111.1 387
with 15 snapshots 62 111.4 387
with 10 snapshots 57 121.9 387
with 15 snapshots 58 122.0 387
Table 4: Results for the elasticity problem with contrast . See Table 1 for the description of the preconditioners.
Preconditioner Iterations Condition Coarse dim.
None >2000
44 140.6 387
using twice many eigenvectors 26 15.1 774
141 387
110 387
58 1.8e2 387
58 140.8 387
with 10 snapshots 69 276.7 387
with 15 snapshots 69 276.5 387
with 10 snapshots 63 141.0 387
with 15 snapshots 63 141.1 387
Table 5: Results for the elasticity problem with contrast . See Table 1 for the description of the preconditioners.
Preconditioner 1
None 292 1583
14 26 53 44
14 28 56 58
25 33 62 69
20 30 58 63
Table 6: PCG iterations for different contrast values. See Table 1 for description of the preconditioners.
Preconditioner 1
None 3.2e3 1.2 1.2 4
4.6 16 113.8 140.6
4.6 18.3 140.8 140.8
15.6 27 111.4 276.5
9.4 22.8 122 141.1
Table 7: Spectral condition number of for different contrast values. See Table 1 for the description of the preconditioners.

6 Conclusions

In this paper we design and implement robust (with respect to the high-contrast and the multiscale structure) two-levels domain decomposition preconditioners for the elasticity equation appearing in topology optimization problems. Our design fits within the framework of the GMsFEM methodology where approximations of locally posed eigenvalue problems are used to construct the coarse space. We present several low cost constructions with similar number of iterations. The computational cost related to the construction of the preconditioners is reduced an order of magnitude. The presented numerical experiments demonstrate the quality and robustness of our iterations.


  • N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund (2017) Giga-voxel computational morphogenesis for structural design. Nature 550 (7674), pp. 84 – 86. External Links: Document Cited by: §1.2.
  • N. Aage and BoyanS. Lazarov (2013) Parallel framework for topology optimization using the method of moving asymptotes. Structural and Multidisciplinary Optimization 47 (4), pp. 493–505. External Links: Document, Link Cited by: §1.2.
  • J. Alexandersen and B. S. Lazarov (2015a) Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Computer Methods in Applied Mechanics and Engineering 290 (0), pp. 156 – 182. External Links: Document Cited by: item 2, §1.3, §1.3, §1.4.
  • J. Alexandersen and B. Lazarov (2015b) Tailoring macroscale response of mechanical and heat transfer systems by topology optimization of microstructural details. In Engineering and Applied Sciences Optimization, N. D. Lagaros and M. Papadrakakis (Eds.), Computational Methods in Applied Sciences, Vol. 38, pp. 267–288. External Links: Document Cited by: §1.3.
  • E. Andreassen, A. Clausen, M. Schevenels, B. Lazarov, and O. Sigmund (2011) Efficient topology optimization in matlab using 88 lines of code. Structural and Multidisciplinary Optimization 43, pp. 1–16. Note: 10.1007/s00158-010-0594-7 External Links: ISSN 1615-147X Cited by: §1.
  • O. Axelsson and I. Gustafsson (1978) Iterative methods for the solution of the navier equations of elasticity. Computer Methods in Applied Mechanics and Engineering 15 (2), pp. 241 – 258. Cited by: §4.
  • M. P. Bendsøe and O. Sigmund (2003) Topology optimization - theory, methods and applications. Springer Verlag, Berlin Heidelberg. Cited by: §1, §1, §1.
  • M. P. Bendsøe (1989) Optimal shape design as a material distribution problem. Structural optimization 1 (4), pp. 193–202. External Links: Document Cited by: §1.
  • R. Blaheta (1994) Displacement decomposition—incomplete factorization preconditioning techniques for linear elasticity problems. Numerical Linear Algebra with Applications 1, pp. 107–128. Cited by: §4.
  • D. Braess (2007) Finite elements: theory, fast solvers, and applications in solid mechanics. Cambridge University Press. Cited by: §1.1.
  • M. Buck, O. Iliev, and H. Andrä (2013) Multiscale finite element coarse spaces for the application to linear elasticity. Central European Journal of Mathematics 11 (4), pp. 680–701. External Links: ISSN 1895-1074, Document Cited by: §1.1, §1.3, §1.3.
  • M. Buck, O. Iliev, and H. Andrä (2014) Multiscale finite elements for linear elasticity: oscillatory boundary conditions. In Domain Decomposition Methods in Science and Engineering XXI, pp. 237–245. Cited by: §1.1, §1.3.
  • V. Calo, Y. Efendiev, J. Galvis, and G. Li (2016) Randomized oversampling for generalized multiscale finite element methods. Multiscale Modeling & Simulation 14 (1), pp. 482–501. External Links: Document Cited by: item 1, §3, §3.
  • Y. Efendiev and J. Galvis (2010) Eigenfunctions and multiscale methods for Darcy problems. Technical report ISC, Texas A& M University. Cited by: §1.3.
  • Y. Efendiev, J. Galvis, and T. Y. Hou (2013a) Generalized multiscale finite element methods (gmsfem). Journal of Computational Physics 251 (0), pp. 116 – 135. External Links: Document Cited by: item 1, §1.3, §1.3, §2.
  • Y. Efendiev, J. Galvis, and T. Y. Hou (2013b) Generalized multiscale finite element methods (gmsfem). Journal of Computational Physics 251 (0), pp. 116 – 135. External Links: Document Cited by: §2.1, §3.
  • Y. Efendiev, J. Galvis, R. Lazarov, and J. Willems (2012) Robust domain decomposition preconditioners for abstract symmetric positive definite bilinear forms. ESAIM: Mathematical Modelling and Numerical Analysis 46, pp. 1175–1199. Cited by: §1.2, §1.3, §2.1, §2.
  • Y. Efendiev, J. Galvis, and P. S. Vassilevski (2011a) Spectral element agglomerate algebraic multigrid methods for elliptic problems with high-contrast coefficients. In Domain Decomposition Methods in Science and Engineering XIX, Y. Huang, R. Kornhuber, O. Widlund, and J. Xu (Eds.), pp. 407–414. Cited by: §2.
  • Y. Efendiev, J. Galvis, and X. Wu (2011b) Multiscale finite element methods for high-contrast problems using local spectral basis functions. Journal of Computational Physics 230 (4), pp. 937 – 955. External Links: ISSN 0021-9991, Document Cited by: §2.1, §2.
  • Y. Efendiev and J. Galvis (2011) A domain decomposition preconditioner for multiscale high-contrast problems. In Domain Decomposition Methods in Science and Engineering XIX, Y. Huang, R. Kornhuber, O. Widlund, J. Xu, T. J. Barth, M. Griebel, D. E. Keyes, R. M. Nieminen, D. Roose, and T. Schlick (Eds.), Lecture Notes in Computational Science and Engineering, Vol. 78, pp. 189–196. External Links: Document Cited by: §1.3.
  • J. Galvis and Y. Efendiev (2010a) Domain decomposition preconditioners for multiscale flows in high contrast media: reduced dimension coarse spaces. Multiscale Modeling & Simulation 8 (5), pp. 1621–1644. External Links: Document Cited by: §1.2, §1.3.
  • J. Galvis, E. T. Chung, Y. Efendiev, and W. T. Leung (2017) On overlapping domain decomposition methods for high-contrast multiscale problems. In International Conference on Domain Decomposition Methods, pp. 45–57. Cited by: §3.
  • J. Galvis and Y. Efendiev (2010b) Domain decomposition preconditioners for multiscale flows in high-contrast media. Multiscale Modeling & Simulation 8 (4), pp. 1461–1483. External Links: Document Cited by: §1.2, §2.
  • M. Griebel, D. Oeltz, and M. Schweitzer (2003) An algebraic multigrid method for linear elasticity. SIAM Journal on Scientific Computing 25 (2), pp. 385–407. External Links: Document,, Link Cited by: §4.2.
  • I. Gustafsson and G. Lindskog (2002) On parallel solution of linear elasticity problems. part ii: methods and some computer experiments. Numerical Linear Algebra with Applications 9 (3), pp. 205–221. External Links: Document Cited by: item 2, §4.
  • I. Gustafsson and G. Lindskog (1998) On parallel solution of linear elasticity problems: part i: theory. Numerical Linear Algebra with Applications 5 (2), pp. 123–139. External Links: Document Cited by: item 2, §4.
  • N. Halko, P. G. Martinsson, and J. A. Tropp (2011) Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review 53 (2), pp. 217–288. External Links: Document Cited by: item 1, §3, §3.
  • B. S. Lazarov, F. Wang, and O. Sigmund (2016) Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics 86 (1), pp. 189–218. External Links: ISSN 1432-0681, Document Cited by: §1.
  • B. S. Lazarov (2014) Topology optimization using multiscale finite element method for high-contrast media. In Large-Scale Scientific Computing, I. Lirkov, S. Margenov, and J. Waśniewski (Eds.), Lecture Notes in Computer Science, pp. 339–346. External Links: Document Cited by: §1.3, §1.3.
  • Y. Saad (2003) Iterative methods for sparse linear systems. Vol. 82, siam. Cited by: §1.2, §5.
  • K. Svanberg (1987) The method of moving asymptotes - a new method for structural optimization. International Journal for Numerical Methods in Engineering 24, pp. 359–373. Cited by: §1.
  • A. Toselli and O. Widlund (2004) Domain decomposition methods - algorithms and theory. Springer. Cited by: §1.2, §1.3, §2.
  • N.P. van Dijk, K. Maute, M. Langelaar, and F. van Keulen (2013) Level-set methods for structural topology optimization: a review. Structural and Multidisciplinary Optimization 48 (3), pp. 437–472. Cited by: §1.
  • P. S. Vassilevski (2008) Multilevel block factorization preconditioners: matrix-based analysis and algorithms for solving finite element equations. Springer, New York. Cited by: §1.2.
  • F. Wang, B. Lazarov, and O. Sigmund (2011) On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization 43 (6), pp. 767–784. Cited by: §1.