Dynamic Adaptive Mesh Refinement for Topology Optimization
aa r X i v : . [ m a t h . NA ] S e p Dynamic Adaptive Mesh Refinement for Topology Optimization ∗ Shun Wang † , Eric de Sturler ‡ , Glaucio H. Paulino § Abstract
We present an improved method for topology optimization with both adaptive mesh refinement andderefinement. Since the total volume fraction in topology optimization is usually modest, after a fewinitial iterations the domain of computation is largely void. Hence, it is inefficient to have many smallelements, in such regions, that contribute significantly to the overall computational cost but contributelittle to the accuracy of computation and design. At the same time, we want high spatial resolutionfor accurate three-dimensional designs to avoid postprocessing or interpretation as much as possible.Dynamic adaptive mesh refinement (AMR) offers the possibility to balance these two requirements. Wediscuss requirements on AMR for topology optimization and the algorithmic features to implement them.The numerical design problems demonstrate (1) that our AMR strategy for topology optimization leadsto designs that are equivalent to optimal designs on uniform meshes, (2) how AMR strategies that donot satisfy the postulated requirements may lead to suboptimal designs, and (3) that our AMR strategysignificantly reduces the time to compute optimal designs. keywords: adaptive mesh refinement, topology optimization, iterative solvers.
Topology optimization is a powerful structural optimization method that combines a numerical solutionmethod, usually the finite element method, with an optimization algorithm to find the optimal materialdistribution inside a given domain [16, 20, 15, 14]. In designing the topology of a structure we determinewhich points in the domain should be material and which points should be void. However, it is known that,in the continuum setting, topology optimization leads to designs with intermediate densities. So continuousvalues between 0 and 1 replace discrete values (0 or 1) to represent the relative densities, and some form ofpenalization is used to obtain designs with almost discrete 0/1 material density distribution [3].In topology optimization, problems are solved most commonly on fixed uniform meshes with a relativelylarge number of elements in order to achieve accurate designs [17, 12]. However, as void and solid regionsappear in the design, it is more efficient to represent the holes with fewer large elements and the solid regions,especially the material surface, with more fine elements. Since the shape and position of holes and solid ∗ This work was supported in part by the National Science Foundation under Grant DMR-03 25939 ITR through the MaterialsComputation Center at the University of Illinois at Urbana-Champaign. † Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, U.S.A., [email protected], ‡ Department of Mathematics, Virginia Tech, Blacksburg, Virginia 24061, U.S.A., [email protected], § Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801,U.S.A., [email protected]. a priori .Therefore, adaptive mesh refinement (AMR) is very suitable for topology optimization.
The purpose of AMRfor topology optimization is to get the design that would be obtained on a uniformly fine mesh, but at a muchlower computational cost by reducing the total number of elements and having fine elements only where andwhen necessary .Highly accurate designs on uniform meshes may require so many elements that the solution of theoptimization problem becomes intractable. However, AMR leads to high resolution in the mesh only whenand where necessary. This makes it possible to obtain accurate designs with a modest number of elementsand hence with a reasonable cost. Even when a design on a uniform mesh is computationally feasible, AMRtends to reduce the computational cost by reducing the required number of elements and by improving theconditioning of linear systems arising from the finite element discretization. Obviously, we do not want theuse of AMR or the AMR procedure to alter the computed designs. However, there is a risk of this, sincethe mesh influences the computed deformations and sensitivities. It is imperative then that the solutionsfrom the finite element analysis using AMR must be as accurate as those obtained on a uniform fine mesh. Moreover, the final design must be governed by accurate sensitivities corresponding to those obtainedon the finest mesh. If coarse mesh solutions drive or limit the design, suboptimal designs may result whendesigns optimal on a coarser mesh differ substantially from the optimal design on a (much) finer mesh. Wewill demonstrate that this occurs in quite simple cases. The early work in this area, though leading toacceptable designs in specific instances, does not satisfy these properties. We will propose relatively simplebut essential changes to these methodologies that lead to AMR-based designs that are equivalent (up tosome small tolerance) to designs on uniform fine meshes. In addition, our approach leads (in principle) to anoverall more efficient method as we reduce the total number of elements further. The topology optimizationmay lead to a sequence of (intermediate) structures requiring high mesh resolution in different parts of thecomputational domain. Therefore, it is important to (1) allow the meshes at all levels to change continually(dynamic) and (2) to allow both mesh refinement and derefinement [23]. Derefinement is important forefficiency when the initial discretization needs to include relatively small elements in certain regions. Thisis important in a number of cases, which are elaborated upon below.In the next section, we provide an assessment of previous AMR strategies, namely the implementations byCosta and Alves [8] and Stainko [22]. In Section 3, we provide a brief introduction to topology optimization.Next, in Section 4, we state the purpose of our AMR strategy for topology optimization and explain therequirements it poses. Based on these requirements, we propose a more robust and dynamic AMR strategy.We describe a number of implementation issues of our AMR strategy in Section 6. We briefly discussthe iterative solution of the large sparse linear systems arising in the finite element analysis in Section 5.In Section 7, we show numerical experiments that demonstrate the robustness and efficiency of our AMRstrategy. The first experiment also explains why the refinement strategies by Costa and Alves [8] and Stainko[22] may lead to suboptimal designs. Finally, in Section 8, we present conclusions about our AMR strategyfor topology optimization algorithms, and we mention some directions for future work.
Little research has been done in applying AMR to topology optimization. So, we start by briefly discussingtwo recent, important, papers in this area. The AMR method by Costa and Alves [8] goes through a prede-termined, fixed sequence of optimizations and subsequent mesh refinements (they do not use derefinements), For comparison everywhere in this paper the element size for the uniform fine mesh is the same as the element size at thehighest level of refinement in the AMR mesh. . . ptimization (Update design variables)InitializationFinite Element AnalysisSensitivity AnalysisOptimal TopologyNo Yes Iteration 0 Iteration 4 Iteration 16 Iteration 106 (Final iteration)Convergence?FilteringAMR Figure 1: Overview of the topology optimization algorithm with dynamic AMR.coarse level designs and refines the mesh to accommodate such changes will be inefficient if fine elements invoid regions cannot be removed.Therefore, a more robust and efficient refinement strategy is needed. Hence, we propose a dynamicmeshing strategy that includes both mesh refinement and derefinement everywhere in the computationaldomain. Our improved AMR strategy has two main components. First, we extend the refinement criteriafrom [8], refining all material elements and elements on the boundary, but with an additional layer ofrefinements around the material boundary (in the void-region). The thickness of the layer is a parameter.This way the fine level design can change shape arbitrarily in optimization steps between mesh refinements.Second, our AMR method updates coarse and fine meshes continually, so that small changes arising in themore accurate computations on finer meshes can change the shape of the design arbitrarily in the courseof the optimization; the fine(r) meshes move with the material boundary. This means that our designs arereally based on the accurate fine mesh computations and are not confined to regions fixed by earlier coarsemesh computations. Since we do continual mesh adaptation, we may have fine elements in regions thathave become void at some point. Derefinement will remove those fine elements. Further details are given inthe next subsection. This approach leads to designs on AMR meshes at greatly reduced cost that are thesame as the designs that would have been obtained on a uniform fine mesh of the highest resolution (withina small tolerance). We will demonstrate this experimentally in section 7. Our approach also allows us tostart with coarser meshes, since the coarse mesh solution does not need to be a good approximation to thefinal solution. Even when we start with a finer mesh for faster convergence or proper functioning of theregularization filter, derefinement allows us to remove fine elements that have become void (which tends tohappen quite quickly). 4
A Brief Topology Optimization Review
In topology optimization we solve for the material distribution in a given computational design domain Ω.The topology optimization problem we consider here is to minimize the compliance of a structure undergiven loads as a function of the material distribution. To solve this problem numerically, we discretize thecomputational domain using finite elements, where we usually use a lower order interpolation for the densityfield (material distribution) than for the displacement field. The most common approach (also employedfor this paper) is to use (bi-,tri-)linear interpolation for the displacement field and constant density in eachelement. The compliance minimization problem after finite element discretization is defined asmin ρ e ∈ [ ρ o , , ∀ e f T u (1)s.t. K ( ρ ) u = f for x ∈ Ω \ Ω , u = u for x ∈ Ω , P e ρ e V e ≤ V , where ρ e is the density in element e , ρ is the vector of element densities, K is the stiffness matrix, a functionof the discretized density field ( ρ ), V e is the volume of element e , V is a maximum volume (fraction) allowedfor the design, and Ω is the part of the domain where the displacement is prescribed. To avoid singularityof the stiffness matrix, we enforce a small positive lower bound ρ o on the element density, typically 10 − .As mentioned in the introduction, our discrete model must drive the continuous material distribution asmuch as possible to a 0/1 solution. We use the Solid Isotropic Material with Penalization (SIMP) methodto make the undesirable intermediate densities between ρ o (replacing 0) and 1 unfavorable [6]. In this case,the elasticity tensor is defined as a function of the element density, E e = ρ pe E , (2)where p is the penalization parameter. With p >
1, intermediate densities are unfavorable as they providerelatively little stiffness compared with their material cost. A common choice is p = 3, which results inintermediate material properties that satisfy the Hashin–Shtrikman bound for composite materials [10]. Toavoid problems with local minima, we usually apply continuation on the parameter p , that is, we start with p = 1 and slowly increase p as the design converges.The general scheme for topology optimization using AMR is illustrated in Figure 1. First, we set up thegeometry, the finite element (FE) mesh, the loading and boundary conditions, and we initialize the densitydistribution ρ . Then, we start the optimization loop. In this loop, we assemble and solve the equilibriumequations K ( ρ ) u = f in (1) using the FE discretization and a linear solver. Next, in the sensitivity analysis,we compute the derivatives of the objective function with respective to the design variables, ∂c/∂ρ e . Afterthis, we can apply an optional low-pass filter to remedy the checkerboard problem [18, 19, 21], which can bealso addressed by an alternative minimum length scale approach [9]. In the next step, we compute an updateof the design variables. There are various optimization algorithms applicable to topology optimization. Forthis paper, we use Optimality Criteria (OC), a simple approach based on a set of intuitive criteria [6, 4].After updating the design variables using a chosen optimization algorithm, we check the convergence of thedesign. Under certain conditions, to be discussed next, dynamic mesh adaptation is carried out before the(next) finite element analysis. 5 A Dynamic AMR Strategy
We base our algorithmic choices on a set of requirements on AMR codes for topology optimization. Asstated above, the purpose of AMR for topology optimization is to get the design that would be obtained ona uniform fine mesh, but at a much lower computational cost by reducing the total number of elements andhaving fine elements only where (and when) necessary.First, since the finite element analysis and the computation of sensitivities drive the changes in materialdistribution, they should be as accurate as on the uniform fine mesh. Therefore, we need a fine mesh thatcovers at least the material region and the boundary. Since the void regions have negligible stiffness they donot influence the (intermediate) linear finite element solutions and resulting sensitivity computations. Thuswe do not need a fine mesh inside the void region, and we can use a refinement criterion similar to that ofCosta and Alves [8]. At this point we focus on refinement and derefinement for shape only. Therefore, we areconservative with respect to accuracy, and we expect that, in future implementations, good error indicatorswill lead to further efficiency gains, in particular because of derefinement in solid material regions.Second, the accurate computations on the finest level should drive the changes in the material distribution.This requires continual mesh adaptation so that computational results after refinements can drive updatesto the material distribution, and designs are not confined by earlier coarse grid results. This also means thatas the material region moves close to the boundary between fine and coarse(r) mesh, additional refinementsallow for further evolution.Third, we need to ensure that the design can change sufficiently in between mesh updates. Therefore,we maintain a layer of additional refinements around the material region (in the void region) and carry outcontinual mesh adaptation. Due to the additional layer of refinements and continual mesh updates, the designcan change arbitrarily following the fine grid computations and resulting sensitivities, and it is not confinedby earlier coarse grid results. To ensure that the design accurately reflects the fine mesh computations, weallow rapid refinements of the mesh early on when voids and material regions (and hence the boundary)develop, rather than delay refinements until later stages, when a suboptimal design might have developed.Fourth, since the design can change substantially from its estimate on a coarse mesh, we may have fineelements in void regions. Those elements must be removed for efficiency, requiring derefinements of the mesh.To facilitate our strategy of continual mesh refinement and derefinement, we use a hierarchical representationof adaptive meshes.We will now state our refinement and derefinement strategy in more detail. We adapt the mesh whenCASE (i):1. the relative change in the compliance is smaller than a given threshold, and
2. a given minimum number of optimization steps have been carried out since the last mesh update, or whenCASE (ii):3. a given maximum number of optimization steps have been carried out without meeting conditions 1and 2.Condition 1 corresponds to a common convergence criterion for topology optimization that the maximumchange in the design variables is smaller than a certain tolerance (which we usually set to 0.01). Thiscondition is satisfied when the solution is near a local minimum, which might be caused by a no longerappropriate mesh. In that case, we must adapt the mesh to allow the design to change further. If the localminimum is not an artifact of the mesh, the design will remain the same after mesh adaptation. Condition 26 a r amr b b r amr Figure 2: Refinement criteria for void element. Element a is marked for refinement, because it has solidelements within distance r AMR ; element b is marked for derefinement.prevents refinement and derefinement from happening too frequently. This is important as the solution needsto adapt to the changed mesh, so that the computed sensitivities reflect the design and not mesh artifacts.This also limits the cost of mesh adaptation. In our experiments this minimum number of optimization stepsis set to five based on experience, and other values are possible. Regarding condition 3, in our experimentswe adapt the mesh at least every ten optimization steps. This condition leads to faster convergence becauseit ensures that the mesh is appropriate for the material distribution. Using these conditions, we can startwith a fairly coarse mesh, and we may carry out mesh (de)refinement before the design converges on anymesh if necessary.We adapt our mesh according to the following procedure.1. Mark all the elements for refinement or derefinement based on the following criteria: • If element e is solid, i.e., ρ e ∈ [ ρ s , ρ s is a chosen density threshold, or element e is withina given distance r AMR from a solid element we mark it for refinement. • If element e is void, i.e., ρ e ∈ [ ρ o , ρ s ], and there are no solid elements within distance r AMR , wemark element e for derefinement. See Figure 2.2. Check compatibility for the mesh that will be generated and make the following adjustments in twosweeps over all elements: • In the first sweep, we unmark elements marked for derefinement, if they have a sibling (an elementgenerated by the same refinement) that is not marked for derefinement. • In the second sweep, we unmark elements marked for derefinement, if derefinement would leadto level two or higher edge incompatibility. We allow only level one incompatibility; see Figure 3and the discussion in Section 6.The above refinement criteria result in a layer of fine elements on the void side of the solid/void interfacethat allows the material to be redistributed locally. If a material boundary moves near the fine/coarseelement interface, mesh refinement creates a new layer of fine elements around the current material surfaceto allow further local redistribution of the material. On the other hand, if some fine elements become void,these fine elements are removed by derefinement to keep the optimization efficient.7
Iterative Solution Scheme
Although AMR significantly reduces the number of DOFs in the finite element simulation, we still have tosolve a sequence of large linear systems, especially for three-dimensional designs. Moreover, because of thelarge difference in density between distinct parts of the computational domain and the elasticity tensor givenby (2), with p = 3 toward the end of the optimization, the linear systems are very ill-conditioned. Hence,proper preconditioning is essential. In [24], we showed how to precondition the linear systems arising intopology optimization, and we also used Krylov subspace recycling [13] to reduce the number of iterationsover multiple linear systems. We briefly mention the main ideas here.To remedy the serious ill-conditioning in topology optimization problems, we explicitly rescale eachstiffness matrix such that the diagonal coefficients are all equal, as is the case for a problem with homogeneousdensity. We rescale the stiffness matrix K by multiplying with a diagonal matrix on both sides,˜ K = D − / KD − / , where D is the diagonal of K . The importance of such scaling and why it helps has been explained for anidealized one-dimensional problem in [24]. We further reduce the condition number of the system matrixand the number of iterations for convergence by applying an incomplete Cholesky preconditioner with zerofill-in to the explicitly rescaled system, ˜ K ≈ ˜ L ˜ L T . The finite element analysis in topology optimization requires the solution of a sequence of (usually)symmetric linear systems. In each optimization step, the algorithm updates the element densities, and afterthe first few optimization steps the changes in the design variables tend to be small from one optimizationstep to the next. Hence, the optimization leads to small changes from one linear system to the next, and thesearch space generated for one linear system provides important information for subsequent linear systems.First, the solution of one system can be used as an initial guess for the next system, reducing the initialresidual. Second, an approximate invariant subspace derived from the Krylov space generated for one linearsystem can be used for subsequent linear systems, improving the convergence rate of the iterative solver.This is the basic idea of Krylov subspace recycling; however, other subspaces may also be used for ’recycling’[13]. Since the linear systems discussed in this paper are symmetric, we use the Recycling MINRES algorithm(RMINRES) for Krylov subspace recycling [24]. Unfortunately, solving a sequence of problems on meshesthat change periodically makes recycling more difficult. Although it is not hard to map relatively smootheigenvectors from a mesh to an updated mesh, the combination of mesh adaptation and preconditioningseems to give accuracy problems. Recycling is still effective for AMR; however, it is not nearly as beneficialas on a static mesh, and its improvement is work in progress.
For the implementation of adaptive mesh refinement, we use the libMesh library [11] developed at theUniversity of Texas at Austin and the Technische Universit¨at Hamburg-Harburg. The libMesh libraryprovides a C++ framework for numerical simulations of partial differential equations on serial and parallelplatforms. It supports one-dimensional, two-dimensional, and three-dimensional finite element and finitevolume simulations on adaptive meshes. The libMesh software uses
PETSc [2, 1] for the solution of linearsystems on both serial and parallel platforms. However, we use our own custom linear solvers with Krylov8ubspace recycling and preconditioners as detailed in Section 5 and references [24, 13]. For compatibilitywith the libMesh package we have implemented the RMINRES method in the
PETSc framework. For theincomplete Cholesky preconditioner we used routines provided by
PETSc .We have developed two-dimensional and three-dimensional topology optimization algorithms on top of libMesh . Currently, we use element-based design variables, the SIMP method for material interpolation [5],the OC method for optimization [6, 4], and Sigmund’s filter technique [18, 19, 21] with some modifications fordealing with adaptive refinement. Following Stainko [22], we make a small modification in Sigmund’s filterfor a nonuniform mesh. Sigmund’s filter takes a distance and density weighted average of the sensitivities ofall elements in a certain radius as d ∂c∂ρ e = 1 ρ e P d H de X d ρ d H de ∂c∂ρ d , (3)where ∂c/∂ρ e is the sensitivity of the compliance with respect to the density of element e , and H de is adistance weight defined as H de = max { r min − dist( d, e ) , } . (4)The parameter r min is a given radius for the filter (for the work reported here we use r min = r AMR ), anddist( d, e ) is the distance between the centers of elements d and e . For a nonuniform mesh, we take thevariation of element size into account by using the element volume to redefine the weight in the filter [22] as d ∂c∂ρ e = 1 ρ e P d H de V d X d ρ d H de V d ∂c∂ρ d . (5)The filter radius r min is often a length scale independent of the mesh representation. Notice that the filterwill be effectively deactivated if its size is smaller than that of the smallest element, i.e., no element has anyneighbors within distance r min .Because of the hierarchical data structure of libMesh , we must allow level-one mesh incompatibility.However, we avoid level-two and higher mesh incompatibility. For example, if the configuration in Figure 3(b)would result from mesh refinement, we refine the gray elements as well. If the configuration, in particularthe gray elements, in Figure 3(b) would result from mesh derefinement, we avoid the derefinement. In thisway, we limit mesh incompatibility to level-one mesh incompatibility. As indicated by the circled nodes inFigure 3(a), level-one mesh incompatibility results in hanging nodes. The libMesh package handles thosehanging nodes by using the projection method to enforce constraints in the stiffness matrix. We divide thedegrees of freedom (DOFs) into two groups. Group one consists of all the unconstrained DOFs, and grouptwo consists of the constrained DOFs on the hanging nodes. The constrained DOFs can be computed bylinear interpolation from unconstrained DOFs. If we define vector ˜ u on the unconstrained DOFs, then u = (cid:18) ˜ uP ˜ u (cid:19) = (cid:20) I P (cid:21) (cid:18) ˜ u (cid:19) (6)is the mapping of ˜ u to all the DOFs, where P is the interpolation matrix. We compute ˆ u by solving theprojected system (cid:20) I P T (cid:21) K (cid:20) I P (cid:21) ˆ u = (cid:20) I P T (cid:21) f . (7)Since libMesh does not drop the constrained DOFs in the linear system, the projected system in (7) issingular when there is any hanging node. Krylov subspace methods can handle such singularities as long as9 c bc bc bc (a) bcbc bc bc bcbc (b)Figure 3: Mesh incompatibility with examples of quad, triangle and hex elements: (a) level-one meshincompatibility marked by circled nodes; (b) level-two mesh incompatibility marked by circled nodes. Weallow level-one mesh incompatibility (see Section 6), but we avoid level-two and higher incompatibility byrefining the gray coarse elements and by not derefining their children elements if these gray elements resultfrom a potential derefinement.the right hand side is consistent, but these singularities may cause problems for preconditioners. To avoidthe singularities, we set the diagonal entries in the matrix corresponding to the constrained DOFs to 1 andsolve (cid:18)(cid:20) I P T (cid:21) K (cid:20) I P (cid:21) + (cid:20) I (cid:21)(cid:19) ˆ u = (cid:20) I P T (cid:21) f . (8)In the end, we recover the constrained DOFs using the interpolation matrix: u = (cid:20) I P (cid:21) ˆ u . (9) We solve three problems to demonstrate the improvements of our new AMR scheme and verify that thecomputed designs using AMR meshes are equivalent to designs on uniform fine meshes.For the first (2D) test problem, we compute the optimal design on a uniform fine mesh and on anadaptively refined mesh with both our AMR scheme and an approach following references [8, 22]. Thehighest level of refinement in the AMR meshes has the same element size as that in the uniform mesh.The results show that our scheme computes a solution equivalent to the optimal design on the uniform finemesh (within a small tolerance), while the alternative AMR approach from [8, 22] does not. Moreover, theexperiments elucidate how this suboptimal design arises from the strategy to only refine the results from afixed coarse mesh.For the second (3D) test problem, we compare the optimal design using an adaptive mesh and our AMRstrategy with the optimal design on a uniform fine mesh for a relatively simple cantilever beam problem.Again the maximum refinement level for the AMR mesh has the same element size as the uniform, fine mesh.10e also use this test to show that our AMR strategy leads to faster convergence in both the linear iterations(finite element solution) and the nonlinear iterations (design optimization) as well as to a significant reductionin runtime (about a factor of three).In the third test problem (also 3D), we compare the optimal design using an adaptive mesh and ourAMR strategy with the optimal design on a uniform, fine mesh for a more complicated design problem. Forall three test problems our AMR strategy leads to essentially the same design as is obtained on a uniformmesh, but at significantly reduced cost.To evaluate the relative difference between two designs, we need a quantitative measure. We define thethe relative difference between two designs as D ( ρ (1) , ρ (2) ) = R Ω | ρ (1) − ρ (2) | d Ω R Ω ρ (1) d Ω . (10)We take ρ (1) to indicate the design on the uniform fine mesh, and ρ (2) to indicate the design on the AMRmesh. This difference can be computed in a number of ways. To simplify comparison, we take the uniformmesh to be the AMR mesh with maximum refinement at every point. So, we refine the AMR mesh to thesame fine mesh level at every point (without changing the design), and then evaluate the 1-norm of thedifference between the designs. Test 1: 2D cantilever beam
We compute the optimal design for the 2D beam problem shown in Figure 4(a). We first compute the designon a uniform, fine mesh. Figure 4(b) shows an intermediate result, and Figure 4(c) the converged design.Note how the truss member at the lower-right corner has risen up noticeably from the intermediate resultto the final design. An effective AMR procedure must be able to capture such an evolution in the design.Next, we solve the same problem following the strategy mentioned in references [8, 22]. We start witha relatively coarse mesh (64 × (b) (c)Figure 4: Topology optimization on a 256 × uniform mesh : (a) problem configuration (the volumeconstraint V is 50% of the domain volume); (b) an intermediate design; (c) final converged design.12 .0010.010.11 (a) (b) (c)Figure 5: Topology optimization on an adaptive mesh with only refinement on each level : (a) convergedresult on the coarsest mesh with 2048 elements; (b) converged result on the intermediate mesh with 5675elements; (c) converged result on the final mesh with 20216 elements. Note the undesirable position of thetruss member near the lower right corner, which remains nearly invariant during the evolution of the design.13 .0010.010.11 (a) (b) (c)Figure 6: Topology optimization on an adaptive mesh with continual dynamic refinement and derefinementon each level : (a)–(b) intermediate designs; (c) final converged design on a nonuniform mesh with 25229elements, whose finest resolution is the same as the resolution of the uniform mesh in Figure 4. Notice thatthe truss member of the lower-right corner moves up as the AMR procedure progresses.14able 1: Mesh adaptation scheme following Costa and Alves [8].opt. step ℓ max L difference80 1 2048 4290 21.16%125 2 5675 11948 19.53%200 3 20216 41824 19.42%300 3 20216 41824 19.60%400 3 20216 41824 19.66%500 3 20216 41824 19.66%600 3 20216 41824 19.65%700 3 20216 41824 19.67%800 3 20216 41824 19.66%900 3 20216 41824 19.65%1000 3 20216 41824 19.67%Table 2: Dynamic AMR schemeopt. step ℓ max L difference27 1 2048 4290 21.13%37 2 6113 12786 19.80%100 3 24707 50560 17.77%200 3 25040 51242 13.00%300 3 25184 51550 5.00%367 3 25229 51654 0.17%15igure 7: 3D cantilever beam example with domain scale 2:1:1.evolution of intermediate designs on the uniform mesh (reference). The figures also demonstrate how themesh changes smoothly with the changes in material distribution. The smallest element size in the AMRmeshes in Figure 6 is the same as the element size for the uniform mesh shown in Figure 4 (reference).Compared with the final solution shown in Figure 5, the solution obtained with our AMR strategy is closerto the solution obtained on the uniform mesh (reference). Indeed, based on the metric in (10), the relativedifference between the designs in Figure 5(c) and Figure 4(c) is 19 . . ℓ max refers to the highest refinement level present in the mesh. Test 2: 3D cantilever beam
We compute the optimal design for the three-dimensional cantilever beam, shown in Figure 7, with a volumeconstraint of 25%. Exploiting symmetry, we discretize only a quarter of the domain. We solve this problemon a (fixed) uniform mesh with 128 × ×
32 B8 elements and also following our AMR strategy. The initialmesh for the AMR-based design has 64 × ×
16 B8 elements. The final results are shown in Figure 8with Figure 8(a) displaying the solution on a uniform, fine mesh, and Figure 8(8) displaying the AMRsolution; note the large blocks in parts of the void region in Figure 8(b). The relative difference betweenthese two designs is only 0 . Test 3: Cross-shaped domain
We compute the optimal design for the more complex three-dimensional test problem shown in Figure 10.For the cross-shaped domain, we compute the optimal design subject to the fixed boundary on the bottomfront and back ends, and two loads on the left and right sides. The maximum volume allowed is 20% ofthe domain volume. We solve this problem both on a uniform mesh and on an adaptive mesh following our16a)(b)Figure 8: Final solutions of the 3D cantilever beam problem (Figure 7) obtained using symmetry on aquarter of the domain as indicated by the mesh: (a) final solution on a fixed uniform mesh with 128 × ×
20 40 60 80 100 120012345 x 10 un k no w n s number of unknowns in FE systems uniformAMR (a) i t e r a t i on s number of iterations for linear solvers uniform − MINRESuniform − RMINRES(200,10)AMR − MINRESAMR − RMINRES(200,10) 0 20 40 60 80 100 120050100150200250300 s e c ond s timing for linear solvers uniform − MINRESuniform − RMINRES(200,10)AMR − MINRESAMR − RMINRES(200,10) (b) (c)Figure 9: Comparison of linear solver statistics for the cantilever beam design problem on a uniform, finemesh and on an adaptive mesh: (a) the number of unknowns in the linear systems arising from the finiteelement discretization; (b) the number of preconditioned MINRES and RMINRES(200,10) iterations (seebelow) for each optimization step; (c) solution times with MINRES and RMINRES(200,10) for the linearsystems arising from finite element discretization at each optimization step. The parameters m and k in RMINRES(m,k) have the following meaning. The method recycles an approximate invariant subspaceassociated with the smallest k eigenvalues from one linear system to the next. In the solution of single linearsystem, the approximate invariant subspace is updated every m iterations.18 b Figure 10: A 3D compliance minimization problem in a cross-shaped domain with the bottom (shaded) frontand back ends fixed and the bottom left and right ends pulled down. The volume constraint V is 20% ofthe domain volume.AMR strategy. The results are shown in Figures 11 and 12, respectively. The uniform mesh consists of 40960B8 elements, while the final adaptive mesh consists of only 19736 B8 elements. Moreover, the optimizationconverges in over 200 steps on the uniform mesh, but in only 106 optimization steps on the adaptive mesh.The adaptive mesh refinement reduces the total solution time by more than a factor of three to about 30%of the solution time for the uniform mesh. Nonetheless, the relative difference between these two designs isonly 2 .
58% (Eq. (10)).
In order to reduce the high computational cost of accurate three-dimensional designs by topology optimiza-tion we use adaptive mesh refinement. We propose several critical improvements to the approaches proposedby Costa and Alves [8] and Stainko [22] in order to attain better designs. In particular, we want to obtainthe same optimal designs that would be obtained on a uniform, fine mesh with AMR discretization hav-ing significantly fewer elements but the same fine mesh resolution. The purpose of AMR is to reduce thecost for the (same) optimal design; we do not want to reduce the quality of designs. For large, complex,three-dimensional design problems we could not possibly use a uniform fine mesh at the desired resolution.Our approach requires a dynamic meshing strategy that involves continual refinements and derefinementsfollowing the strategy laid out in Section 4. Derefinements should also lead to further efficiency improve-ment by reducing the number of elements in void regions, especially for three-dimensional problems. Usingthree test problems, we demonstrate that our AMR algorithm achieves the desired designs that are withina small tolerance of those obtained on a uniform, fine mesh with the same finest level. Our AMR strategysignificantly reduces the total runtime, nonlinear, and linear iterations with respect to using uniform meshes.Important future work includes error estimation in the finite element analysis and mesh refinement andderefinement governed by both considerations of accurate design and error estimation. In addition, we planto work on preconditioners that can be adapted with the mesh (rather than recomputed) and to improve theconvergence rate of Krylov methods with subspace recycling [24, 13]. We also intend to extend the presentAMR technique to multiphysics problems [7]. 19igure 11: The optimal solution to the design problem shown in Figure 10 on a uniform finite element meshwith 40960 B8 elements. The quarter-mesh discretization is shown on the bottom figure.20igure 12: The optimal solution to the design problem shown in Figure 10 on an adaptively refined mesh.The final mesh consists of 19736 B8 elements. The quarter-mesh discretization is shown on the bottomfigure. 21
Acknowledgements
We are indebted to Hong Zhang and Mat Knepley from Argonne National Laboratory for their help withthe
PETSc library and to Roy Stogner, John Peterson, and Benjamin Kirk from the University of Texas atAustin for their help with the libMesh library. We also thank Cameron Talischi for insigthful discussionswhich contributed to improve the present manuscript.
References
Structural Optimization , (4):193–202, 1989.[4] M. P. Bendsøe and N. Kikuchi. Generating optimal topologies in structural design using a homogeniza-tion method. Computer Methods in Applied Mechanics and Engineering , (2):197–224, 1988.[5] M. P. Bendsøe and O. Sigmund. Material interpolation schemes in topology optimization. Archives ofApplied Mechanics , (9–10):635–654, 1999.[6] M. P. Bendsøe and O. Sigmund. Topology Optimization: Theory, Methods and Applications . Springer-Verlag, Berlin, 2003. ISBN 3-540-42992-1.[7] R. C. Carbonari, E. C. N. Silva, and G. H. Paulino. Topology optimization design of functionallygraded bimorph-type piezoelectric actuators.
Smart Materials and Structures , :2605–2620, 2007.DOI 10.1088/0964-1726/16/6/065.[8] J. C. A. Costa Jr. and M. K. Alves. Layout optimization with h -adaptivity of structures. InternationalJournal for Numerical Methods in Engineering , (1):83–102, 2003.[9] J. K. Guest, J. H. Pr´evost, and T. Belytschko. Achieving minimum length scale in topology optimizationusing nodal design variables and projection functions. International Journal for Numerical Methods inEngineering , (2):238–254, 2004.[10] Z. Hashin and S. Shtrikman. A variational approach to the theory of the elastic behaviour of multiphasematerials. Journal of the Mechanics and Physics of Solids , :127–140, 1963.[11] B. S. Kirk, J. W. Peterson, R. H. Stogner, and G. F. Carey. libMesh : a C++ library for paralleladaptive mesh refinement/coarsening simulations. Engineering with Computers , (3):237–254, 2006.[12] J. Mackerle. Topology and shape optimization of structures using FEM and BEM: A bibliography(1999–2001). Finite Elements in Analysis and Design , (3):243–253, 2003.[13] M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti. Recycling Krylov subspaces forsequences of linear systems. SIAM Journal on Scientific Computing , (5):1651–1674, 2006.2214] G. H. Paulino and C. H. Le. A modified Q4/Q4 element for topology optimization. Structural andMultidisciplinary Optimization , 2008. In press.[15] G. H. Paulino, R. C. Page III, and E. C. N. Silva. A Java-based topology optimization program with webaccess: nodal design variable approach. In
Proceedings of the 6th World Congress on Structural and Mul-tidisciplinary Optimization (WCSMO6), Rio de Janeiro, Brazil, May 30 – June 3, 2005 . InternationalSociety for Structural and Multidisciplinary Optimization, 2005.[16] S. F. Rahmatalla and C. C. Swan. A Q4/Q4 continuum structural topology optimization implementa-tion.
Structural and Multidisciplinary Optimization , :130–135, 2004. DOI 10.1007/s00158-003-0365-9.[17] G. I. N. Rozvany. Aims, scope, methods, history and unified terminology of computer-aided topol-ogy optimization in structural mechanics. Journal for Structural and Multidisciplinary Optimization , (2):90–108, 2001.[18] O. Sigmund. Design of Material Structures Using Topology Optimization . Ph.D. thesis, Department ofSolid Mechanics, Technical University of Denmark, 1994.[19] O. Sigmund. On the design of compliant mechanisms using topology optimization.
Mechanics ofStructures and Machines , (4):495–526, 1997.[20] O. Sigmund. Topology optimization: a tool for the tailoring of structures and materials. A special issueof the Philosophical Transactions of the Royal Society: Science into the next Millennium (Issue III,Mathematics, Physics and Engineering) , (1765):211–228, 2000.[21] O. Sigmund and J. Petersson. Numerical instabilities in topology optimization: a survey on proce-dures dealing with checkerboards, mesh-dependencies and local minima. Journal for Structural andMultidisciplinary Optimization , (1):68–75, 1998.[22] R. Stainko. An adaptive multilevel approach to the minimal compliance problem in topology optimiza-tion. Communications in Numerical Methods in Engineering , (2):109–118, 2006.[23] S. Wang. Krylov Subspace Methods for Topology Optimization and Adaptive Meshes . Ph.D. thesis,University of Illinois at Urbana-Champaign, Department of Computer Science, Urbana, Illinois 61801,September 2007.[24] S. Wang, E. de Sturler, and G. H. Paulino. Large-scale topology optimization using preconditionedKrylov subspace methods with recycling.
International Journal for Numerical Methods in Engineering ,69