Structural and Multidisciplinary Optimization manuscript No.

(will be inserted by the editor)

Efficient topology optimization in MATLAB using 88 lines of code

Erik Andreassen · Anders Clausen · Mattias Schevenels · Boyan S. Lazarov · Ole
Sigmund

Received: date / Accepted: date

Abstract The paper presents an efficient 88 line MATLAB line topology optimization code (Sigmund, 2001). The 99
code for topology optimization. It has been developed using line code is intended for educational purposes and serves
the 99 line code presented by Sigmund (2001) as a starting as an introductory example to topology optimization for
point. The original code has been extended by a density students and newcomers to the field. The use of MATLAB,
filter, and a considerable improvement in efficiency has been with its accessible syntax, excellent debugging tools, and
achieved, mainly by preallocating arrays and vectorizing extensive graphics handling opportunities, allows the user
loops. A speed improvement with a factor of 100 is obtained to focus on the physical and mathematical background
for a benchmark example with 7500 elements. Moreover, of the optimization problem without being distracted by
the length of the code has been reduced to a mere 88 technical implementation issues. Other examples of simple
lines. These improvements have been accomplished without MATLAB code used to provide insight in finite element
sacrificing the readability of the code. The 88 line code can analysis or topology optimization include a finite element
therefore be considered as a valuable successor to the 99 line code for the solution of elliptic problems with mixed
code, providing a practical instrument that may help to ease boundary conditions on unstructured grids (Alberty et al,
the learning curve for those entering the field of topology 1999), a similar code for problems in linear elasticity
optimization. The paper also discusses simple extensions of (Alberty et al, 2002), a topology optimization code for
the basic code to include recent PDE-based and black-and- compliant mechanism design and for heat conduction
white projection filtering methods. The complete 88 line problems (Bendsøe and Sigmund, 2003), a code for Pareto-
code is included as an appendix and can be downloaded optimal tracing in topology optimization (Suresh, 2010),
from the web site www.topopt.dtu.dk. a discrete level-set topology optimization code (Challis,
2010), and a Scilab code for two-dimensional optimization
Keywords Topology optimization · MATLAB · Education ·
problems based on the level set method (Allaire, 2009).
Computational efficiency
Compared to high performance programming languages
such as C++ and Fortran, MATLAB is generally perceived
1 Introduction to be far behind when it comes to computational power.
This can partly be explained by (1) the fact that many users
MATLAB is a high-level programming language that allows apply the same programming strategies as in Fortran or
for the solution of numerous scientific problems with a C++, such as the extensive use of for and while loops,
minimum of coding effort. An example is Sigmund’s 99 and (2) the fact that MATLAB is relatively tolerant towards
bad programming practices, such as the use of dynamically
E. Andreassen, A. Clausen, O. Sigmund∗ , B.S. Lazarov growing variable arrays. In both cases the potential of
Department of Mechanical Engineering, Solid Mechanics,
Technical University of Denmark, Nils Koppels Alle, B. 404,
MATLAB is far from optimally utilized. Efficient use
DK-2800 Lyngby, Denmark of MATLAB implies loop vectorization and memory
∗ E-mail: sigmund@mek.dtu.dk
preallocation (The MathWorks, 2010). Loop vectorization is
M. Schevenels the use of vector and matrix operations in order to avoid
Department of Civil Engineering, K.U.Leuven, for and while loops. Memory preallocation means that
Kasteelpark Arenberg 40, B-3001 Leuven, Belgium the maximum amount of memory required for an array is

Preprint submitted to Structural and Multidisciplinary Optimization 27 October 2010

Published version: http://dx.doi.org/10.1007/s00158-010-0594-7
2

reserved a priori, hence avoiding the costly operation of in MATLAB convolution operator function conv2. This
reallocating memory and moving data as elements are added modification implies a further reduction of the code to 71
to the array. Loop vectorization and memory preallocation lines and leads to a reduction of the memory footprint,
are used in combination with a number of more advanced but this comes at the expense of the code’s readability
performance improving techniques in the MILAMIN code, for those unfamiliar with the conv2 function. The second
a MATLAB program capable of solving two-dimensional alternative is based on the application of a Helmholtz type
finite element problems with one million unknowns in one partial differential equation to the density or sensitivity
minute on a desktop computer (Dabrowski et al, 2008). field (Lazarov and Sigmund, 2010). This approach allows
In the 99 line topology optimization code, the perfor- for the use of a finite element solver to perform the
mance of several operations (such as the filtering procedure filtering operation, which reduces the complexity of the
and the assembly of the finite element matrices) can implementation for serial and parallel machines, as well
be increased dramatically. Partly by properly exploiting as the computation time for large problems and complex
the strengths of MATLAB (using loop vectorization and geometries. Section 5 shows how to extend the 88 line
memory preallocation), partly by restructuring the program code to problems involving different boundary conditions,
(moving portions of code out of the optimization loop so multiple load cases, and passive elements. Furthermore, the
that they are only executed once), a substantial increase in inclusion of a Heaviside filter in order to obtain black-and-
efficiency has been achieved: for an example problem with white solutions is elaborated. In section 6, the performance
7500 elements, the total computation time has been reduced of the 88 line code and its variants is examined. The
by a factor 100. In addition, the original code has been computation time is analyzed for three benchmark examples
extended by the inclusion of density filtering, while reducing solved with both the original 99-line code and the new
the length of the code to only 88 lines. versions of the code. The memory usage of the new code
The aim of this paper is to present the 88 line code. It is also briefly discussed.
should be considered as a successor to the 99 line code,
and it is published with the same objective: to provide
an educational instrument for newcomers to the field of 2 Problem formulation
topology optimization. The main improvements with respect
The MBB beam is a classical problem in topology opti-
to the original code are the increased speed and the inclusion
mization. In accordance with the original paper (Sigmund,
of a density filter. These are relevant improvements, as the
2001), the MBB beam is used here as an example. The
99 line code has been downloaded by a more than 8000
design domain, the boundary conditions, and the external
people since 1999 and is still used as a basis for new
load for the MBB beam are shown in figure 1. The aim
developments in the field of topology optimization. The
of the optimization problem is to find the optimal material
density filter is a useful addition as it paves the way for the
distribution, in terms of minimum compliance, with a
implementation of more modern filters such as the Heaviside
constraint on the total amount of material.
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
notation 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
optimization problem is formulated in section 2. As in
the original paper, the focus is restricted to minimum
compliance problems with a constraint on the amount Fig. 1 The design domain, boundary conditions, and external load for
the optimization of a symmetric MBB beam.
of material available. The 88 line code is explained in
section 3. Special 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 remaining sections have a supplementary 2.1 Modified SIMP approach
character, addressing variants of and extensions of the
88 line code and discussing its performance. Section 4 The design domain is discretized by square finite elements
presents two alternative implementations of the filtering and a “density-based approach to topology optimization”
operation. The first alternative is based on the built- is followed (Bendsøe, 1989; Zhou and Rozvany, 1991); i.e.
 3

each element e is assigned a density xe that determines its condition as:

Young’s modulus Ee :
∂c
−
Ee (xe ) = Emin + xpe (E0 − Emin ), xe ∈ [0, 1] (1) ∂xe
Be = (4)
∂V
λ
where E0 is the stiffness of the material, Emin is a ∂xe
very small stiffness assigned to void regions in order to where the Lagrangian multiplier λ must be chosen so that
prevent the stiffness matrix from becoming singular, and the volume constraint is satisfied; the appropriate value can
p is a penalization factor (typically p = 3) introduced to be found by means of a bisection algorithm.
ensure black-and-white solutions. Equation (1) corresponds The sensitivities of the objective function c and the
to the modified SIMP approach, which differs from the material volume V with respect to the element densities xe
classical SIMP approach used in the original paper in are given by:
the occurrence of the term Emin . In the classical SIMP
approach, elements with zero stiffness are avoided by ∂c
= −pxp−1
e (E0 − Emin )uT
e k0 u (5)
imposing a lower limit slightly larger than zero on the ∂xe
densities xe . The modified SIMP approach has a number of ∂V
=1 (6)
advantages (Sigmund, 2007), most importantly that it allows ∂xe
for a straightforward implementation of additional filters, as
illustrated in section 5. Equation (6) is based on the assumption that each element
The mathematical formulation of the optimization has unit volume.
problem reads as follows:
N
X 2.3 Filtering
min: c(x) = UT KU = Ee (xe )uT
e k0 ue
x
e=1 In order to ensure existence of solutions to the topology
subject to: V (x)/V0 = f (2)
optimization problem and to avoid the formation of checker-
KU = F board patterns (Dı́az and Sigmund, 1995; Jog and Haber,
0≤x≤1 1996; Sigmund and Petersson, 1998), some restriction on
the design must be imposed. A common approach is
where c is the compliance, U and F are the global
the application of a filter to either the sensitivities or
displacement and force vectors, respectively, K is the global
the densities. A whole range of filtering methods is
stiffness matrix, ue is the element displacement vector, k0
thoroughly described by Sigmund (2007). In addition to the
is the element stiffness matrix for an element with unit
sensitivity filter (Sigmund, 1994, 1997), which is already
Young’s modulus, x is the vector of design variables (i.e.
implemented in the 99 line code, the new 88 line code
the element densities), N is the number of elements used to
also includes density filtering (Bruns and Tortorelli, 2001;
discretize the design domain, V (x) and V0 are the material
Bourdin, 2001).
volume and design domain volume, respectively, and f is
The sensitivity filter modifies the sensitivities ∂c/∂xe as
the prescribed volume fraction.
follows:
d
∂c 1 X ∂c
= X Hei xi (7)
2.2 Optimality criteria method ∂xe max(γ, xe ) Hei i∈Ne ∂xi
i∈Ne
The optimization problem (2) is solved by means of a
standard optimality criteria method. A heuristic updating where Ne is the set of elements i for which the center-to-
scheme identical to the scheme used in the original paper center distance ∆(e, i) to element e is smaller than the filter
is followed: radius rmin and Hei is a weight factor defined as:

 η Hei = max (0, rmin − ∆(e, i)) (8)
max(0, xe − m) if xe Be ≤ max(0, xe − m)
xnew
e = min(1, xe + m) if xe Beη ≥ min(1, xe − m) The term γ (= 10−3 ) in equation (7) is a small positive


xe Beη otherwise number introduced in order to avoid division by zero. This
(3) is a difference as compared to the original paper, where
the classical SIMP approach is used. In the classical SIMP
where m is a positive move limit, η (= 1/2) is a numerical approach, the density variables cannot become zero, and the
damping coefficient, and Be is obtained from the optimality term γ is not required.
4

The density filter transforms the original densities xe as used to assemble the finite element matrices, to compute the
follows: compliance, and to perform the filtering operation have been
1 X vectorized, (2) the remaining arrays constructed by means
x̃e = X Hei xi (9) of a for loop are properly preallocated, (3) a maximum
Hei i∈Ne amount of code is moved out of the optimization loop to
i∈Ne
ensure that it is only executed once, (4) a distinction is made
In the following, the original densities xe are referred to as between the design variables x and the physical densities
the design variables. The filtered densities x̃e are referred to xPhys in order to facilitate the application of a density filter,
as the physical densities. This terminology is used to stress and (5) all subroutines have been integrated in the main
the fact that the application of a density filter causes the program.
original densities xe to loose their physical meaning. One The 88 line code consists of three parts: the finite
should therefore always present the filtered density field x̃e element analysis, the sensitivity or density filter, and the
rather than the original density field xe as the solution to the optimization loop. These parts are discussed in detail in
optimization problem (Sigmund, 2007). subsections 3.1 to 3.3. Subsection 3.4 presents some results
In the case where a density filter is applied, the obtained with the 88 line code.
sensitivities of the objective function c and the material
volume V with respect to the physical densities x̃e are still
given by equations (5) and (6), provided that the variable 3.1 Finite element analysis
xe is replaced with x̃e . The sensitivities with respect to the
design variables xj are obtained by means of the chain rule: The design domain is assumed to be rectangular and
X ∂ψ ∂ x̃e X discretized with square elements. A coarse example mesh
∂ψ 1 ∂ψ
= = X Hje (10) consisting of 12 elements with four nodes per element and
∂xj ∂ x̃e ∂xj Hei ∂ x̃e two degrees of freedom (DOFs) per node is presented in
e∈Nj e∈Nj
i∈Ne figure 2. Both nodes and elements are numbered column-
where the function ψ represents either the objective function wise from left to right, and the DOFs 2n − 1 and 2n
c or the material volume V . correspond 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
3 MATLAB implementation computational effort in the optimization loop to a minimum.

In this section the 88 line MATLAB code (see appendix) is

explained. The code is called from the MATLAB prompt by 1,2 9,10 17,18 25,26 33,34
means of the following line:
1 4 7 10
top88(nelx,nely,volfrac,penal,rmin,ft) 3,4 11,12 19,20 27,28 35,36

where nelx and nely are the number of elements in the 2 5 8 11

5,6 13,14 21,22 29,30 37,38
horizontal and vertical direction, respectively, volfrac is
the prescribed volume fraction f , penal is the penalization 3 6 9 12
power p, rmin is the filter radius rmin (divided by the 7,8 15,16 23,24 31,32 39,40
element size), and the additional argument (compared to
the 99 line code) ft specifies whether sensitivity filtering
Fig. 2 The design domain with 12 elements.
(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 The finite element preprocessing part starts with the
optimized MBB beam shown in figure 1 of the original definition of the material properties (lines 4-6): E0 is the
paper by Sigmund (2001) can be reproduced by means of Young’s modulus E0 of the material, Emin is the artificial
the following function call: Young’s modulus Emin assigned to void regions (or the
top88(60,20,0.5,3,1.5,1) Young’s modulus of the second material in a two-phase
design problem), and nu is the Poisson’s ratio ν.
The most obvious differences between the 88 line code Next the element stiffness matrix k0 for an element with
and the 99 line code are the following: (1) the for loops unit Young’s modulus is computed (lines 8-12). This matrix
1 The slight difference which can be observed between the 88-line is denoted as KE. Due to the regularity of the mesh, this
and the 99-line code is due to the difference in the SIMP formulation. matrix is identical for all elements.
 5

In order to allow for an efficient assembly of the The third vector, containing the entries of the sparse
stiffness matrix in the optimization loop, a matrix edofMat stiffness matrix, is computed in the optimization loop (line
is constructed (lines 13-15). The i-th row of this matrix 54), as it depends on the physical densities x̃. This vector
contains the eight DOF indices corresponding to the i- sK is obtained by reshaping the element stiffness matrix KE
th element (in a similar way as the edof vector in the to obtain a column vector, multiplying this vector with the
original 99 line code). The matrix edofMat is constructed appropriate Young’s modulus Ee (x̃e ) for each element, and
in three steps. First, a (nely + 1) × (nelx + 1) matrix concatenating the results for all elements. The multiplication
nodenrs with the node numbers is defined. The MATLAB and concatenation are implemented as a matrix product
function reshape is used; this function returns a matrix with followed by a reshaping operation.
the size specified by the second and third input argument, The actual assembly of the stiffness matrix K is
whose elements are taken column-wise from the first input performed on line 55 by means of the sparse function,
argument (which is in this case a vector containing the node using the index vectors iK and jK and the vector with non-
numbers). Next, the matrix nodenrs is used to determine zero entries sK. This procedure could be further improved
the first DOF index for all elements, which are stored in by using the sparse2 function from CHOLMOD (Davis,
a vector edofVec. Finally, the matrix edofVec is used to 2008), which is faster than the standard MATLAB sparse
determine the eight DOF indices for each element. To this function due to the use of a more efficient sorting algorithm
end, the MATLAB function repmat is called twice. This for the indices, but this is beyond the scope of the present
function copies a matrix the specified number of times in paper. The second statement on line 55 ensures that the
the vertical and horizontal direction. The first call to the stiffness matrix is perfectly symmetric. This is important as
repmat function returns a matrix with eight columns which it determines the algorithm used by MATLAB to solve the
are all copies of the vector edofVec. The second call returns system of finite element equations. If the stiffness matrix is
a matrix of the same size where all rows are identical; this sparse, symmetric, and has real positive diagonal elements,
matrix relates the indices of the eight DOFs of an element to Cholesky factorization is used. If the stiffness matrix is
the index of its first DOF stored in the vector edofVec. The not symmetric (due to rounding errors in the assembly
results are added up and collected in the matrix edofMat. procedure), LU factorization is used instead, resulting in a
For the example mesh shown in figure 2, this procedure longer computation time.
yields the following result: The boundary conditions and the load vector are defined
  on lines 18-23. These lines are almost identical to those in
3 4 11 12 9 10 1 2 ← Element 1 the original 99 line code and are therefore not discussed in
 5 6 13 14 11 12 3 4  ← Element 2 the present paper. The main difference with the original code
 
 7 8 15 16 13 14 5 6  ← Element 3 is that these lines are moved out of the optimization loop.
 
edofMat =  11 12 19 20 17 18 9 10  ← Element 4
  The system of finite element equations is finally solved
 . . . . 
 .. .. .. ..  on line 56.
31 32 39 40 37 38 29 30 ← Element 12

In each iteration of the optimization loop, the assembly 3.2 Filtering

of the global stiffness matrix K is efficiently performed by
means of the sparse function in MATLAB, so avoiding the The application of a sensitivity filter according to
use of for loops. The procedure followed here is inspired equation (7) involves a weighted average over different
by the approach described by Davis (2007). The sparse elements. This is a linear operation; it can therefore be
function takes three vectors as input arguments: the first implemented as a matrix product of a coefficient matrix
and second contain the row and column indices of the non- and a vector containing the original sensitivities ∂c/∂xi
zero matrix entries, which are collected in the third vector. (multiplied with the design P variables xi ). Dividing the
Specifying the same row and column indices multiple times result by a factor max(γ, xe ) i∈Ne Hei yields the filtered
results in a summation of the corresponding entries. d e . This operation is performed on line 64.
sensitivities ∂c/∂x
The row and colums index vectors (iK and jK, The matrix H and the vector Hs contain
P the coefficients Hei
respectively) are created in lines 16-17 using the edofMat and the normalization constants i∈Ne Hei , respectively.
matrix. Use is made of a Kronecker matrix product with a The use of a density filter not only implies filtering of
unit vector of length 8, followed by a reshaping operation. the densities according to equation (9) but also a chain rule
The resulting vectors iK and jK are structured so that the modification of the sensitivities of the objective function
indices iK(k) and jK(k) correspond to the (i, j)-th entry of and the volume constraint according to equation (10).
the stiffness matrix for element e, where k = i + 8(j − 1) + Both operations involve a weighted average over different
64(e − 1). elements. The density filtering is performed on line 77, the
6

modification of the sensitivities on lines 66-67. Use is made The sensitivities are subsequently filtered (if sensitivity
of the same coefficients H and normalization constants Hs as filtering is used) or modified (if density filtering is used) as
described above. explained in subsection 3.2 (lines 63-68).
Both the matrix H and the vector Hs remain invariant On lines 70-82, an optimality criteria method is used to
during the optimization and are computed a priori. The update the design variables according to equation (3). The
(nelx × nely) × (nelx × nely) coefficient matrix H update is performed in a similar way as in the original 99
establishes a relationship between all elements. However, line code, except that (1) the sensitivity dv of the volume
as the filter kernel defined in equation (8) has a bounded constraint is explicitly taken into account, (2) the Lagrange
support, only neighboring elements affect one another. As a multiplier lmid is determined using the physical densities
consequence, the majority of the coefficients is zero and the instead of the design variables, and (3) the stop condition
matrix H is sparse. It is constructed by means of the built-in is specified in relative terms. The first change is made for
sparse MATLAB function. Row and column index vectors the sake of the density filter: in the sensitivity filtering
iH and jH as well as a vector sH with non-zero entries are approach, the sensitivities dv are identical for all elements
assembled by means of four nested for loops on lines 25- and can therefore be omitted from the definition of the
42. In order to avoid continuous resizing of these vectors as heuristic updating factor Be , but in the density filtering
entries are added, a sufficient (but slightly too high) amount approach, this is no longer true due to the modification of
of memory is preallocated. The entries that remain unused the sensitivities performed on line 67. The second change is
in the vectors iH, jH, and sH have no effect: they preserve strictly speaking not necessary: the density filter is volume-
their initial value (1, 1, and 0, respectively) and result in preserving, which means that the volume constraint can
the addition of a zero term to the first element of the sparse equally well be evaluated in terms of the design variables.
matrix H. The assembly of the matrix H from the vectors When another (non-volume-preserving) filter is applied,
iH, jH, and sH is performed on line 43. The vector Hs is however, it is absolutely necessary to evaluate the volume
subsequently computed on line 44. constraint in terms of the physical densities. The third
change is simply made to optimize the balance between
3.3 Optimization loop accuracy and computation speed.
Finally, the intermediate results are printed (lines 84-85)
The main part of the 88 line code is the optimization and plotted (line 87) in the same way as in the original 99
loop. The loop is initialized on lines 46-49. All design line code.
variables xe are initially set equal to the prescribed volume
The optimization loop is terminated when the L∞ norm
fraction f . The corresponding physical densities x̃e are
of the difference between two consecutive designs (in terms
identical to the design variables xe : in the sensitivity
of design variables) is less than 1 percent.
filtering approach, this equality always holds, while in the
density filtering approach, it holds as long as the design
variables represent a homogeneous field. For other types
of filters (especially non-volume-preserving filters), it may
be necessary to compute the initial physical densities x̃e by
explicit application of the filter to the initial design variables
3.4 Results
xe , and to adjust the initial design variables in such a way
that the volume constraint is satisfied (as this constraint is
specified in terms of the physical densities x̃e ). The 88 line code is used to optimize the MBB beam. Three
Each iteration of the optimization loop starts with the different mesh sizes are considered, consisting of 60 × 20
finite element analysis as described in subsection 3.1 (lines elements, 150 × 50 elements, and 300 × 100 elements. The
54-56). volume constraint is set to 50 % and the usual value p = 3
Next, the objective function (the compliance) c is is used for the penalization exponent. The problem is solved
computed, as well as the sensitivities dc and dv of using sensitivity and density filtering. The filter radius rmin
the objective function and the volume constraint with equals 0.04 times the width of the design domain, i.e. 2.4, 6,
respect to the physical densities (lines 58-61). Compared and 16 for the different meshes.
to the original 99 line code, efficient use is made of the Figure 3 shows the optimized design and the cor-
edofMat matrix to compute the compliance for all elements responding compliance c. The figures demonstrate that
simultaneously: the edofMat matrix is used as an index into both sensitivity filtering and density filtering suppress
the displacement vector U, resulting in a matrix with the size checkerboard patterns and lead to mesh independent
of edofMat that contains the displacements corresponding designs; refining the mesh only leads to a refinement of the
to the DOFs listed in edofMat. solution, not to a different topology.
 7

c = 216.81 c = 219.52 c = 222.29

c = 233.71 c = 235.73 c = 238.31

Fig. 3 Optimized design of the MBB beam and corresponding compliance c obtained with the 88 line code using sensitivity filtering (top) and
density filtering (bottom). A mesh with 60 × 20 elements (left), 150 × 50 elements (middle), and 300 × 100 elements (right) has been used.

4 Alternative implementations The density filter, defined in equation (9), is reformu-

lated as follows:
This section presents two alternatives to the 88 line code X
discussed in the previous section. The focus is on the H(m, n)x(k−m,l−n)
implementation of the filters. m,n
x̃(k,l) = X (11)
The first alternative makes use of the built-in MATLAB H(m, n)
function conv2. This approach is mathematically equivalent m,n
to the implementation presented in the previous section, and
where x(i,j) and x̃(i,j) denote the design variable and the
it allows for a reduction of the code to 71 lines. It also
physical density, respectively, for the element in the i-th row
leads to a reduction of the memory requirements, as will
and the j-th column. The filter kernel H(m, n) is a function
be discussed in subsection 6.2. A possible disadvantage of
of the discrete variables m and n:
this approach is that it may obfuscate the filtering procedure
for readers unfamiliar with the conv2 function and that its H(m, n) = max (0, rmin − δ(m, n)) (12)
applicability is limited to regular meshes.
The second alternative presents the use of filtering where δ(m, n) represents the center-to-center distance
based on a Helmholtz type partial differential equation between two elements separated by m rows and n columns.
(PDE). This approach allows for the use of a finite element Both sums in equation (11) must be taken over all indices
solver to perform the filtering operation, which speeds m and n for which the kernel H(m, n) is non-zero and for
up significantly the filtering process for three-dimensional which (k − m, l − n) refers to an existing element.
problems and simplifies parallel implementations of filtered The non-zero part of the filter kernel H(m, n) can be
topology optimization problems. The results obtained with expressed as an M × N matrix h defined as:
the PDE filter are similar to those obtained using an
exponentially decaying filter kernel (Bruns and Tortorelli, h(m+ M +1 ,n+ N +1 ) = H(m, n) (13)
2 2
2001).
Introducing equation (13) in equation (11) yields the
following expression:
4.1 Filtering using the CONV2 function X
h(m+ M +1 ,n+ N +1 ) x(k−m,l−n)
2 2
The optimization problem discussed in the previous m,n
sections has two properties that allow for a more concise x̃(k,l) = X (14)
h(m+ M +1 ,n+ N +1 )
implementation. First, a rectangular mesh consisting of m,n
2 2

rectangular (square) finite elements is used. Second, the

filter kernel is invariant in space (or, loosely speaking, the The sum in the numerator corresponds to the (k, l)-th ele-
filter radius rmin is the same at all positions in the design ment of the central part of the two-dimensional convolution
domain). As a consequence, the filtering operation can be of the matrices x and h, which is obtained in MATLAB as
interpreted as a two-dimensional discrete convolution. In conv2(x,h,’same’). The sum in the denominator must
the following paragraphs, the convolution based approach is be taken over the same indices, which is most easily
elaborated for the filtering of the densities. The filtering or accomplished by using the same MATLAB code as for the
modification of the sensitivities can be addressed in a similar numerator, substituting the matrix x with a unit matrix of
way. the same size.
8

∂c
Using the convolution based approach for the density field in (15) replaced by ψ = x ∂x and the output field given
filter, the modification of the sensitivities, and the sensitivity ˜
∂c
by ψ̃ = x (Lazarov and Sigmund, 2009).
∂x
filter allows for a reduction of the 88 line code to 71 lines. The filter properties have been discussed extensively by
Three modifications are required. Lazarov and Sigmund (2010), and here only the main ad-
First, the preparation of the filter (lines 25-44) is vantages with respect to memory usage and computational
replaced with the following lines: cost are highlighted. The classical filter requires information
[dy,dx] = meshgrid(-ceil(rmin)+1:ceil(rmin)-1, ... about the neighbor elements, which for irregular meshes and
-ceil(rmin)+1:ceil(rmin)-1); complex geometries is obtained by a relatively expensive
h = max(0,rmin-sqrt(dx.^2+dy.^2));
search. Clearly the approach presented in subsection 3.2
Hs = conv2(ones(nely,nelx),h,’same’);
speeds up the filtering process if the search procedure is
where the matrix h is the non-zero part of the filter kernel performed only once as a preprocessing step, however, the
and the matrix Hs represents the sum in the denominator computational complexity and the memory utilization are
on the right hand side of equation (14). This sum does 2 3
proportional to rmin in two dimensions and to rmin in three
not change during the optimization loop and is therefore dimensions, respectively. The PDE filter approach utilizes
computed in advance. Note that the matrix Hs computed the mesh used for the state problem and does not require
here is identical to the matrix Hs in the 88 line code. any additional information, which avoids excessive memory
Second, the filtering or modification of the sensitivities usage. Furthermore, the computational cost depends linearly
(lines 63-68) is replaced with the following code: on the length parameter Rmin if the solution of the PDE
if ft == 1 (15) is obtained by an iterative method. Therefore, for
dc = conv2(dc.*xPhys,h,’same’)./Hs./max(1e-3,xPhys); a large filter radius, especially in three dimensions, the
elseif ft == 2 PDE filtering scheme should be the preferred choice. In
dc = conv2(dc./Hs,h,’same’);
the presented two-dimensional examples with a regular
dv = conv2(dv./Hs,h,’same’);
end mesh, the concept will not result in improved performance,
however, we include it here for educational reasons and
Third, the filtering of the densities (line 77) is replaced inspiration.
with the following line: FE discretization of equation (15) leads to the following
xPhys = conv2(xnew,h,’same’)./Hs; system of linear equations:
KF x̃N = TF x (18)

4.2 Filtering based on Helmholtz type differential equations where KF is the standard FE stiffness matrix for scalar
problems, TF is a matrix which maps the element design
The density filter given by (9) can be implicitly rep- values x to a vector with nodal values, and x̃N is the
resented by the solution of a Helmholtz type PDE nodal representation of the filtered field. The element-wise
(Lazarov and Sigmund, 2010) with homogeneous Neumann representation of the filtered field is obtained as:
boundary conditions: x̃ = TT
F x̃N (19)
2
−Rmin ∇2 ψ̃ + ψ̃ = ψ (15) The PDE filter requires minor changes of the 88 line
∂ ψ̃ code and reduces it to 82 lines. The preparation of the filter
=0 (16) (lines 25-44) is replaced with the following lines:
∂n
Rmin = rmin/2/sqrt(3);
where ψ is a continuous representation of the unfiltered
KEF = Rmin^2*[4 -1 -2 -1; -1 4 -1 -2; ...
design field, and ψ̃ is the filtered field. The solution -2 -1 4 -1; -1 -2 -1 4]/6 + ...
of the above PDE can be written in the form of a [4 2 1 2; 2 4 2 1; ...
convolution integral which is equivalent to the classical 1 2 4 2; 2 1 2 4]/36;
edofVecF = reshape(nodenrs(1:end-1,1:end-1),nelx*nely,1);
filter. The parameter Rmin in (15) plays a similar role as
edofMatF = repmat(edofVecF,1,4) + ...
rmin in (8). An approximate relation between the length repmat([0 nely+[1:2] 1],nelx*nely,1);
scales for the classical and the PDE filter is given by iKF = reshape(kron(edofMatF,ones(4,1))’,16*nelx*nely,1);
(Lazarov and Sigmund, 2010): jKF = reshape(kron(edofMatF,ones(1,4))’,16*nelx*nely,1);
sKF = reshape(KEF(:)*ones(1,nelx*nely),16*nelx*nely,1);
√
Rmin = rmin /2 3 (17) KF = sparse(iKF,jKF,sKF);
LF = chol(KF,’lower’);
iTF = reshape(edofMatF,4*nelx*nely,1);
The PDE filter is volume preserving, i.e. the volume of the
jTF = reshape(repmat([1:nelx*nely],4,1)’,4*nelx*nely,1);
input field is equal to the volume of the filtered field. The sTF = repmat(1/4,4*nelx*nely,1);
same idea can be applied as a sensitivity filter with the input TF = sparse(iTF,jTF,sTF);
 9

c = 218.79 c = 217.88 c = 219.44

c = 237.60 c = 235.36 c = 236.62

Fig. 4 Optimized design of the MBB beam and corresponding compliance c obtained with the variant of the 88 line code using PDE based
sensitivity filtering (top) and density filtering (bottom). A mesh with 60 × 20 elements (left), 150 × 50 elements (middle), and 300 × 100 elements
(right) has been used.

where KF corresponds to the tangent filter matrix and TF original paper are reconsidered, now starting from the 88
corresponds to the transformation matrix on the right hand line code. In addition, the implementation of a black-and-
side of equation (18). In order to keep the MATLAB code white projection filter is also addressed.
readable, the linear system obtained by FE discretization
of equation (15) is solved by factorization instead of an
iterative method, which hides some of the filter advantages. 5.1 Other boundary conditions
The second change is a replacement of the filtering or
modification of the sensitivities (lines 63-68) with the
following code:
if ft == 1
dc(:) = (TF’*(LF’\(LF\(TF*(dc(:).*xPhys(:)))))) ...
./max(1e-3,xPhys(:));
elseif ft == 2
dc(:) = TF’*(LF’\(LF\(TF*dc(:))));
dv(:) = TF’*(LF’\(LF\(TF*dv(:))));
end
Fig. 5 The design domain, boundary conditions, and external load for
Finally, the filtering of the densities (line 77) is replaced with the optimization of a cantilever beam (left) and the optimized design
the following line: obtained with a variant of the 88 line code using sensitivity filtering
(right).
xPhys(:) = (TF’*(LF’\(LF\(TF*xnew(:)))));

Figure 4 shows the optimized MBB beam and the

Changing load and support conditions in order to solve other
corresponding compliance c obtained with the PDE filter,
optimization problems is very straightforward. In order to
using the same input parameters as in subsection 3.4. The
solve the short cantilever example shown in figure 5, line 19
figure shows that the PDE filter leads to a mesh independent
of the 88 line code must be changed to:
design without checkerboard patterns. The optimized design
and the corresponding compliance c are similar to those F = sparse(2*(nely+1)*(nelx+1),1,-1, ...
2*(nely+1)*(nelx+1),1);
obtained with the standard density and sensitivity filters
shown in figure 3. They are not identical, however. The Line 21 must be changed to:
difference is due to the fact that the PDE filter is based on fixeddofs = [1:2*nely+1];
an exponentially decaying filter kernel, while the standard
With these changes, the optimized design shown in figure 5
filters are based on a linearly decaying filter kernel.
is obtained by means of the following function call:
top88(160,100,0.4,3,6,1)
5 Extensions

Sigmund (2001) describes how to extend the 99 line code 5.2 Multiple load cases
to account for different boundary conditions, multiple load
cases, and passive elements, and how to replace the optimal- It is also very simple to extend the algorithm to account for
ity criteria based optimizer with a more general optimization multiple load cases. As an example, the problem outlined in
scheme. In this section, the extensions discussed in the figure 6 is considered.
10

with a density fixed to be zero or one. As an example, the

optimization problem defined in figure 7 is addressed. A
circular region of the design domain with radius nely/3 and
center (nely/2, nelx/3) is fixed to be void.

Fig. 6 The design domain, boundary conditions, and external loads

for the optimization of a cantilever beam with two load cases (left and
middle) and the optimized design obtained with a variant of the 88 line
code using sensitivity filtering (right).
Fig. 7 The design domain, boundary conditions, and external load for
the optimization of a cantilever beam with a fixed hole (left) and the
In the case of two load cases, the force and displacement optimized design obtained with a variant of the 88 line code using
vectors must be defined as two-column vectors, which sensitivity filtering (right).
means that lines 19 and 20 are changed to:
F = sparse([2*(nely+1)*nelx+2,2*(nely+1)*(nelx+1)], ...
The load vector (line 19) and the support conditions
[1 2],[1 -1],2*(nely+1)*(nelx+1),2);
U = zeros(2*(nely+1)*(nelx+1),2); (line 21) in this example are defined in the same way as in
subsection 5.1. In order to distinguish between active and
The support conditions (line 21) are defined in the same way passive elements, a nely × nelx matrix passive is defined
as in the previous subsection. The equilibrium equations with 0 at elements free to change, 1 at elements fixed to be
must be solved for both load cases, which is accomplished void, and 2 at elements fixed to be solid:
by changing line 56 as follows:
passive = zeros(nely,nelx);
U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:); for i = 1:nelx
for j = 1:nely
The objective function is now defined as the sum of two if sqrt((j-nely/2)^2+(i-nelx/3)^2) < nely/3
compliances: passive(j,i) = 1;
end
2
X end
c(x) = UT
i KUi (20) end
i=1
These lines must be inserted in the 88 line code before the
Lines 58-60 are thus replaced with the following code: start of the optimization loop. The optimality criteria method
c=0;
must be modified by adding the following code between
dc=0; lines 78 and 79:
for i = 1:size(F,2) xPhys(passive==1) = 0;
Ui = U(:,i); xPhys(passive==2) = 1;
ce = reshape(sum((Ui(edofMat)*KE).*Ui(edofMat),2), ...
nely,nelx); With these modifications, the 88 line code can be used to
c = c + sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce));
generate the optimized design shown in figure 7 by means
dc = dc - penal*(E0-Emin)*xPhys.^(penal-1).*ce;
end of the following function call:
top88(150,100,0.5,3,5,1)
The optimized design shown in figure 6 can now be obtained
by means of the following function call:
top88(150,150,0.4,3,6,1)
5.4 Heaviside projection filter

This subsection focuses on the implementation of a

5.3 Passive elements black-and-white projection filter. As an example, the
implementation of the filter proposed by Guest et al (2004)
In some cases, certain areas of the design domain may be is explained. This filter is referred to as the Heaviside
required to be void or solid (e.g. to allow for the passage projection filter in the present paper. The aim of the
of a pipe or to support a secondary structure). This can Heaviside projection filter is (1) to achieve a minimum
be easily accomplished by means of the 88 line code length scale in the optimized design, and (2) to obtain black-
through the definition of passive elements, i.e. elements and-white solutions. Guest et al (2004) apply this filter using
 11

nodal design variables, but as shown by Sigmund (2007), it This code will lead to initial physical densities x̄e that do
is equally applicable when element design variables are used not satisfy the volume constraint, which could be avoided
(which is the case in the present paper). by adjusting the initial values of the design variables xe .
The Heaviside filter is a modification of the original In the present code, however, the optimality criteria update
density filter (9) with a Heaviside step function that projects scheme is relied upon to correct the violation of the volume
the density x̃e (from now on called the intermediate density) constraint.
to a physical density x̄e . The physical density x̄e equals one Second, the modification of the sensitivities is accom-
if x̃e > 0 and zero if x̃e = 0. In order to allow for the plished by inserting the following supplementary elseif
use of a gradient-based optimization scheme, the Heaviside statement on line 68:
function is replaced with the following smooth function: elseif ft == 3
dx = beta*exp(-beta*xTilde)+exp(-beta);
x̄e = 1 − e−β x̃e + x̃e e−β (21) dc(:) = H*(dc(:).*dx(:)./Hs);
dv(:) = H*(dv(:).*dx(:)./Hs);

The parameter β controls the smoothness of the approxi- Third, the application of the Heaviside filter to the
mation: for β equal to zero, the Heaviside filter is identical densities is realized by means of the following additional
to the original density filter; for β approaching infinity, the elseif statement, to be inserted on line 78:
approximation approaches a true Heaviside step function. In
elseif ft == 3
order to avoid local minima and to ensure differentiability xTilde(:) = (H*xnew(:))./Hs;
in the optimization, a continuation scheme is used where xPhys = 1-exp(-beta*xTilde)+xTilde*exp(-beta);
the parameter β is gradually increased from 1 to 512 by
doubling its value every 50 iterations or when the change in Finally, the continuation scheme for the regularization
terms of design variables between two consecutive designs parameter β is implemented by inserting the following block
becomes less than 0.01. of code at the end of the optimization loop:
It should be noted that Guest et al (2004) include an if ft == 3 && beta < 512 && ...
extra term in equation (21) to ensure that the lower bound on (loopbeta >= 50 || change <= 0.01)
beta = 2*beta;
the densities x̄e is satisfied; this term is not necessary here
loopbeta = 0;
due to the use of the modified SIMP approach (Sigmund, change = 1;
2007). fprintf(’Parameter beta increased to %g.\n’,beta);
The sensitivities of a function f (x̄e ) with respect to the end
intermediate densities x̃e are obtained by means of the chain The additional counter loopbeta must be initialized and
rule: incremented in the same way as the existing counter loop.
∂f ∂f ∂ x̄e The modified code is used to optimize the MBB beam.
= (22) The same parameter values are used as in subsection 3.4,
∂ x̃e ∂ x̄e ∂ x̃e
except for the filter radius, which is reduced to 0.03 times
where the derivative of the physical density x̄e with respect the width of the design domain. The motivation for this
to the intermediate density x̃e is given by: reduction is that the material resource constraint prohibits
the transformation of the topology obtained in the initial
∂ x̄e phase of the continuation scheme (which is similar to the
= βe−β x̃e + e−β (23)
∂ x̃e topology obtained in subsection 3.4) into a black-and-white
design consisting of bars with a large thickness.
The implementation of the Heaviside filter in the 88 line Figure 8 shows the optimized design obtained with the
code as a third filter option (ft = 3) involves the following three meshes. The optimized design is almost perfectly
modifications. black-and-white and does not exhibit structural details
First, the β parameter (beta) must be defined and the smaller than the filter radius rmin . The Heaviside projection
densities must be filtered before the start of the optimization filter relies on the compact support of the classical filter
loop. To this end, line 47 is replaced with the following function to impose length scale in the solid regions, and
lines: therefore the Heaviside projection cannot be directly applied
beta = 1; with the PDE filter. It can be observed also that the minimum
if ft == 1 || ft == 2 length scale imposed on the material distribution does
xPhys = x; not prevent the occurrence of very small holes. This can
elseif ft == 3
be avoided by using a more advanced filter such as the
xTilde = x;
xPhys = 1-exp(-beta*xTilde)+xTilde*exp(-beta); morphological close-open or open-close filter (Sigmund,
end 2007) or by following a robust approach in the formulation
12

c = 189.14 c = 191.49 c = 193.24

Fig. 8 Optimized design of the MBB beam and corresponding compliance c obtained with the 88 line code extended by a Heaviside filter. A mesh
with 60 × 20 elements (left), 150 × 50 elements (middle), and 300 × 100 elements (right) has been used.

of the optimization problem (Sigmund, 2009; Wang et al, over the first ten iterations of the optimization loop, using
2010). a Lenovo Thinkpad X301 laptop with an Intel Core2 Duo
U9400 processor, 2 GB memory, Windows XP with SP3
(32-bit x86), and MATLAB R2010a. It is clear from the
5.5 Other extensions table that the new 88 line implementation is significantly
faster than the original 99 line code. For the mesh with
The implementation of additional extensions to the 88 line 150 × 50 elements, a factor of 100 speed improvement is
code should be relatively straightforward. The extension accomplished. The 99 line code has not been tested using
to three-dimensional problems may require a lot of the mesh with 300 × 100 elements as the computation time
modifications, but the general structure of the code would becomes excessively large. The computation time for the
remain unchanged. Following the guidelines given by alternative implementations using conv2 based filtering and
Bendsøe and Sigmund (2003), the 88 line code can be PDE based filtering is almost equal to the computation time
converted to a code for mechanism synthesis or for for the 88 line code.
heat conduction problems. The optimality criteria based
optimizer can be replaced with a more versatile optimization
scheme, such as the method of moving asymptotes (MMA) Table 2 Computation time in seconds per iteration for the optimization
of the MBB beam using density filtering.
introduced by Svanberg (1987), in order to enable the
solution of problems with more than one constraint. Mesh size 60 × 20 150 × 50 300 × 100

88 line code 0.12 0.94 5.67

CONV2 based filtering code 0.16 0.78 3.30
6 Performance PDE based filtering code 0.19 1.79 10.08

6.1 Computation time

Table 2 shows the results obtained with a density filter,
In this subsection, the computation time for the 88 line code
using the same configuration as for the sensitivity filter. No
(and the variants presented sections 4 and 5.4) is compared
results are given for the 99 line code as it does not include
with the original 99 line code. The MBB beam optimization
a density filter. The computation time is slightly higher than
problem introduced earlier is considered as a benchmark
for the sensitivity filter, due to the application of the density
problem, using the same parameter values as in subsection
filter in every iteration of the bisection algorithm used to
3.4.
determine the Lagrangian multiplier λ. This is especially
true for the problem with the largest mesh, and for the
Table 1 Computation time in seconds per iteration for the optimization code using PDE based filtering, where the application of
of the MBB beam using sensitivity filtering. the density filter involves a relatively costly backsubstitution
Mesh size 60 × 20 150 × 50 300 × 100 operation. As the PDE based filter is volume-preserving,
this could be avoided using the design variables instead of
99 line code 0.65 75.19 -
88 line code 0.15 0.72 1.85
the physical densities to check the volume constraint in the
CONV2 based filtering code 0.13 0.69 1.98 bisection algorithm. Moreover, the computational cost can
PDE based filtering code 0.13 0.78 2.18 be significantly reduced by employing an iterative solver
(Lazarov and Sigmund, 2010).
Finally, the performance of the code extended by a
Table 1 gives an overview of the computation time in Heaviside filter is described in table 3. Compared to the
seconds per iteration for the optimization of the MBB beam standard density filtering code, the additional Heaviside
using sensitivity filtering. Results are given for four variants projection has no significant impact on the computation
of the optimization code and for three different mesh sizes. time. The computation time per iteration is slightly lower
The computation times have been determined as the average due to the use of a smaller filter radius rmin in the
 13

Table 3 Computation time in seconds per iteration for the optimization to the use of a smaller filter radius rmin in the example
of the MBB beam using Heaviside filtering. problem, the coefficient matrix H becomes sparser, and a
Mesh size 60 × 20 150 × 50 300 × 100 problem with 350 × 117 = 40950 elements could be
solved. The original 99 line code has not been tested, as the
Code with Heaviside filter 0.13 0.86 4.65
computation time becomes prohibitively large for problems
of these dimensions.
example problem, which leads to a sparser coefficient
matrix H, reducing the number of operations required to
multiply the coefficient matrix H with the design variables 7 Conclusion
or the sensitivities. It should be noted, however, that the
use of the Heaviside filter requires the application of a This paper presents a MATLAB code for topology
continuation scheme, which implies that the number of optimization. The code is considered as the successor to the
iterations becomes considerably larger. 99 line code presented by Sigmund (2001). It is published
with the same objective: to provide students and newcomers
to the field with a very simple implementation of a topology
6.2 Memory usage
optimization algorithm that can serve as an introductory
example and as a basis for further developments. The code
While the use of the sparse function to assemble the
can be download from the web site www.topopt.dtu.dk.
stiffness matrix leads to a vast improvement in terms of
computation time, it also increases the program’s memory The major difference with respect to the original 99
footprint. The index vectors iK and jK and the vector sK line code is the computational efficiency. An improvement
with non-zero entries are relatively large and remain in in speed with a factor of 100 has been measured for an
memory throughout the entire optimization loop. Each of example problem with 7500 elements. This has mainly been
these vectors has the size of the element stiffness matrix accomplished by means of loop vectorization and memory
times the number of elements, which is considerably more preallocation.
than the size of the stiffness matrix itself. Moreover, the In addition, the code has been extended by a density
sparse function requires the index vectors iK and jK to filter. The inclusion of a density filter has an important ed-
be defined as double precision arrays, prohibiting the use ucational value, as it paves the road for the implementation
of a more memory efficient integer type. In contrast to of more sophisticated filters such as the Heaviside filter also
the sparse function built-in in MATLAB, the sparse2 discussed in the paper.
function from CHOLMOD does accept integer type index Special care has been taken not to compromise the
vectors. simplicity of the code. As a result, the new code is
In order to get a rough idea of the memory requirements characterized by the same readability as the original 99 line
for the 88 line code and its variants, an informal test has code, although the number of lines has been reduced to 88.
been conducted: the example problem from the previous The paper also presents two alternative implementations.
subsection has been solved multiple times, each time The first alternative takes advantage of the conv2 function
incrementing the mesh size, until the computer ran out built-in in MATLAB to filter densities and sensitivities, so
of memory. The same computer has been used as for reducing the number of lines to 71 without affecting the
the determination of the computation time. The test has computational cost or the readability of the code (for those
been performed using both sensitivity and density filtering, familiar with the conv2 function). The second alternative
leading to identical results. uses a filter based on a Helmholtz type differential equation,
The largest problem that could be solved with the 88 allowing for the use of a finite element solver to perform
line code consisted of 300 × 100 = 30000 elements. The the filtering operation. This is beneficial for problems with
code using conv2 based filtering requires less memory a complex geometry or when the optimization problem is
as it avoids the definition of the (sparse but nonetheless solved in parallel.
large) coefficient matrix H, using a small matrix h instead
that represents the non-zero part of the filter kernel. As a Acknowledgements This work was financially supported by the
consequence, this code allowed for the solution of a problem Eurohorcs/ESF European Young Investigator Award (EURYI), by a
with 700 × 233 = 163100 elements. The code based Center of Advanced User Support (CAUS) grant from the Danish
Center of Scientific Computing (DCSC), and by an Elite Research
on PDE filtering performs in between and allowed us to
Prize from the Danish Minister of Research. The third author is a
solve a problem with 600 × 200 = 120000 elements. The postdoctoral fellow of the Research Foundation - Flanders and a
extension of the 88 line code with a Heaviside filter has member of K.U.Leuven-BOF PFV/10/002 OPTEC-Optimization in
no noticeable influence on memory usage; however, due Engineering Center.
14

References Lazarov B, Sigmund O (2010) Filters in topology opti-

mization based on Helmholtz type differential equations
Alberty J, Carstensen C, Funken S (1999) Remarks around (submitted for publication)
50 lines of Matlab: short finite element implementation. Sigmund O (1994) Design of material structures using
Numerical Algorithms 20(2–3):117–137 topology optimization. PhD thesis, DCAMM S-report
Alberty J, Carstensen C, Funken S, Klose R (2002) Matlab S69, Department of Solid Mechanics, Technical Univer-
implementation of the finite element method in elasticity. sity of Denmark
Computing 69(3):239–263 Sigmund O (1997) On the design of compliant mechanisms
Allaire G (2009) Shape and topology using topology optimization. Mechanics of Structures and
optimization by the level set method. URL Machines 25(4):493–524
http://www.cmap.polytechnique.fr/~allaire Sigmund O (2001) A 99 line topology optimization
Bendsøe M, Sigmund O (2003) Topology Optimization. code written in Matlab. Structural and Multidisciplinary
Theory, Methods and Applications. Springer Optimization 21(2):120–127
Bendsøe M (1989) Optimal shape design as a material Sigmund O (2007) Morphology-based black and white
distribution problem. Structural optimization 1:193–202 filters for topology optimization. Structural and Multidis-
Bourdin B (2001) Filters in topology optimization. Inter- ciplinary Optimization 33(4-5):401–424
national Journal for Numerical Methods in Engineering Sigmund O (2009) Manufacturing tolerant topology opti-
50(9):2143–2158 mization. Acta Mechanica Sinica 25(2):227–239
Bruns TE, Tortorelli DA (2001) Topology optimization of Sigmund O, Petersson J (1998) Numerical instabilities
non-linear elastic structures and compliant mechanisms. in topology optimization: A survey on procedures
Computer Methods in Applied Mechanics and Engineer- dealing with checkerboards, mesh-dependencies and
ing 190(26-27):3443–3459 local minima. Structural Optimization 16(1):68–75
Challis VJ (2010) A discrete level-set topology optimization Suresh K (2010) A 199-line Matlab code for Pareto-
code written in Matlab. STRUCTURAL AND MULTI- optimal tracing in topology optimization. Structural and
DISCIPLINARY OPTIMIZATION 41(3):453–464 Multidisciplinary Optimization Published online
Dabrowski M, Krotkiewski M, Schmid D (2008) MIL- Svanberg K (1987) Method of moving asymptotes - a new
AMIN: MATLAB-based finite element method solver for method for structural optimization. International Journal
large problems. Geochemistry Geophysics Geosystems for Numerical Methods in Engineering 24(2):359–373
9(4) The MathWorks (2010) MATLAB Programming Funda-
Davis T (2007) Creating sparse finite-element matrices in mentals
MATLAB. Guest blog in Loren on the Art of MATLAB, Wang F, Lazarov B, Sigmund O (2010) On projection meth-
http://blogs.mathworks.com/loren/2007/03/01/creating- ods, convergence and robust formulations in topology
sparse-finite-element-matrices-in-matlab/. URL optimization (submitted for publication)
http://blogs.mathworks.com/loren/ Zhou M, Rozvany G (1991) The COC algorithm, part
Davis T (2008) User Guide for CHOLMOD: a sparse II: Topological, geometrical and generalized shape
Cholesky factorization and modification package. De- optimization. Computer Methods in Applied Mechanics
partment of Computer and Information Science and and Engineering 89(1–3):309–336
Engineering, University of Florida, Gainesville, FL, USA
Dı́az A, Sigmund O (1995) Checkerboard patterns in layout
optimization. Structural Optimization 10(1):40–45
Guest J, Prevost J, Belytschko T (2004) Achieving minimum
length scale in topology optimization using nodal design
variables and projection functions. International Journal
for Numerical Methods in Engineering 61(2):238–254
Jog C, Haber R (1996) Stability of finite element models
for distributed-parameter optimization and topology
design. Computer Methods in Applied Mechanics and
Engineering 130(3–4):203–226
Lazarov B, Sigmund O (2009) Sensitivity filters in topology
optimisation as a solution to helmholtz type differential
equation. In: In Proc. of the 8th World Congress on
Structural and Multidisciplinary Optimization
 15

Appendix - MATLAB code

1 %%%% AN 88 LINE TOPOLOGY OPTIMIZATION CODE %%%%

2 function top88(nelx,nely,volfrac,penal,rmin,ft)
3 %% MATERIAL PROPERTIES
4 E0 = 1;
5 Emin = 1e-9;
6 nu = 0.3;
7 %% PREPARE FINITE ELEMENT ANALYSIS
8 A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
9 A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
10 B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
11 B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
12 KE = 1/(1-nu^2)/24*([A11 A12;A12’ A11]+nu*[B11 B12;B12’ B11]);
13 nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);
14 edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,nelx*nely,1);
15 edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1);
16 iK = reshape(kron(edofMat,ones(8,1))’,64*nelx*nely,1);
17 jK = reshape(kron(edofMat,ones(1,8))’,64*nelx*nely,1);
18 % DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
19 F = sparse(2,1,-1,2*(nely+1)*(nelx+1),1);
20 U = zeros(2*(nely+1)*(nelx+1),1);
21 fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
22 alldofs = [1:2*(nely+1)*(nelx+1)];
23 freedofs = setdiff(alldofs,fixeddofs);
24 %% PREPARE FILTER
25 iH = ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2,1);
26 jH = ones(size(iH));
27 sH = zeros(size(iH));
28 k = 0;
29 for i1 = 1:nelx
30 for j1 = 1:nely
31 e1 = (i1-1)*nely+j1;
32 for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
33 for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
34 e2 = (i2-1)*nely+j2;
35 k = k+1;
36 iH(k) = e1;
37 jH(k) = e2;
38 sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2));
39 end
40 end
41 end
42 end
43 H = sparse(iH,jH,sH);
44 Hs = sum(H,2);
45 %% INITIALIZE ITERATION
46 x = repmat(volfrac,nely,nelx);
47 xPhys = x;
48 loop = 0;
49 change = 1;
50 %% START ITERATION
51 while change > 0.01
52 loop = loop + 1;
53 %% FE-ANALYSIS
54 sK = reshape(KE(:)*(Emin+xPhys(:)’.^penal*(E0-Emin)),64*nelx*nely,1);
55 K = sparse(iK,jK,sK); K = (K+K’)/2;
56 U(freedofs) = K(freedofs,freedofs)\F(freedofs);
57 %% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
58 ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),nely,nelx);
59 c = sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce));
60 dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce;
61 dv = ones(nely,nelx);
62 %% FILTERING/MODIFICATION OF SENSITIVITIES
63 if ft == 1
64 dc(:) = H*(x(:).*dc(:))./Hs./max(1e-3,x(:));
65 elseif ft == 2
66 dc(:) = H*(dc(:)./Hs);
67 dv(:) = H*(dv(:)./Hs);
68 end
69 %% OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES AND PHYSICAL DENSITIES
 16

70 l1 = 0; l2 = 1e9; move = 0.2;

71 while (l2-l1)/(l1+l2) > 1e-3
72 lmid = 0.5*(l2+l1);
73 xnew = max(0,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
74 if ft == 1
75 xPhys = xnew;
76 elseif ft == 2
77 xPhys(:) = (H*xnew(:))./Hs;
78 end
79 if sum(xPhys(:)) > volfrac*nelx*nely, l1 = lmid; else l2 = lmid; end
80 end
81 change = max(abs(xnew(:)-x(:)));
82 x = xnew;
83 %% PRINT RESULTS
84 fprintf(’ It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n’,loop,c, ...
85 mean(xPhys(:)),change);
86 %% PLOT DENSITIES
87 colormap(gray); imagesc(1-xPhys); caxis([0 1]); axis equal; axis off; drawnow;
88 end