top88 line code

Struct Multidisc Optim (2011) 43:1–16 DOI 10.1007/s00158-010-0594-7 EDUCATIONAL ARTICLE Efficient topology optimization in MATLAB using 88 lines of code Erik Andreassen · Anders Clausen · Mattias Schevenels · Boyan S. Lazarov · Ole Sigmund Received: 19 September 2010 / Revised: 27 October 2010 / Accepted: 28 October 2010 / Published online: 20 November 2010 c© Springer-Verlag 2010 Abstract The paper presents an efficient 88 line MATLAB code for topology optimization. It has been developed using the 99 line code presented by Sigmund (Struct Multidisc Optim 21(2):120–127, 2001) as a start- ing point. The original code has been extended by a density filter, and a considerable improvement in efficiency has been achieved, mainly by preallocating arrays and vector- izing loops. A speed improvement with a factor of 100 is obtained for a benchmark example with 7,500 elements. Moreover, the length of the code has been reduced to a mere 88 lines. These improvements have been accomplished without sacrificing the readability of the code. The 88 line code can therefore be considered as a valuable successor to the 99 line code, providing a practical instrument that may help to ease the learning curve for those entering the field of topology optimization. The paper also discusses simple extensions of the basic code to include recent PDE-based and black-and-white projection filtering methods. The com- plete 88 line code is included as an appendix and can be downloaded from the web site www.topopt.dtu.dk. Keywords Topology optimization · MATLAB · Education · Computational efficiency E. Andreassen · A. Clausen · B. S. Lazarov · O. Sigmund (B) Department of Mechanical Engineering, Solid Mechanics, Technical University of Denmark, Nils Koppels Alle, B. 404, 2800, Lyngby, Denmark e-mail: sigmund@mek.dtu.dk M. Schevenels Department of Civil Engineering, K.U. Leuven, Kasteelpark Arenberg 40, 3001 Leuven, Belgium 1 Introduction MATLAB is a high-level programming language that allows for the solution of numerous scientific problems with a min- imum of coding effort. An example is Sigmund’s 99 line topology optimization code (Sigmund 2001). The 99 line code is intended for educational purposes and serves as an introductory example to topology optimization for stu- dents and newcomers to the field. The use of MATLAB, with its accessible syntax, excellent debugging tools, and extensive graphics handling opportunities, allows the user to focus on the physical and mathematical background of the optimization problem without being distracted by tech- nical implementation issues. Other examples of simple MATLAB code used to provide insight in finite element analysis or topology optimization include a finite element code for the solution of elliptic problems with mixed bound- ary conditions on unstructured grids (Alberty et al. 1999), a similar code for problems in linear elasticity (Alberty et al. 2002), a topology optimization code for compliant mecha- nism design and for heat conduction problems (Bendsøe and Sigmund 2003), a code for Pareto-optimal tracing in topol- ogy optimization (Suresh 2010), a discrete level-set topol- ogy optimization code (Challis 2010), and a Scilab code for two-dimensional optimization problems based on the level set method (Allaire 2009). Compared to high performance programming languages such as C++ and Fortran, MATLAB is generally perceived to be far behind when it comes to computational power. This can partly be explained by (1) the fact that many users apply the same programming strategies as in Fortran or C++, such as the extensive use of for and while loops, and (2) the fact that MATLAB is relatively tolerant towards bad programming practices, such as the use of dynami- cally growing variable arrays. In both cases the potential 2 E. Andreassen et al. of MATLAB is far from optimally utilized. Efficient use of MATLAB implies loop vectorization and memory preallo- cation (The MathWorks 2010). Loop vectorization is the use of vector and matrix operations in order to avoid for and while loops. Memory preallocation means that the max- imum amount of memory required for an array is reserved a priori, hence avoiding the costly operation of reallocating memory and moving data as elements are added to the array. Loop vectorization and memory preallocation are used in combination with a number of more advanced performance improving techniques in the MILAMIN code, a MATLAB program capable of solving two-dimensional finite element problems with one million unknowns in one minute on a desktop computer (Dabrowski et al. 2008). In the 99 line topology optimization code, the perfor- mance of several operations (such as the filtering proce- dure and the assembly of the finite element matrices) can be increased dramatically. Partly by properly exploiting the strengths of MATLAB (using loop vectorization and memory preallocation), partly by restructuring the program (moving portions of code out of the optimization loop so that they are only executed once), a substantial increase in efficiency has been achieved: for an example problem with 7,500 elements, the total computation time has been reduced by a factor 100. In addition, the original code has been extended by the inclusion of density filtering, while reducing the length of the code to only 88 lines. The aim of this paper is to present the 88 line code. It should be considered as a successor to the 99 line code, and it is published with the same objective: to provide an edu- cational instrument for newcomers to the field of topology optimization. The main improvements with respect to the original code are the increased speed and the inclusion of a density filter. These are relevant improvements, as the 99 line code has been downloaded by more than 8,000 unique users since 1999 and is still used as a basis for new devel- opments in the field of topology optimization. The density filter is a useful addition as it paves the way for the imple- mentation of more modern filters such as the Heaviside filters proposed by Guest et al. (2004) and Sigmund (2007). The present text is conceived as an extension of the paper by Sigmund (2001). Large parts of the 88 line code are identical to the original 99 line code, and the same nota- tion is adopted. This approach is followed in an attempt to minimize the effort required to upgrade to the new implementation. The paper is organized as follows. The topology opti- mization problem is formulated in Section 2. As in the original paper, the focus is restricted to minimum compli- ance problems with a constraint on the amount of material available. The 88 line code is explained in Section 3. Spe- cial attention is paid to the portions of the code that have changed with respect to the original 99 line code. These two sections constitute the core of the paper. The remain- ing sections have a supplementary character, addressing variants of and extensions of the 88 line code and dis- cussing its performance. Section 4 presents two alternative implementations of the filtering operation. The first alterna- tive is based on the built-in MATLAB convolution operator function conv2. This modification implies a further reduc- tion of the code to 71 lines and leads to a reduction of the memory footprint, but this comes at the expense of the code’s readability for those unfamiliar with the conv2 function. The second alternative is based on the application of a Helmholtz type partial differential equation to the den- sity or sensitivity field (Lazarov and Sigmund 2010). This approach allows for the use of a finite element solver to per- form the filtering operation, which reduces the complexity of the implementation for serial and parallel machines, as well as the computation time for large problems and com- plex geometries. Section 5 shows how to extend the 88 line code to problems involving different boundary conditions, multiple load cases, and passive elements. Furthermore, the inclusion of a Heaviside filter in order to obtain black-and- white solutions is elaborated. In Section 6, the performance of the 88 line code and its variants is examined. The compu- tation time is analyzed for three benchmark examples solved with both the original 99-line code and the new versions of the code. The memory usage of the new code is also briefly discussed. 2 Problem formulation The MBB beam is a classical problem in topology opti- mization. In accordance with the original paper (Sigmund 2001), the MBB beam is used here as an example. The design domain, the boundary conditions, and the external load for the MBB beam are shown in Fig. 1. The aim of the optimization problem is to find the optimal material distri- bution, in terms of minimum compliance, with a constraint on the total amount of material. Fig. 1 The design domain, boundary conditions, and external load for the optimization of a symmetric MBB beam Efficient topology optimization in MATLAB using 88 lines of code 3 2.1 Modified SIMP approach The design domain is discretized by square finite elements and a “density-based approach to topology optimization” is followed (Bendsøe 1989; Zhou and Rozvany 1991); i.e. each element e is assigned a density xe that determines its Young’s modulus Ee: Ee(xe) = Emin + x pe (E0 − Emin), xe ∈ [0, 1] (1) where E0 is the stiffness of the material, Emin is a very small stiffness assigned to void regions in order to prevent the stiffness matrix from becoming singular, and p is a penalization factor (typically p = 3) introduced to ensure black-and-white solutions. Equation (1) corresponds to the modified SIMP approach, which differs from the classical SIMP approach used in the original paper in the occurrence of the term Emin. In the classical SIMP approach, elements with zero stiffness are avoided by imposing a lower limit slightly larger than zero on the densities xe. The modified SIMP approach has a number of advantages (Sigmund 2007), most importantly that it allows for a straightfor- ward implementation of additional filters, as illustrated in Section 5. The mathematical formulation of the optimization prob- lem reads as follows: min x : c(x) = UTKU = N∑ e=1 Ee(xe)uTe k0ue subject to : V (x)/V0 = f KU = F 0 ≤ x ≤ 1 (2) where c is the compliance, U and F are the global dis- placement and force vectors, respectively, K is the global stiffness matrix, ue is the element displacement vector, k0 is the element stiffness matrix for an element with unit Young’s modulus, x is the vector of design variables (i.e. the element densities), N is the number of elements used to discretize the design domain, V (x) and V0 are the material volume and design domain volume, respectively, and f is the prescribed volume fraction. 2.2 Optimality criteria method The optimization problem (2) is solved by means of a standard optimality criteria method. A heuristic updating scheme identical to the scheme used in the original paper is followed: xnewe = ⎧ ⎪⎨ ⎪⎩ max(0, xe − m) if xeBηe ≤ max(0, xe − m) min(1, xe + m) if xeBηe ≥ min(1, xe − m) xeB η e otherwise (3) where m is a positive move limit, η (= 1/2) is a numerical damping coefficient, and Be is obtained from the optimality condition as: Be = − ∂c ∂xe λ ∂V ∂xe (4) where the Lagrangian multiplier λ must be chosen so that the volume constraint is satisfied; the appropriate value can be found by means of a bisection algorithm. The sensitivities of the objective function c and the mate- rial volume V with respect to the element densities xe are given by: ∂c ∂xe = −px p−1e (E0 − Emin)uTe k0u (5) ∂V ∂xe = 1 (6) Equation (6) is based on the assumption that each element has unit volume. 2.3 Filtering In order to ensure existence of solutions to the topology optimization problem and to avoid the formation of checker- board patterns (Díaz and Sigmund 1995; Jog and Haber 1996; Sigmund and Petersson 1998), some restriction on the design must be imposed. A common approach is the appli- cation of a filter to either the sensitivities or the densities. A whole range of filtering methods is thoroughly described by Sigmund (2007). In addition to the sensitivity filter (Sigmund 1994, 1997), which is already implemented in the 99 line code, the new 88 line code also includes density filtering (Bruns and Tortorelli 2001; Bourdin 2001). The sensitivity filter modifies the sensitivities ∂c/∂xe as follows: ∂̂c ∂xe = 1 max(γ, xe) ∑ i∈Ne Hei ∑ i∈Ne Hei xi ∂c ∂xi (7) where Ne is the set of elements i for which the center-to- center distance �(e, i) to element e is smaller than the filter radius rmin and Hei is a weight factor defined as: Hei = max (0, rmin − �(e, i)) (8) The term γ (= 10−3) in (7) is a small positive number intro- duced in order to avoid division by zero. This is a difference as compared to the original paper, where the classical SIMP approach is used. In the classical SIMP approach, the 4 E. Andreassen et al. density variables cannot become zero, and the term γ is not required. The density filter transforms the original densities xe as follows: x˜e = 1∑ i∈Ne Hei ∑ i∈Ne Hei xi (9) In the following, the original densities xe are referred to as the design variables. The filtered densities x˜e are referred to as the physical densities. This terminology is used to stress the fact that the application of a density filter causes the original densities xe to loose their physical meaning. One should therefore always present the filtered density field x˜e rather than the original density field xe as the solution to the optimization problem (Sigmund 2007). In the case where a density filter is applied, the sensitiv- ities of the objective function c and the material volume V with respect to the physical densities x˜e are still given by (5) and (6), provided that the variable xe is replaced with x˜e. The sensitivities with respect to the design variables x j are obtained by means of the chain rule: ∂ψ ∂x j = ∑ e∈N j ∂ψ ∂ x˜e ∂ x˜e ∂x j = ∑ e∈N j 1 ∑ i∈Ne Hei Hje ∂ψ ∂ x˜e (10) where the function ψ represents either the objective func- tion c or the material volume V . 3 MATLAB implementation In this section the 88 line MATLAB code (see Appendix) is explained. The code is called from the MATLAB prompt by means of the following line: top88(nelx,nely,volfrac,penal,rmin,ft) where nelx and nely are the number of elements in the horizontal and vertical direction, respectively, volfrac is the prescribed volume fraction f , penal is the penaliza- tion power p, rmin is the filter radius rmin (divided by the element size), and the additional argument (compared to the 99 line code) ft specifies whether sensitivity filtering (ft = 1) or density filtering (ft = 2) should be used. When sensitivity filtering is chosen, the 88 line code yields practically1 the same results as the 99 line code; e.g. the 1The slight difference which can be observed between the 88-line and the 99-line code is due to the difference in the SIMP formulation. optimized MBB beam shown in Fig. 1 of the original paper by Sigmund (2001) can be reproduced by means of the following function call: top88(60,20,0.5,3,1.5,1) The most obvious differences between the 88 line code and the 99 line code are the following: (1) the for loops used to assemble the finite element matrices, to compute the compliance, and to perform the filtering operation have been vectorized, (2) the remaining arrays constructed by means of a for loop are properly preallocated, (3) a maxi- mum amount of code is moved out of the optimization loop to ensure that it is only executed once, (4) a distinction is made between the design variables x and the physical densi- ties xPhys in order to facilitate the application of a density filter, and (5) all subroutines have been integrated in the main program. The 88 line code consists of three parts: the finite ele- ment analysis, the sensitivity or density filter, and the optimization loop. These parts are discussed in detail in Sections 3.1–3.3. Section 3.4 presents some results obtained with the 88 line code. 3.1 Finite element analysis The design domain is assumed to be rectangular and dis- cretized with square elements. A coarse example mesh consisting of 12 elements with four nodes per element and two degrees of freedom (DOFs) per node is presented in Fig. 2. Both nodes and elements are numbered column-wise from left to right, and the DOFs 2n − 1 and 2n corre- spond to the horizontal and vertical displacement of node n, respectively. This highly regular mesh can be exploited in several ways in order to reduce the computational effort in the optimization loop to a minimum. The finite element preprocessing part starts with the definition of the material properties (lines 4–6): E0 is the Young’s modulus E0 of the material, Emin is the Fig. 2 The design domain with 12 elements Efficient topology optimization in MATLAB using 88 lines of code 5 artificial Young’s modulus Emin assigned to void regions (or the Young’s modulus of the second material in a two-phase design problem), and nu is the Poisson’s ratio ν. Next the element stiffness matrix k0 for an element with unit Young’s modulus is computed (lines 8–12). This matrix is denoted as KE. Due to the regularity of the mesh, this matrix is identical for all elements. In order to allow for an efficient assembly of the stiffness matrix in the optimization loop, a matrix edofMat is con- structed (lines 13–15). The i-th row of this matrix contains the eight DOF indices corresponding to the i-th element (in a similar way as the edof vector in the original 99 line code). The matrix edofMat is constructed in three steps. First, a (nely+ 1) × (nelx+ 1) matrix nodenrs with the node numbers is defined. The MATLAB func- tion reshape is used; this function returns a matrix with the size specified by the second and third input argument, whose elements are taken column-wise from the first input argument (which is in this case a vector containing the node numbers). Next, the matrix nodenrs is used to determine the first DOF index for all elements, which are stored in a vector edofVec. Finally, the matrix edofVec is used to determine the eight DOF indices for each element. To this end, the MATLAB function repmat is called twice. This function copies a matrix the specified number of times in the vertical and horizontal direction. The first call to the repmat function returns a matrix with eight columns which are all copies of the vector edofVec. The second call returns a matrix of the same size where all rows are identical; this matrix relates the indices of the eight DOFs of an element to the index of its first DOF stored in the vec- tor edofVec. The results are added up and collected in the matrix edofMat. For the example mesh shown in Fig. 2, this procedure yields the following result: edofMat = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎣ 3 4 11 12 9 10 1 2 5 6 13 14 11 12 3 4 7 8 15 16 13 14 5 6 11 12 19 20 17 18 9 10 ... ... ... ... 31 32 39 40 37 38 29 30 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎦ ← Element 1 ← Element 2 ← Element 3 ← Element 4 ← Element 12 In each iteration of the optimization loop, the assembly of the global stiffness matrix K is efficiently performed by means of the sparse function in MATLAB, so avoid- ing the use of for loops. The procedure followed here is inspired by the approach described by Davis (2007). The sparse function takes three vectors as input arguments: the first and second contain the row and column indices of the non-zero matrix entries, which are collected in the third vector. Specifying the same row and column indices multiple times results in a summation of the corresponding entries. The row and colums index vectors (iK and jK, respec- tively) are created in lines 16 and 17 using the edofMat matrix. Use is made of a Kronecker matrix product with a unit vector of length 8, followed by a reshaping oper- ation. The resulting vectors iK and jK are structured so tha