Educational article Struct Multidisc Optim 21, 120–127  Springer-Verlag 2001

A 99 line topology optimization code written in Matlab

O. Sigmund

Abstract The paper presents a compact Matlab im- fraction and then the program shows the correct optimal
plementation of a topology optimization code for com- topology for comparison.
pliance minimization of statically loaded structures. The In the literature, one can find a multitude of approaches
total number of Matlab input lines is 99 including opti- for the solving of topology optimization problems. In the
mizer and Finite Element subroutine. The 99 lines are original paper Bendsøe and Kikuchi (1988) a so-called
divided into 36 lines for the main program, 12 lines for the microstructure or homogenization based approach was
Optimality Criteria based optimizer, 16 lines for a mesh- used, based on studies of existence of solutions.
independency filter and 35 lines for the finite element The homogenization based approach has been adopted
code. In fact, excluding comment lines and lines associ- in many papers but has the disadvantage that the deter-
ated with output and finite element analysis, it is shown mination and evaluation of optimal microstructures and
that only 49 Matlab input lines are required for solving their orientations is cumbersome if not unresolved (for
a well-posed topology optimization problem. By adding noncompliance problems) and furthermore, the resulting
three additional lines, the program can solve problems structures cannot be built since no definite length-scale
with multiple load cases. The code is intended for edu- is associated with the microstructures. However, the ho-
cational purposes. The complete Matlab code is given in mogenization approach to topology optimization is still
the Appendix and can be down-loaded from the web-site important in the sense that it can provide bounds on the
http://www.topopt.dtu.dk. theoretical performance of structures.
An alternative approach to topology optimization is
Key words topology optimization, education, optimal- the so-called “power-law approach” or SIMP approach
ity criteria, world-wide web, Matlab code (Solid Isotropic Material with Penalization) (Bendsøe
1989; Zhou and Rozvany 1991; Mlejnek 1992). Here, ma-
terial properties are assumed constant within each elem-
ent used to discretize the design domain and the variables
are the element relative densities. The material proper-
1 ties are modelled as the relative material density raised
Introduction to some power times the material properties of solid ma-
terial. This approach has been criticized since it was ar-
The Matlab code presented in this paper is intended gued that no physical material exists with properties de-
for engineering education. Students and newcomers to scribed by the power-law interpolation. However, a recent
the field of topology optimization can down-load the paper by Bendsøe and Sigmund (1999) proved that the
code from the web-page http://www.topopt.dtu.dk. power-law approach is physically permissible as long as
The code may be used in courses in structural optimiza- simple conditions on the power are satisfied (e.g. p ≥ 3
tion where students may be assigned to do extensions for Poisson’s ratio equal to 13 ). To ensure existence of so-
such as multiple load-cases, alternative mesh-independ- lutions, the power-law approach must be combined with
ency schemes, passive areas, etc. Another possibility is to a perimeter constraint, a gradient constraint or with fil-
use the program to develop students’ intuition for optimal tering techniques (see Sigmund and Petersson 1998, for
design. Advanced students may be asked to guess the op- an overview). The power-law approach to topology op-
timal topology for given boundary condition and volume timization has been applied to problems with multiple
constraints, multiple physics and multiple materials.
Received October 22, 1999
Whereas the solution of the above mentioned ap-
O. Sigmund proaches is based on mathematical programming tech-
Department of Solid Mechanics, Building 404, Technical Uni- niques and continuous design variables, a number of pa-
versity of Denmark, DK-2800 Lyngby, Denmark pers have appeared on solving the topology optimization
e-mail: sigmund@fam.dtu.dk problem as an integer problem. Beckers (1999) success-
 121

fully solved large-scale compliance minimization prob- A topology optimization problem based on the power-
lems using a dual-approach but other approaches based law approach, where the objective is to minimize compli-
on genetic algorithms or other semi-random approaches ance can be written as
require thousands of function evaluations even for small N 


number of elements and must be considered impractical. min T
: c(x) = U KU = (xe ) ue k0 ue 
p T 

x 

Apart from above mentioned approaches, which all e=1 

solve well defined problems (e.g. minimization of com- V (x) 
subject to : V0 = f ,
pliance) a number of heuristic or intuition based ap- 



proaches have been shown to decrease compliance or : KU = F 



other objective functions. Among these methods are so- 
called evolutionary design methods (see e.g. Xie and : 0 < xmin ≤ x ≤ 1
Steven 1997; Baumgartner et al. 1992). Apart from be- (1)
ing very easy to understand and implement (at least
where U and F are the global displacement and force
for the compliance minimization case), the main moti-
vectors, respectively, K is the global stiffness matrix, ue
vation for the evolutionary approaches seems to be that
and ke are the element displacement vector and stiffness
mathematically based or continuous variable approaches
matrix, respectively, x is the vector of design variables,
“involve some complex calculus operations and mathe-
xmin is a vector of minimum relative densities (non-zero
matical programming” (citation from Li et al. 1999) and
to avoid singularity), N (= nelx × nely) is the number
they contain “mathematical methods of some complex-
of elements used to discretize the design domain, p is the
ity” (citation from Zhao et al. 1998) whereas the evo-
penalization power (typically p = 3), V (x) and V0 is the
lutionary approach “takes advantage of powerful com-
material volume and design domain volume, respectively
puting technology and intuitive concepts of evolution
and f (volfrac) is the prescribed volume fraction.
processes in nature” (citation from Li et al. 1999). Two
The optimization problem (1) could be solved using
things can be argued against this. First, the evolutionary
several different approaches such as Optimality Criteria
approaches become complicated themselves, once more
(OC) methods, Sequential Linear Programming (SLP)
complex objectives than compliance minimization are
methods or the Method of Moving Asymptotes (MMA by
considered and second, as shown in this paper, the “math-
Svanberg 1987) and others. For simplicity, we will here
ematically based” approaches for compliance minimiza-
use a standard OC-method.
tion are simple to implement as well and are compu-
Following Bendsøe (1995) a heuristic updating scheme
tationally equally efficient. Furthermore, mathematical
for the design variables can be formulated as
programming based methods can easily be extended to
other non-compliance objectives such as non-self-adjoint xnew
e =
and multiphysics problems and to problems with multiple

constraints (e.g. Sigmund 1999). Extensions of the evolu-
max(xmin , xe − m)


tionary approach to such cases seem more questionable. 
 if xe Beη ≤ max(xmin , xe − m) ,


The complete Matlab code is given in the Appendix. 

The remainder of the paper consists of definition and 
x B η
e e
discussion of the optimization problem (Sect. 2), com-

 if max(xmin , xe − m) < xe Beη < min(1, xe + m) ,
ments about the Matlab implementation (Sect. 3) fol- 



lowed by a discussion of extensions (Sect. 4) and a conclu- 
min(1, xe + m)


sion (Sect. 5). 
if min(1, xe + m) ≤ xe Beη ,
(2)
2
where m (move) is a positive move-limit, η (= 1/2) is
The topology optimization problem
a numerical damping coefficient and Be is found from the
optimality condition as
A number of simplifications are introduced to simplify the
Matlab code. First, the design domain is assumed to be ∂c
−
rectangular and discretized by square finite elements. In ∂xe
this way, the numbering of elements and nodes is simple Be = , (3)
∂V
(column by column starting in the upper left corner) and λ
the aspect ratio of the structure is given by the ratio of ∂xe
elements in the horizontal (nelx) and the vertical direc- where λ is a Lagrangian multiplier that can be found by
tion (nely).1 a bi-sectioning algorithm.
The sensitivity of the objective function is found as
1
Names in type-writer style refer to Matlab variable names
∂c
that differ from the obvious (see the Matlab code in the = −p(xe )p−1 uTe k0 ue . (4)
Appendix) ∂xe
122

For more details on the derivation and implementation of

the optimality criteria method, the reader is referred to
the literature (e.g. Bendsøe 1995).
In order to ensure existence of solutions to the top-
ology optimization problem (1), some sort of restriction ?
on the resulting design must be introduced (see Sigmund 111
000
and Petersson 1998, for an overview). Here we use a fil- 111
000 000
111
tering technique (Sigmund 1994, 1997). It must be em- 1
0
0
1
phasized that this filter has not yet been proven to en- 0
1
0
1
sure existence of solutions, but numerous applications
by the author have proven the the filter produces mesh-
0
1
0
1
0
1
0
1
?
independent designs in practice. 0
1
The mesh-independency filter works by modifying the 0
1 111
000
element sensitivities as follows:
N


∂c 1 ∂c
= Ĥf xf . (5) 111
000
∂xe N

f =1
∂xf 111
000 000
111
xe Ĥf
Fig. 1 Topology optimization of the MBB-beam. Top: full
f =1
design domain, middle: half design domain with symmetry
boundary conditions and bottom: resulting topology opti-
The convolution operator (weight factor) Ĥf is written as mized beam (both halves)

Ĥf = rmin − dist(e, f ) ,

The default boundary conditions correspond to half of the
{f ∈ N | dist(e, f ) ≤ rmin }, e = 1, . . . , N , (6) “MBB-beam” (Fig. 1). The load is applied vertically in
the upper left corner and there is symmetric boundary
where the operator dist(e, f ) is defined as the distance be- conditions along the left edge and the structure is sup-
tween centre of element e and centre of element f . The ported horizontally in the lower right corner.
convolution operator Ĥf is zero outside the filter area. Important details of the Matlab code are discussed in
The convolution operator decays linearly with the dis- the following subsections.
tance from element f . Instead of the original sensitivities
(4), the modified sensitivities (5) are used in the Optimal- 3.1
ity Criteria update (3). Main program (lines 1–36)

The main program (lines 1–36) starts by distributing the

3 material evenly in the design domain (line 4). After some
Matlab implementation other initializations, the main loop starts with a call to
the Finite Element subroutine (line 12) which returns the
The Matlab code (see the Appendix), is built up as a stan- displacement vector U. Since the element stiffness matrix
dard topology optimization code. The main program is for solid material is the same for all elements, the elem-
called from the Matlab prompt by the line ent stiffness matrix subroutine is called only once (line
14). Following this, a loop over all elements (lines 16–
24) determines objective function and sensitivities (4).
top(nelx,nely,volfrac,penal,rmin)
The variables n1 and n2 denote upper left and right
element node numbers in global node numbers and are
where nelx and nely are the number of elements in the used to extract the element displacement vector Ue from
horizontal and vertical directions, respectively, volfrac the global displacement vector U. The sensitivity analy-
is the volume fraction, penal is the penalization power sis is followed by a call to the mesh-independency filter
and rmin is the filter size (divided by element size). Other (line 26) and the Optimality Criteria optimizer (line 28).
variables as well as boundary conditions are defined in The current compliance as well as other parameters are
the Matlab code itself and can be edited if needed. For printed by lines 30–33 and the resulting density distri-
each iteration in the topology optimization loop, the code bution is plotted (line 35). The main loop is terminated
generates a picture of the current density distribution. if the change in design variables (change determined in
Figure 1 shows the resulting density distribution obtained line 30) is less than 1 percent2 . Otherwise above steps are
by the code given in the Appendix called with the input repeated.
line
2
this is a rather “sloppy” convergence criterion and could
top(60,20,0.5,3.0,1.5) be decreased if needed
 123

3.2 1
0 1
0
0
1 0
1
Optimality criteria based optimizer (lines 37–48) 0
1 0
1
The updated design variables are found by the opti-
0
1
0
1
0
1
0
1
? 0
1
0
1
0
1
0
1
mizer (lines 37–48). Knowing that the material volume 0
1 0
1
(sum(sum(xnew))) is a monotonously decreasing func-
tion of the Lagrange multiplier (lag), the value of the
Fig. 2 Topology optimization of a cantilever beam. Left: de-
Lagrangian multiplier that satisfies the volume constraint sign domain and right: topology optimized beam
can be found by a bi-sectioning algorithm (lines 40-48).
The bi-sectioning algorithm is initialized by guessing
a lower l1 and an upper l2 bound for the Lagrangian
multiplier (line 39). The interval which bounds the La- the difference between all degrees of freedom and the fixed
grangian multiplier is repeatedly halved until its size is degrees of freedom (line 82).
less than the convergence criteria (line 40). The element stiffness matrix is calculated in lines 86–
99. The 8 by 8 matrix for a square bi-linear 4-node elem-
ent was determined analytically using a symbolic manip-
3.3 ulation software. The Young’s modulus E and the Pois-
Mesh-independency filtering (lines 49–64) son’s ratio nu can be altered in lines 88 and 89.

Lines 49–64 represent the Matlab implementation of (5).

Note that not all elements in the design domain are
4
searched in order to find the elements that lie within
Extensions
the radius rmin but only those within a square with side
lengths two times round(rmin) around the considered
element. By selecting rmin less than one in the call of the The Matlab code given in the Appendix solves the prob-
routine, the filtered sensitivities will be equal to the ori- lem of optimizing the material distribution in the MBB-
ginal sensitivities making the filter inactive. beam (Fig. 1) such that its compliance is minimized.
A number of extensions and changes in the algorithm
can be thought of, a few of which are mentioned in the
3.4 following.
Finite element code (lines 65–99)

The finite element code is written in lines 65–99. Note 4.1

that the solver makes use of the sparse option in Mat- Other boundary conditions
lab. The global stiffness matrix is formed by a loop over
all elements (lines 70–77). As was the case in the main It is very simple to change boundary conditions and sup-
program, variables n1 and n2 denote upper left and right port conditions in order to solve other optimization prob-
element node numbers in global node numbers and are lems. In order to solve the short cantilever example shown
used to insert the element stiffness matrix at the right in Fig. 2, only lines 79 and 80 must be changed to
places in the global stiffness matrix.
As mentioned before, both nodes and elements are 79 F(2*(nelx+1)*(nely+1),1) = -1;
numbered column wise from left to right. Furthermore, 80 fixeddofs = [1:2*(nely+1)];
each node has two degrees of freedom (horizontal and ver-
tical), thus the command F(2,1)=-1. (line 79) applies
a vertical unit force force in the upper left corner. With these changes, the input line for the case shown
Supports are implemented by eliminating fixed de- in Fig. 2 is
grees of freedom from the linear equations. Matlab can do
this very elegantly with the line top(32,20,0.4,3.0,1.2)

84 U(freedofs,:) = K(freedofs,freedofs) \
F(freedofs,:); 4.2
Multiple load cases
where freedofs indicate the degrees of freedom which
are unconstrained. Mostly, it is easier to define the de- It is also very simple to extend the algorithm to account
grees of freedom that are fixed (fixeddofs) thereafter the for multiple load cases. In fact, this can be done by adding
freedofs are found automatically using the Matlab oper- only three additional lines and making minor changes to
ator setdiff which finds the free degrees of freedoms as another 4 lines.
124

In the case of two load cases, force and displacement

1 1 1
vectors must be defined as two-column vectors which
means that line 69 is changed to 1
0
0
1 1
0
0
1 1
0
0
1
0
1 0
1 0
1
69 F = sparse(2*(nely+1)*(nelx+1),2);
U = sparse(2*(nely+1)*(nelx+1),2);
0
1
0
1
0
1
0
1
? 0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1 0
1 0
1
2 1 2
The objective function is now the sum of two compli-
ances, i.e. Fig. 3 Topology optimization of a cantilever beam with two
load-cases. Left: design domain, middle: topology optimized
beam using one load case and right: topology optimized beam
2

using two load cases
c(x) = UTi KUi (7)
i=1

0110 0110
thus lines 20–22 are substituted with the lines 1010 1010
19b dc(ely,elx) = 0.;
1010
1010
? 1010
1010
19c for i = 1:2
20 Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2;
2*n2+1;2*n2+2;2*n1+1;2*n1+2],i);
Fig. 4 Topology optimization of a cantilever beam with
21 c = c + x(ely,elx)^penal*Ue’*KE*Ue;
a fixed hole. Left: design domain and right: topology opti-
22 dc(ely,elx) = dc(ely,elx) - mized beam
penal*x(ely,elx)^(penal-1)*Ue’*KE*Ue;
22b end

Figure 4 shows the resulting structure obtained with

To solve the two-load problem indicated in Fig. 3, a unit the input
upward load in the top-right corner is added to line 79,
which then becomes top(45,30,0.5,3.0,1.5),

79 F(2*(nelx+1)*(nely+1),1) = -1.;
F(2*(nelx)*(nely+1)+2,2) = 1.; when the following 10 lines were added to the main
program (after line 4) in order to find passive elem-
ents within a circle with radius nely/3. and center
The input line for Fig. 3 is (nely/2., nelx/3.)

top(30,30,0.4,3.0,1.2). for ely = 1:nely

for elx = 1:nelx
if sqrt((ely-nely/2.)^2+(elx-nelx/3.)^2) <
nely/3.
4.3 passive(ely,elx) = 1;
Passive elements x(ely,elx) = 0.001;
else
In some cases, some of the elements may be required to passive(ely,elx) = 0;
take the minimum density value (e.g. a hole for a pipe). end
An nely × nelx array passive with zeros at elements end
free to change and ones at elements fixed to be zero can end
be defined in the main program and transferred to the OC
subroutine (adding passive to the call in lines 28 and 38).
The added line 4.4
Alternative optimizer
42b xnew(find(passive)) = 0.001;
Admittedly, the optimality criteria based optimizer im-
plemented here is only good for a single constraint and it
in the OC subroutine looks for passive elements and sets is based on a heuristic fixed point type updating scheme.
their density equal to the minimum density (0.001). In order to install a better optimizer, one can obtain (free
 125

of charge for academic purposes) the Matlab version of http: //www.mekanik.ikp.liu.se/andridiv/matlab/

the MMA-algorithm (Svanberg 1987) from Krister Svan- theory.html.
berg, KTH, Sweden. The MMA code is called with the The code was intentionally kept compact in order
following input line to keep the total number of lines below 100. If users
of the code should find ways to further compactify or
mmasub(INPUT-variables, ... , simplify the code, the author would be happy to re-
OUTPUT-variables) ceive suggested modifications that can be implemented
in the public domain code (the author’s e-mail address is
sigmund@fam.dtu.dk).
where the total number of input/output variables is 20, Since its first publication on the World Wide Web in
including objective function, constraints, old and new October 1999, the Matlab code has been down-loaded
densities, etc. Implementing the MMA-optimizer is fairly more than 500 times by different users (as of August
simple, but requires the definition of several auxiliary 2000). Among other positive feedbacks, several profes-
variables. However, it allows for the solving of more com- sors reported that they have used the code in courses on
plex design problems with more than one constraint. The structural optimization and have let their students imple-
Matlab optimizer will solve the standard topology op- ment alternative boundary conditions and multiple load
timization problem using less iterations at the cost of cases.
a slightly increased CPU-time per iteration.
Acknowledgements This work was supported by the Dan-
ish Technical Research Council through the THOR/Talent-
4.5
programme: Design of MicroElectroMechanical Systems
Other extensions
(MEMS). The author would also like to thank Thomas Buhl,
Technical University of Denmark, for his inputs to an earlier
Extensions to three dimensions should be straight for- version of the code.
ward whereas more complex problems such as compliant
mechanism design (Sigmund 1997) requires the imple-
mentation of the MMA optimizer and the definition of ex-
tra constraints. The simplicity of the Matlab commands References
allow for easy extensions of the graphical output, interac-
tive input etc. Baumgartner, A., Harzheim, L.; Mattheck, C. 1992: SKO
(Soft Kill Option): The biological way to find an optimum
structure topology. Int. J. Fatigue 14, 387–393
5
Conclusions Beckers, M. 1999: Topology optimization using a dual method
with discrete variables. Struct. Optim. 17, 14–24

This paper has presented a very simple implementation Bendsøe, M.P. 1989: Optimal shape design as a material dis-
of a mathematical programming base topology optimiza- tribution problem. Struct. Optim. 1, 193–202
tion algorithm. The code is implemented using only 99
Bendsøe, M.P. 1995: Optimization of structural topology,
Matlab input lines and includes optimizer, mesh-indepe- shape and material . Berlin, Heidelberg, New York: Springer
ndency filtering and Finite Element code.
The Matlab code can be down-loaded from the web- Bendsøe, M.P.; Kikuchi, N. 1988: Generating optimal topolo-
page http://www.topopt.dtu.dk and is intended for ed- gies in optimal design using a homogenization method. Comp.
ucational purposes. The code can easily be extended to Meth. Appl. Mech. Engrg. 71, 197–224
include multi load problems and the definition of passive Bendsøe, M.P.; Sigmund, O. 1999: Material interpolations in
areas. topology optimization. Arch. Appl. Mech. 69, 635–654
Running the code in Matlab is rather slow compared
to a Fortran implementation of the same code which can Li, Q.; Steven, G.P.; Xie, Y. M. 1999: On equivalence between
be tested at the web-site http://www.topopt.dtu.dk. stress criterion and stiffness criterion in evolutionary struc-
However, an add-on package to Matlab (MATLAB Com- tural optimization. Struct. Optim. 18, 67–73
piler) allows for the generation of more efficient C-code Mlejnek, H.P. 1992: Some aspects of the genesis of structures.
that can be optimized for run-time (this option, how- Struct. Optim. 5, 64–69
ever, has not been tested by the author). It should be
noted that speed can be gained by modifying the Mat- Sigmund, O. 1994: Design of material structures using top-
lab code itself, however the speed is gained on the cost of ology optimization. Ph.D. Thesis, Department of Solid Me-
chanics, Technical University of Denmark
simplicity of the program. The modification is suggested
by Andreas Rietz from Linköping University who uses Sigmund, O. 1997: On the design of compliant mechanisms
sparsity options in the assembly of the global stiffness ma- using topology optimization. Mech. Struct. Mach. 25, 495–
trix. The reader may down-load his code at the web-page: 526
126

Sigmund, O. 1999: Topology optimization of multi-physics, sprintf(’%10.4f’,c) ...

multi-material structures. Proc. WCSMO-3 (held in Buffalo, 32 ’ Vol.: ’ sprintf(’%6.3f’,sum(sum(x))/
NY) (nelx*nely)) ...
33 ’ ch.: ’ sprintf(’%6.3f’,change )])
Sigmund, O.; Petersson, J. 1998: Numerical instabilities in 34 % PLOT DENSITIES
topology optimization: a survey on procedures dealing with 35 colormap(gray); imagesc(-x); axis equal; axis
checkerboards, mesh-dependencies and local minima. Struct. tight; axis off;pause(1e-6);
Optim. 16, 68–75 36 end
37 %%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%
Svanberg, K. 1987: The method of moving asymptotes – a new 38 function [xnew]=OC(nelx,nely,x,volfrac,dc)
method for structural optimization. Int. J. Numer. Meth. En- 39 l1 = 0; l2 = 100000; move = 0.2;
grg. 24, 359–373 40 while (l2-l1 > 1e-4)
41 lmid = 0.5*(l2+l1);
Xie, Y.M.; Steven, G.P. 1997: Evolutionary structural opti-
42 xnew = max(0.001,max(x-move,min(1.,min(x+move,x.
mization. Berlin, Heidelberg, New York: Springer *sqrt(-dc./lmid)))));
43 if sum(sum(xnew)) - volfrac*nelx*nely > 0;
Zhao, C.; Hornby, P.; Steven, G.P.; Xie, Y.M. 1998: A gener-
44 l1 = lmid;
alized evolutionary method for numerical topology optimiza-
45 else
tion of structures under static loading conditions. Struct. Op-
46 l2 = lmid;
tim. 15, 251–260
47 end
48 end
Zhou, M.; Rozvany, G.I.N. 1991: The COC algorithm, part
II: Topological, geometry and generalized shape optimization. 49 %%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%
Comp. Meth. Appl. Mech. Engrng. 89, 197–224 50 function [dcn]=check(nelx,nely,rmin,x,dc)
51 dcn=zeros(nely,nelx);
52 for i = 1:nelx
53 for j = 1:nely
6 54 sum=0.0;
Appendix – Matlab code 55 for k = max(i-round(rmin),1):
min(i+round(rmin),nelx)
56 for l = max(j-round(rmin),1):
min(j+round(rmin), nely)
1 %%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE
57 fac = rmin-sqrt((i-k)^2+(j-l)^2);
SIGMUND, OCTOBER 1999 %%%
58 sum = sum+max(0,fac);
2 function top(nelx,nely,volfrac,penal,rmin);
59 dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)
3 % INITIALIZE
*dc(l,k);
4 x(1:nely,1:nelx) = volfrac;
60 end
5 loop = 0;
61 end
6 change = 1.;
62 dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
7 % START ITERATION
63 end
8 while change > 0.01
64 end
9 loop = loop + 1;
10 xold = x; 65 %%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%
11 % FE-ANALYSIS 66 function [U]=FE(nelx,nely,x,penal)
12 [U]=FE(nelx,nely,x,penal); 67 [KE] = lk;
13 % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS 68 K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*
14 [KE] = lk; (nely+1));
15 c = 0.; 69 F = sparse(2*(nely+1)*(nelx+1),1); U =
16 for ely = 1:nely sparse(2*(nely+1)*(nelx+1),1);
17 for elx = 1:nelx 70 for ely = 1:nely
18 n1 = (nely+1)*(elx-1)+ely; 71 for elx = 1:nelx
19 n2 = (nely+1)* elx +ely; 72 n1 = (nely+1)*(elx-1)+ely;
20 Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1; 73 n2 = (nely+1)* elx +ely;
2*n2+2; 2*n1+1;2*n1+2],1); 74 edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1;
21 c = c + x(ely,elx)^penal*Ue’*KE*Ue; 2*n2+2;2*n1+1; 2*n1+2];
22 dc(ely,elx) = -penal*x(ely,elx)^(penal-1)* 75 K(edof,edof) = K(edof,edof) +
Ue’*KE*Ue; x(ely,elx)^penal*KE;
23 end 76 end
24 end 77 end
25 % FILTERING OF SENSITIVITIES 78 % DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
26 [dc] = check(nelx,nely,rmin,x,dc); 79 F(2,1) = -1;
27 % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD 80 fixeddofs = union([1:2:2*(nely+1)],
28 [x] = OC(nelx,nely,x,volfrac,dc); [2*(nelx+1)*(nely+1)]);
29 % PRINT RESULTS 81 alldofs = [1:2*(nely+1)*(nelx+1)];
30 change = max(max(abs(x-xold))); 82 freedofs = setdiff(alldofs,fixeddofs);
31 disp([’ It.: ’ sprintf(’%4i’,loop) ’ Obj.: ’ 83 % SOLVING
 127

84 U(freedofs,:) = K(freedofs,freedofs) \ 92 KE = E/(1-nu^2)*

F(freedofs,:); [ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
85 U(fixeddofs,:)= 0; 93 k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
86 %%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%% 94 k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
87 function [KE]=lk 95 k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
88 E = 1.; 96 k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
89 nu = 0.3; 97 k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
90 k=[ 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 98 k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
91 -1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8]; 99 k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];