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
(1)  
(2)  
(3)  
(4) 
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(5) 
where is the vector with the external forces applied on the system and is the socalled 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 socalled 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
(6) 
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 socalled robust formulation in topology optimization Wang et al. (2011), and other problems with highcontrast and multiscale coefficients like levelset 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
(7) 
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 NavierCauchy partial differential equation
(8) 
where
is the stress tensor,
is the strain tensor given by(9) 
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
(10) 
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
(11) 
with bilinear form and linear form
(12) 
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 illconditioned 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 Highcontrast coefficients in topology optimization
The focus here is on domain decomposition preconditioners, in particular, on twolevels domain decomposition techniques. The classical variants of these preconditioners Toselli and Widlund (2004) do not perform well for highcontrast problems Efendiev et al. (2012)
. The condition number estimates for the traditional domain decomposition case depends on the contrast
if the highconductivity 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 highconductivity 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 highconductivity 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 highcontrast 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 highcontrast, 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. Highcontrast 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 recomputation 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 highcontrast 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 twodimensions 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.
2 GMsFEM twolevels domain decomposition for the elasticity equation
The focus In this section is on the utilization of GMsFEM coarse spaces in the construction of twolevels domain decomposition preconditioners. In particular, we show that the proposed preconditioners yield a contrastindependent 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 nonoverlapping blocks . Based on the above decompotions an overlapping decomposition, shown in Figure 3, denoted as is obtained from an original nonoverlapping decomposition by enlarging each subdomain to
(13) 
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
(14) 
The part corresponding to the first level is
(15) 
where , , and the part corresponding to the second (or coarse) level is
where . The matrix consists of vectors defining the socalled 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 .
The finescale linear system Equation 5 is solved iteratively with the preconditioned conjugate gradient (PCG) method. The application of the preconditioner involves solving a coarsescale 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)
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
(16) 
with homogeneous Neumann boundary condition on and Dirichlet boundary condition on if is not empty. In matrix form we have
(17) 
with eigenvalues and the eigenvectors denoted as and respectively. The eigenvalues are ordered as
(18) 
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(19) 
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 highstiffness 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 highstiffness regions and its contrast dependent modes. There are three modes for each highstiffness 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 highcontrast 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 highstiffness 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.
For the elasticity problem the GMsFEM approximate the solution on the coarse grid as,
where are the unknown constants. The coarse matrix is constructed,
where
and the multiscale finite element solution is the finite element projection of the finescale 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 lowrank 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:
(20)  
(21) 
where every column of is a vector from . The approximations of the eigenvalues and the eigenvectors are computed as
(22)  
(23) 
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 ydirections, i.e., . Following the above decomposition, the stiffness matrix and the load vector , from Equation 5, can be written in a block form as
(24) 
The displacement split preconditioner matrix is constructed by keeping only the diagonal block matrices.
(25) 
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
(26) 
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 blockdiagonal preconditioner constructed for the diffusion equation and a coarselevel 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
(27) 
where
(28) 
and is the stiffness matrix obtained by a finite element discretization of the NavierStokes 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 componentwise (xand ydirection displacements) and each of those uses the local eigenvalue problem
(29) 
where or some other quantity that captures the highcontrast 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
(30) 
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 contrastdependent 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 finegrid 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
(31) 
where is the enriched coarse space and
(32) 
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,
(33) 
and
(34) 
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 
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 highcontrast material distribution as shown in Figure 6. The computational domain is partitioned applying coarse mesh and each coarseelement is further partitioned utilizing finemesh. 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 solidmaterial 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 levelone 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 25 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 contrastasymptoticallyvanishing eigenvalues and the corresponding eigenvectors.
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 
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 
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 
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 
Contrast  

Preconditioner  1  
None  292  1583  
14  26  53  44  
14  28  56  58  
25  33  62  69  
20  30  58  63 
Contrast  

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 
6 Conclusions
In this paper we design and implement robust (with respect to the highcontrast and the multiscale structure) twolevels 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.
References
 Gigavoxel computational morphogenesis for structural design. Nature 550 (7674), pp. 84 – 86. External Links: Document Cited by: §1.2.
 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.
 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.
 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.
 Efficient topology optimization in matlab using 88 lines of code. Structural and Multidisciplinary Optimization 43, pp. 1–16. Note: 10.1007/s0015801005947 External Links: ISSN 1615147X Cited by: §1.
 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.
 Topology optimization  theory, methods and applications. Springer Verlag, Berlin Heidelberg. Cited by: §1, §1, §1.
 Optimal shape design as a material distribution problem. Structural optimization 1 (4), pp. 193–202. External Links: Document Cited by: §1.
 Displacement decomposition—incomplete factorization preconditioning techniques for linear elasticity problems. Numerical Linear Algebra with Applications 1, pp. 107–128. Cited by: §4.
 Finite elements: theory, fast solvers, and applications in solid mechanics. Cambridge University Press. Cited by: §1.1.
 Multiscale finite element coarse spaces for the application to linear elasticity. Central European Journal of Mathematics 11 (4), pp. 680–701. External Links: ISSN 18951074, Document Cited by: §1.1, §1.3, §1.3.
 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.
 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.
 Eigenfunctions and multiscale methods for Darcy problems. Technical report ISC, Texas A& M University. Cited by: §1.3.
 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.
 Generalized multiscale finite element methods (gmsfem). Journal of Computational Physics 251 (0), pp. 116 – 135. External Links: Document Cited by: §2.1, §3.
 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.
 Spectral element agglomerate algebraic multigrid methods for elliptic problems with highcontrast 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.
 Multiscale finite element methods for highcontrast problems using local spectral basis functions. Journal of Computational Physics 230 (4), pp. 937 – 955. External Links: ISSN 00219991, Document Cited by: §2.1, §2.
 A domain decomposition preconditioner for multiscale highcontrast 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.
 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.
 On overlapping domain decomposition methods for highcontrast multiscale problems. In International Conference on Domain Decomposition Methods, pp. 45–57. Cited by: §3.
 Domain decomposition preconditioners for multiscale flows in highcontrast media. Multiscale Modeling & Simulation 8 (4), pp. 1461–1483. External Links: Document Cited by: §1.2, §2.
 An algebraic multigrid method for linear elasticity. SIAM Journal on Scientific Computing 25 (2), pp. 385–407. External Links: Document, https://doi.org/10.1137/S1064827502407810, Link Cited by: §4.2.
 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.
 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.
 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.
 Length scale and manufacturability in densitybased topology optimization. Archive of Applied Mechanics 86 (1), pp. 189–218. External Links: ISSN 14320681, Document Cited by: §1.
 Topology optimization using multiscale finite element method for highcontrast media. In LargeScale 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.
 Iterative methods for sparse linear systems. Vol. 82, siam. Cited by: §1.2, §5.
 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.
 Domain decomposition methods  algorithms and theory. Springer. Cited by: §1.2, §1.3, §2.
 Levelset methods for structural topology optimization: a review. Structural and Multidisciplinary Optimization 48 (3), pp. 437–472. Cited by: §1.
 Multilevel block factorization preconditioners: matrixbased analysis and algorithms for solving finite element equations. Springer, New York. Cited by: §1.2.
 On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization 43 (6), pp. 767–784. Cited by: §1.
Comments
There are no comments yet.