To Refine or Not to Refine: Topology Optimization of Adaptively
Refined Infill Structures for Additive Manufacturing
Jun Wu

arXiv:1712.04523v1 [math.NA] 12 Dec 2017

Depart. of Design Engineering, TU Delft, The Netherlands

Abstract
We present a novel method for optimizing structures by adaptively refining the structural details.
The adaptively refined structures topologically resemble the graph composed of edges in the
quadtree mesh. However, while adaptive mesh refinement is employed in numerical analysis for
reducing computational complexity, in this work we interpret the edges as structural elements
carrying mechanical loads. The adaptivity and full coverage over the design domain make adaptively refined structures well-suited as infill for 3D printed parts, where uniform infill structures
have been typically used in practice.
The topology optimization of adaptively refined structures is realized based on two novel
ideas. First, we relax the binary design variables (i.e., to refine or not to refine) by using continuous variables, giving rise to sensitivity analysis for an efficient gradient-based optimization.
Second, we propose a refinement filter to encode the dependence of design variables among multiple levels in the structural hierarchy. The refinement filter thus enables considering the design
variables on all levels simultaneously in the optimization. Our numerical results demonstrate
optimized structures spanning multiple levels in the quadtree, with cell sizes smoothly varying
in the closed-walled design domain.
Keywords: Topology optimization, adaptive structures, additive manufacturing, infill

1. Introduction
1.1. Background: Topology Optimization for Additive Manufacturing
Topology optimization is a promising design method for additive manufacturing, as it fully
exploits the manufacturing flexibility enabled by the additive process. It distributes a limited
amount of material in the design domain with the objective of maximizing the structural performance [1]. Reviews of topology optimization techniques can be found in [2, 3]. Additive
manufacturing (AM, also known as 3D printing) provides an effective means to fabricate such
optimal structures, which are often too complex to be produced by conventional manufacturing
such as machining. Nevertheless, AM also poses a few specific constraints and features which
call for innovation in topology optimization [4, 5].
The focus of most research works in topology optimization targeting AM is on manufacturability constraints. Building upon density-based topology optimization, Langelaar [6], Qian [7],
Email address: j.wu-1@tudelft.nl (Jun Wu)
Preprint submitted to Elsevier

December 14, 2017

Figure 1: Illustration of uniform infill that is commonly used in 3D printing (left) and the proposed adaptively refined
infill structures (right). Under the illustrated boundary conditions, the adaptive infill is four times stiffer than the uniform
infill; The compliances of the uniform and adaptive infill are 566.2 and 143.4, respectively. The two designs use the same
amount of material volume. Please refer to the attached GIF animation for the optimization sequence.

and Gaynor and Guest [8] proposed methods to ensure that optimized structures satisfy a constraint on the overhang angle, and thus the usage of additional support structures is reduced or
completely eliminated. The overhang constraint has also been addressed in level set-based topology optimization [9, 10], and topology optimization using moving morphable components [11]
or using a virtual skeleton [12]. Of interest to AM is also the constraint on length scale [13, 14]. A
thorough examination of minimum and maximum length scale is presented in [15]. Among others, Zhou et al. [16] presented geometric constraints to ensure the minimum length scale, while
Zhang et al. [17] developed the minimum length scale based on moving morphable components.
Lazarov and Wang [18] proposed methods to impose maximum length scale in density-based
topology optimization.
Closely related to our current work is research considering infill structures in topology optimization for AM. Rather than printing solid structures, engineering practices in additive manufacturing seem to favour porous structures [19, 20]. For instance, in fused deposition modelling
(FDM), the interior of 3D models is often filled with repetitive triangular or rectangular grids
known as infill. See Figure 1 (left) for an illustration of uniform, rectangular infill pattern within
a 2D bracket-shaped domain. The concept of porous infill is initially introduced to roughly balance cost and mechanical properties: A larger infill volume percentage leads to a stronger print
at the cost of more material and longer printing time. Recent research found out that infill structures can obtain significantly increased stability with respect to buckling [21] and unpredicted
loading conditions [22] at the expense of a minor decrease in stiffness. However, these works
have assumed a constant infill volume percentage across the design domain: In [21, 23] the base
material is interpreted as uniform, isotropic infill, and in [22, 24] anisotropic infill is generated
yet the effective local volume percentage is constant.
1.2. Overview of Adaptive Infill Optimization
In this paper we propose a topology optimization method to generate adaptively refined infill
structures. An optimized quadtree for the bracket shape is shown in Figure 1 (right). Adaptively
refined structures are a good option for 3D printed infill. On the one hand, the adaptivity can
be exploited to enhance regions where denser structures are beneficial to support the mechanical
loads. For instance, in the bracket example the optimized quadtree infill performs four times
stiffer than the uniform infill. On the other hand, the structures on a uniform coarse grid filling
2

the part’s interior make the printed part robust for uncertain and smaller loads, such as forces
during transportation and manipulation that possibly happen to any surface region of the printed
model.
Doing topology optimization by adaptive structure refinement allows for an intuitive control
over some structural features. For instance, the coarsest and finest grids effectively determine
the maximum and minimum void sizes, respectively. The interior structures have a uniform
thickness, which simplifies tool-path generation in 3D printing. The uniform thickness and the
minimum void size have positive implications on thermal effects in the additive manufacturing
process [19]. Going beyond additive manufacturing, structures with a certain level of regularity
facilitate the production of large engineering structures by assembling a few basic sub-structures.
In this case, the final quadtree can be interpreted as being assembled from frames. The frames
have only a few variations in size, and thus each variation can be massively produced.
In topology optimization of adaptive structures we consider structural refinement (i.e., to refine or not to refine) as design variables. This type of design variables was suggested in [25],
where rhombic cells to selectively refined based on a heuristic criterion. Starting with a coarse
grid, it refines the grid cells in a greedy manner: In each iteration, the cells with a high strain
energy density are subdivided. As cells are refined new edges are added into the structure. This
leads to an increase in the structural stiffness. The process is terminated once the material volume
reaches the prescribed budget. This refinement strategy can be enhanced by combing selective
coarsening, which revokes cells that are low in strain energy density. However, accurately solving discrete optimization problem is generally difficult. The greedy approach results in a local
optimal solution, and it might be away from the global optimum.
In this work we address the adaptive structure refinement problem using gradient-based continuous optimization. To do so, we solve two challenging issues. First, for each cell, to refine
it or not to refine it essentially represents a binary design variable. To enable efficient continuous optimization, we relax the binary design variables to continuous values between 0 - not to
refine and 1 - to refine. Following projection approaches in density-based topology optimization, a Heaviside projection [13, 26] is embedded within the optimization in order to improve the
convergence of obtaining discrete 0-1 values.
The continuous representation of refinements also gives rise to a solution to the second challenge — Recursive refinement implies that the design variable of refinement for a fine cell is
not activated until its parent coarse cell has been refined. A consequence is that during the optimization iterations the number of design variables is not constant: The coarse cells that have
already been subdivided are removed from the list of subdivision candidates, and meanwhile
some new subdivision candidates corresponding to newly created fine cells come into action.
Using the continuous design variables, we propose a method to simultaneously take into account
the refinements on all levels. The dependency of the refinement of a fine cell on refinements of
cells on successively coarser levels is encoded by a novel operator, which we call a refinement
filter. Consequently, this results in a constant number of design variables in the whole optimization process, and thus facilitates the implementation. More importantly, considering the design
variables on all levels altogether enlarges the solution space, compared with the greedy approach
where in each iteration only a part of the design variables is activated.
It shall be noted that obtaining spatially adaptive structures clearly distinguishes from previous use of adaptive mesh refinement in topology optimization, where the motivation is to reduce the intensive finite element analysis by reducing the number of elements. Spatial adaptive
meshes (i.e., quadtree in 2D and octree in 3D) are often employed in numerical analysis to attain
a required accuracy for a minimum amount of computation [27]. In the context of topology op3

timization, Maute and Ramm [28] are among the first to employ adaptive techniques to decrease
the number of design variables and to generate smooth structures. Some recent developments
along this direction include adaptive polygonal elements by Nguyen-Xuan [29], and error control in adaptive meshing by Wang et al. [30] and Lambe and Czekanski [31]. The adaptive mesh
used in the literature represents the finite elements for elasticity analysis, and its refinement criterion is based on estimated errors. The obtained structures are similar, if not identical, to the
structures optimized with uniform fine meshes. In contrast, here the adaptive mesh, its edges
in particular, represents the physical structures to be obtained from optimization. The structural
refinement respects global measures on the structural stiffness and the material volume. The
structures our method is aiming at cover the entire closed-walled design domain (see Figure 1
right).
2. Formulation
This section starts with the discretization of the design domain, including the design variables representing mesh refinements, the finite elements for linear elasticity, and the mapping
from design variables to finite elements (Section 2.1). We then proceed to the mathematical formulation of the optimization problem (Section 2.2). The proposed refinement filter is presented
in Section 2.3, and the consequent sensitivity analysis is given in Section 2.4.
2.1. Hybrid Discretization
We maintain a hybrid discretization. As illustrated in Figure 2, it is composed of a quadtree
grid (a) which encodes the topology information, and a discretization by uniform square elements (b) upon which elasticity analysis is performed. On the left, the red lines form the coarse
grid where refinements start from. The coarse cells (red) are mapped to blocks of square elements in the finite element grid on the right. In this example, each coarse cell corresponds to a
block of 162 elements. Within each block, elements in the first boundary layer are assigned solid
(i.e., a density value of 1.0), indicated by black colour. Adjacent coarse cells share a common
edge, and this edge in the finite element grid has a thickness of two elements — one element on
each side.
For illustrative purpose two coarse cells in the bottom right of the domain are refined, creating
two cross-shaped sub-structures (blue). One quadrant is further refined, indicated by the green
cross. The crosses are mapped to solid elements on the right, with different grey intensities
indicating the different levels of refinement. The mapping from the two-level refinements to finite
elements is highlighted in (c) and (d). Here for clarity the mapping of the coarse grid (red) to
finite elements is not included. Note that there is a gap between the two bold crosses in (d). This
is because the refinement does not affect the elements that have already been assigned a value.
This arrangement ensures the mapping from the density field back to the refinement crosses is
one-to-one.
In this example two refinement levels are applied. A further refinement level would merge the
edges and form fully solid blocks in the finite element grid. While such fully solid blocks cause
no problem for the optimization, here we exclude this situation for the purpose of generating a
uniform thickness inside the design domain. The number of supported refinement levels depends
on the resolution within the coarse cell. For instance, increasing from the current 162 within each
coarse cell (red) to 642 would allow two more refinement levels.
4

i
j

(a)

(b)

(c)

(d)

Figure 2: Illustration of the quadtree grid (left) and the corresponding finite elements (right).

2.1.1. Multi-level Design Variables (xk )
A design variable xi, j ∈ [0, 1] is assigned to each of the cells in the quadtree grid, with
xi, j = 1 (resp. xi, j = 0) indicating a refinement (resp. non-refinement) for the cell with the 2D
indices (i, j). Design variables are assigned on each level in the quadtree, denoted by xi,k j , where
k indicates the k-th refinement level. The initial coarse grid is denoted as the first level. From the
first level downwards, the resolution of cells, which is also the resolution of the design field, is
doubled on each dimension from the previous level. Denoting the resolution on the k-th level by
(nkx , nky ), and the resolution of the initial coarse grid by (n0x , n0y ), the relation can be written as
(nkx , nky ) = (2k−1 n0x , 2k−1 n0y ),

k = 1, ..., k̄,

(1)

where k̄ is the maximum allowed refinement level. For a coarse block of 2m × 2m elements, k̄ is
calculated by
k̄ = m − 2.
(2)
When the refinement level reaches k̄, it creates voids in the size of 22 finite elements.
2.1.2. Mapping (xk → ρ)
The mapping from the quadtree grid to the density field (ρ) has been illustrated in Figure 2.
Solid finite elements can be mapped either from the design variables xk , k = 1, ..., k̄, corresponding to cross-shaped sub-structures, or from the initial coarse cells (red). For compactness in the
following formulation, let us denote the initial coarse cells by x0 which has the same resolution
as x1 , i.e., (n0x , n0y ) = (n1x , n1y ). Note that x0 has a constant value xe0 = 1 for all coarse cells. This
corresponds to a coarse frame for ensuring the minimum stiffness of the optimized structure.
Reshaping the multi-level design fields xk |nkx ×nky , k = 0, ..., k̄, and the density field ρ|2m n0x ×2m n0y
into column vectors xk and ρ, respectively, the mapping can be written compactly as a transformation by
k̄
P
(3)
ρ = T k xk
k=0

5

The transformation T k is a sparse matrix with a dimension of 22m n0x n0y ×nkx nky . The non-zero entries
in T k have a constant value of 1, and concern the specific finite elements that will be assigned by
refinements xi,k j . Since the non-zero entries have a unit value, the transformation assigns the same
value of the refinement to the affected finite elements. For instance, an intermediate refinement
by xi,k j = 0.5 leads to intermediate densities for the corresponding finite elements, ρ = 0.5.
2.2. Optimization Problem
With the design variables (xk ) and the physical density field (ρ) defined in the previous subsection, we can write the commonly-studied compliance minimization problem (e.g., in [32, 33])
as follows:
min c = UT K(ρ)U,
x

s.t.

K(ρ)U = F,
P
ρe ve ≤ V ∗ ,

V(ρ) =

(4)
(5)
(6)

∀e

xi,k j ∈ [0, 1].

(7)

The objective is to minimize the compliance represented by c. The density vector ρ is calculated
from the design variable field xk by Eq. (3). The first constraint is the static equation of linear
elasticity. K, U, and F are stiffness matrix, displacement vector, and force vector, respectively.
The global stiffness matrix K is assembled from element stiffness matrix ke . The modified SIMP
interpolation model [33] is used to compute the element stiffness matrix:
ke = Ee (ρe )k0 ,
Ee (ρe ) = Emin + ρep (E0 − Emin ),

(8)
(9)

where k0 is the element stiffness matrix for an element with unit Young’s modulus, E0 is the
stiffness of the material, Emin is a very small stiffness assigned to void elements in order to prevent
the stiffness matrix from becoming singular, and p is a penalization factor (typically p = 3)
introduced to ensure black-and-white solutions. The second constraint measures the volume
occupied by solid elements. ve is the unit volume for each element, and V ∗ is the maximum
allowed volume. The last constraint restricts the design variables to take a value between 0
(non-refinement) and 1 (refinement).
Solving this optimization problem for the cantilever beam presents the design shown in Figure 3. It generates many crossing sub-structures which correspond to the design variables. However, many of such sub-structures (pointed out by the red arrows) are suspended. This is due to
the fact that the (imaginary) fine cells are refined, even though the fine cells themselves have not
yet been created from refining their parent cells. The suspended sub-structures don’t effectively
carry mechanical loads, and disobey the intention to create quadtree-like mechanical structures.
In the following, a refinement filter is proposed to encode the recursive refinement rules, and thus
to effectively filter out such suspended sub-structures.
2.3. The Refinement Filter (xk → x̃k )
In creating a quadtree the refinement is only applicable to a cell that already exists. This
implies dependency of the refinement of a fine cell on the refinement of its parent cell. In the
binary setting, the refinement of a cell on the k−th level, xi,k j , will not happen, as long as the
6

Figure 3: When the refinement rules are not applied, the optimized structure consists of many suspended structures,
indicated by the red arrows. The compliance value is 136.1.

−1 −1
refinement of its parent cell has not yet taken place, i.e., xik−1
, j ) refers to
−1 , j−1 = 0. Here (i
the indices of the parent cell. To encode this dependency in the setting of continuous design
variables, we introduce a filtered design variable x̃i,k j :

x̃i,k j = min(xi,k j , x̃ik−1
−1 , j−1 ).

(10)

In case of xik−1
−1 , j−1 = 0, i.e., the parent cell is not refined, the filtered value is
x̃i,k j = min(xi,k j , 0) = 0.

(11)

This min function thus restricts the design variable on the fine level when the parent cell has not
been refined. In case of xik−1
−1 , j−1 = 1, the filtered value is unaffected,
x̃i,k j = min(xi,k j , 1) = xi,k j .

(12)

Applying refinement recursively creates a quadtree with multiple levels. The dependency
thus goes from the fine cell through all the intermediate levels to one of the coarse cells on the
3
first level. For instance, consider the design variable x4,1
illustrated in Figure 4. The filtered
3
3
2
2
2
1
value is x̃4,1 = min(x4,1 , x̃2,1 ), where x̃2,1 = min(x2,1 , x̃1,1 ). At first glance this requests multiple
filters executed sequentially, leading to much complex sensitive analysis afterwards. However,
the recursive min functions can be consolidated into an equivalent min function by taking in as
input all the involved arguments in the individual min functions, i.e.,


1
x̃i,k j = min xi,k j , xik−1
(13)
−1 , j−1 , ..., xi−k+1 , j−k+1 .
This consolidation does not introduce approximate error, and significantly simplifies sensitivity
analysis.
After the design variables have been filtered, the density vector, ρ, is updated. This is done
by replacing design vector xk by the filtered design vector x̃k in Eq. (3), i.e.,
ρ=

k̄
P

T k x̃k .

k=0

7

(14)

i
j
Level 1
1
x1,1

Level 2
2
x1,1

2
x2,1

2
x1,2

2
x2,2

3
x3,1

3
x4,1

3
x3,2

3
x4,2

Level 3

3 depends on x2 which depends on x1 .
Figure 4: Illustration of dependency in recursive refinement. For instance, x4,1
2,1
1,1

Figure 5: Left: A randomly initialized design field shows many suspended structures. Right: The suspended structures
disappear when the refinement filter is applied.

The effect of the recursive refinement filter is visualized in Figure 5. On the left, the multilevel design variables are initialized randomly with a value of 0 or 1. A few suspended substructures can be observed on the left, and they disappear on the right. Figure 6 shows the result
of integrating this filter into the optimization for the cantilever beam. The effect can be observed
by comparing to Figure 3 which contains many suspended sub-structures.
2.3.1. Smooth Approximation
Gradient-based numerical optimization necessitates differentiable functions. To this end, the
non-differentiable min function (Eq. (13)) is approximated by a smooth p-norm function, with a

8

Figure 6: A quadtree structure obtained from optimizing a cantilever beam. The compliance value is 143.3.

negative exponent,


1
x̃i,k j = min xi,k j , xik−1
−1 , j−1 , ..., xi−k+1 , j−k+1
!1
P  k−l  pn pn
1 k−1
≈
.
x
k l=0 i−l , j−l

(15)
(16)

Here pn represents the p-norm exponent, to be distinguished from the penalty p in SIMP interpolation (Eq. (9)). As pn approaches negative infinity, the p-norm equals the minimum value. The
p-norm is normalized by a factor considering the number of arguments, 1k in this formulation.
The normalization reduces the approximation error when a practical pn value is not infinite [22].
In our tests, we use pn = −16.
The indices of parent cells are needed in Eq. (16). When using a regular discretization, the
indices can be found by a recursive function,
i−l = h(i, −l) = h(h(i, −l + 1), −1),

(17)

where the function h is defined as
!
1
h(i, −1) = floor (i + 1) .
2

(18)

2.4. Sensitivity Analysis
By using the chain rule, the sensitivity of the objective function c with respect to the design
variable xk is given by:
!
k̄
P
∂c
∂c ∂ρ ∂ x̃l
=
.
(19)
∂xk l=1 ∂ρ ∂ x̃l ∂xk
The first term ∂c/∂ρ is derived from adjoint analysis. The individual entry is
∂c
= −p(ρe ) p−1 (E0 − Emin )uTe k0 ue .
∂ρe

(20)

The second term can be evaluated from Eq. (14),
∂ρ
= Tl .
∂ x̃l
9

(21)

The third term comes from the smoothed refinement filter, Eq. (16). The non-zero entries
(i.e., if l ≥ k) are calculated by
∂ x̃i,l j
∂xikk−l , jk−l

P  l−m  pn
1 1 l−1
=
x −m −m
l l m=0 i , j

! p1 −1 
n

xikk−l , jk−l

 pn −1

.

(22)

The sensitivity of the volume V with respect to the design variable xk is derived in a similar
way:
!
k̄ ∂V ∂ρ ∂ x̃l
P
∂V
=
.
(23)
∂xk l=k ∂ρ ∂ x̃l ∂xk
Individual entries in the first term are given by
∂V
= ve ,
∂ρe

(24)

where ve is the unit volume for each element.
3. Extensions
3.1. Balanced Quadtree
We provide an option to restrict the quadtree such that the level difference between neighbouring cells is at most one, known as balanced or restricted quadtree. A balanced quadtree
maintains a smooth transition from fine cells to coarse cells. For 3D printed infill, a smooth
transition is expected to make the part more stable for uncertain loads.
In an unbalanced quadtree, the refinement of a cell depends only on its parent cell (Eq. (10)).
In the balanced quadtree, the refinement additionally depends on the neighbours of its parent
cell. The refinement filter (Eq. (10)) is thus updated by
k−1
x̃i,k j = min(xi,k j , x̃ik−1
−1 , j−1 , x̃i−1 ±1, j−1 ±1 ).

(25)

−1
k−1
Here, it is assumed that the indices i−1 ± 1 belong to [1, nk−1
x ] and j ± 1 belong to [1, ny ], otherwise the invalid arguments are excluded. A consolidation of this recursive function is obtained
similar to Eq. (13). The recursive function can point to the same coarser cell multiple times,
due to the introduction of neighbours. Only one copy of each dependent cell is included in the
consolidated filter. Since the encoding of neighbourhood and hierarchy remains constant during
the optimization, finding the dependency is only performed once in the initialization.
The effect of the balanced refinement filter is visualized in Figure 7, by comparing the unbalanced refinement (left) and the balanced refinement (right). An optimized balanced quadtree for
the cantilever beam is shown in Figure 8.

3.2. Integration with other Filters
In density-based topology optimization using square elements, regulations such as density
filter [13] or sensitivity filter [32] are often applied to get rid of checkerboard patterns (i.e.,
regions of alternating black and white elements). In our formulation, the refinement itself can be
regarded as a regulation, and an additional density/sensitivity filter is no longer necessary.
In topology optimization with continuous design variables, projection filters are often applied
to improve the convergence of approaching a binary white-or-black design. Our continuous
10

Figure 7: Left: An unbalanced quadtree. Right: A balanced quadtree obtained by filtering out fine cells which are at least
two levels smaller than one of its neighbours.

Figure 8: A balanced quadtree structure obtained from optimizing a cantilever beam. The mechanical structure is
adaptively refined for achieving the highest stiffness for the illustrated load, under the constraint that 40% of rectangular
design domain is occupied by the structure. The compliance value is 165.7.

reformulation of refinement also benefits from such a projection. In particular, we integrate the
projection filter proposed in [26]:

 
tanh (βη) + tanh β xi,k j − η
k
.
(26)
x̄i, j =
tanh (βη) + tanh (β (1 − η))
Here β controls the sharpness of projection. In the limit of β → ∞ the projection approaches
a discontinuous step function. We employ a continuation of this parameter. η ∈ [0, 1] is the
threshold value, and we choose η = 0.5.
4. Numerical Examples
The method has been implemented based on the Matlab code provided by Andreassen et
al. [33]. The optimization problem is solved by using the method of moving asymptotes (MMA) [34].
In the following we report numerical results tested on regular design domains (Section 4.1)
and curved design domains (section 4.2), and compare results from continuous and discrete optimization (Section 4.3). The following parameters are used in all tests: The SIMP penalization
p = 3, the p-norm pn = −16, the threshold for projection η = 0.5. The sharpness of projection β
starts with 1 and is increased up to 32, by doubling its value every 60 iterations. The maximum
number of iterations is 400.
11

4.1. Regular Domain
4.1.1. Cantilever Beam

(a) It. 40, c = 317.2, s = 0.631

(b) It. 100, c = 249.8, s = 0.500

(c) It. 160, c = 215.4, s = 0.362

(d) It. 220, c = 177.9, s = 0.151

(e) It. 280, c = 169.2, s = 0.070

(f) It. 360, c = 165.7, s = 0.040

Figure 9: A sequence of the balanced quadtree during the optimization process. The number of iterations (It.), the
compliance value (c), and the sharpness (s) are reported. Please refer to the attached GIF animation for the entire
optimization sequence.
0.5

1

350
300
250

0.8

0.45

Sharpness

Volume fraction

Compliance

400

0.4

0.6
0.4

0.35
0.2

200

0.3
0

50

100

150

200

Iteration

250

300

350

0

0

50

100

150

200

Iteration

250

300

350

0

50

100

150

200

250

300

350

Iteration

Figure 10: The convergence plots for the cantilever beam design (Figure 9). From left to right: Convergence plot of
compliance, of volume fraction, and of sharpness. The jumps in curves are due to the β-continuation in the Heaviside
projection.

The first example is a cantilever beam that has appeared multiple times during the exposition
of the method. The boundary condition is illustrated in Figure 3. The beam is fixed on the left
side, while a force is applied on the right side in the middle. The prescribed material volume is
40% of the design domain. The rectangular design domain is initialized with a coarse grid of
12

8 × 4. Each grid cell contains 64 × 64 square elements, leading to a finite element resolution of
512 × 256 for the entire domain. The maximum refinement level is 4.
Figure 9 shows a sequence of intermediate structures during the optimization of a balanced
quadtree structure. As the optimization progresses, the quadtree structure emerges from the grey
density distribution. To measure how close the continuous density field is to a binary field, the
sharpness factor [35] is defined as
s=

4X
(ρe (1 − ρe )),
n e

(27)

where n is the number of finite elements. The sharpness factor becomes 0.0 when the density
values are converged to a strict 0-1 solution, and it becomes 1.0 if all elements take a value
of 0.5. The optimized quadtree (Figure 9 bottom right) has a sharpness factor of 0.04, which
corresponds the field where density values on average take 0.01 or 0.99.
The convergence plots are shown in Figure 10. The plot on the left shows that the compliance
value decreases almost monotonically, except for a few jumps corresponding to the increase of
β in the Heaviside projection. The middle shows the volume fraction during the process. The
design variables are initialized homogeneously with a value of 0.49 such that the volume fraction
of 0.4 has been satisfied from the beginning. The plot on the right shows the convergence of
sharpness. The optimization is terminated at the iteration of 379 as the maximum change in the
density field becomes smaller than 1e−4.
4.1.2. MBB Beam
The second example is the MBB design problem shown in Figure 11. Due to symmetry only
one half of the design domain is optimized. The prescribed material volume is 30% of the design
domain. The half domain is initialized with a coarse grid of 12 × 4, where each grid cell is
composed of 64 × 64 elements, leading to a total of 768 × 256.

L

6L
Figure 11: Illustration of the MBB beam design problem.

Figure 12 compares the quadtree structures obtained by using different refinement levels. The
three rows from top to bottom show optimized quadtree using a maximum refinement level of
three, four, and five, respectively. In each row, the left corresponds to the unbalanced refinement,
while the right is obtained in the setting of balanced refinement. Comparing left and right, it
can be found that the balanced quadtree (right) has a larger compliance and a larger sharpness
value than in the unbalanced setting (left). This is due to the fact that the balanced refinement
additionally involves the neighbours of the parent cells. Comparing the three rows, we observe
smaller compliances when a larger refinement level is applied. This is because a larger refinement
level increases the solution space. At a refinement level of five (bottom), small cavities from the
13

(a) Unbalanced quadtree, k̄ = 3, c = 1083.5, s = 0.002

(b) Balanced quadtree, k̄ = 3, c = 1127.9, s = 0.008

(c) Unbalanced quadtree, k̄ = 4, c = 571.6, s = 0.010

(d) Balanced quadtree, k̄ = 4, c = 711.5, s = 0.028

(e) Unbalanced quadtree, k̄ = 5, c = 449.0, s = 0.030

(f) Balanced quadtree, k̄ = 5, c = 586.1, s = 0.070

Figure 12: A comparison of optimized structures using a refinement level of, from top to bottom, three, four, and five.
The unbalanced and balanced quatrees are shown on the left and right, respectively.
0.5

1500

Volume fraction

Compliance

2000

1000

0.4
0.35

0.25

0

0.2
50

100

150

200

250

300

350

Iteration

400

level 3, balanced
level 3, unbalanced
level 4, balanced
level 4, unbalanced
level 5, balanced
level 5, unbalanced

0.8

0.3

500

0

1
level 3, balanced
level 3, unbalanced
level 4, balanced
level 4, unbalanced
level 5, balanced
level 5, unbalanced

0.45

Sharpness

level 3, balanced
level 3, unbalanced
level 4, balanced
level 4, unbalanced
level 5, balanced
level 5, unbalanced

2500

0.6
0.4
0.2
0

0

50

100

150

200

Iteration

250

300

350

400

0

50

100

150

200

250

300

350

400

Iteration

Figure 13: The convergence of compliance (left), volume fraction (middle), and sharpness (right) for the MBB designs
shown in Figure 12. The jumps in curves are due to the β-continuation in the Heaviside projection.

previous refinement (middle) are filled with material. The optimizations with these different
settings all converge well, as can be seen from the convergence plots in Figure 13.
4.2. Curved Design Domain
The adaptive structure refinement works as well on curved design domains. For a curved
design domain, for simplicity in implementation we simulate the smallest axis-aligned bounding
box enclosing the curved domain. The voids outside the curved domain are assigned as passive
void elements, while a boundary layer with a thickness of 2 elements are passive solid elements.
This requires updating the corresponding entries in the transformation matrix in Eq. (14). In the
following we show two examples.

14

r
2r

r
r
4r

Figure 14: Illustration of the bracket design domain and boundary conditions.

(a) c = 3789.0

(b) c = 1002.3

Figure 15: The adaptively refined infill is 3.8 times stiffer than the uniform infill.

4.2.1. Bracket
Figure 14 illustrates the dimension and boundary conditions of the bracket design problem.
It contains curved boundary where a quadtree mesh is not conformal. The tight bounding box of
the bracket is discretized using a coarse grid of 12 × 8, where each grid cell has 32 × 32 elements,
i.e., a total of 384 × 256 elements for the bounding box. First we uniformly refine the coarse
grid by two levels. The uniform infill is shown in Figure 1 (left). Its volume fraction is 33.02%.
Using this volume fraction as a constraint we perform adaptive structure optimization, allowing
three levels of refinement. The optimized quadtree infill is shown in Figure 1 (right). Numerical
analysis suggests that the adaptive infill is four times stiffer than the uniform infill.
Figure 15 further compares the structural performance using a different volume fraction. With
one refinement level the uniform infill takes a volume fraction of 18.53% and results in a compliance value of 3789.0. Using the same volume fraction and allowing two levels of refinement
the optimized infill has a compliance value of 1002.3.
4.2.2. Curved Beam
The second curved design domain is shown in Figure 16 (left). On top of the curved beam
there are two loads. The design domain has a resolution of 8 × 8 initial cells, each of which has
64 × 64 finite elements, leading to a finite element grid of 512 × 512. The reference uniform infill
shown in Figure 16 (right) is obtained by refining the initial cells with two refinement levels. The
volume fraction is 14.66%.
Figure 17 shows the optimized balanced quadtree under different force specifications. A
refinement level of three is applied in the optimization. The improved stiffness of using adaptive
15

r

2r

Figure 16: Left: Illustration of the bracket design domain and boundary conditions. Right: A uniformly refined infill
within the curved design domain.

Figure 17: The optimized infill under different force specification. On the right, the two forces are applied independently.

refinement is also confirmed with this example. From left to right, the compliances are 141.7,
611.9, and 420.9 (the averaged value considering both loads), respectively. As a reference, the
compliances of the uniform infill under the respective boundary conditions are 585.1, 1448.4,
and 1016.8, respectively.
4.3. Comparison to Greedy Approach
The quadtree structures optimized by the continuous formulation are compared with those
obtained by the greedy approach [25]. The greedy approach starts refining from the initial coarse
grid, and performs selective refinement at odd iterations and selective coarsening at even iterations. In each odd iteration, following finite element analysis, the leaf cells in the quadtree are
sorted by their strain energy density. The first few cells in this sorted list are then successively
refined, until the volume increase exceeds a small volume change, 0.4% in our examples. In
each even iteration, following finite element analysis, blocks of 22 leave cells, i.e., the cells on
the second finest level, are sorted by the averaged strain energy density. The last few cells in
this sorted list are coarsened, until the volume is decreased by 0.1%. The combination of refinement and coarsening is found more optimal than purely refinement, and thus serves a reference
for comparison. For simplicity the greedy approach is only applied in the setting of unbalanced
quadtree.
16

v f = 0.20, c = 1919.9

v f = 0.20, c = 1266.7

v f = 0.30, c = 742.5

v f = 0.30, c = 571.6

v f = 0.40, c = 449.6

v f = 0.40, c = 402.8

v f = 0.50, c = 328.2

v f = 0.50, c = 312.0

Figure 18: Quadtree structures obtained by a greedy approach in the binary formulation (left) and optimized in the
continuous formulation (right). From top to bottom the comparisons are made under increasing volume fractions.

Figure 18 show the optimized MBB structures by the greedy approach (left) and by the continuous formulation (right). The structures are compared under four different volume fractions.
The structures on the right form some distinct dense and sparse regions, to some extent similar to
the standard topology optimization, while structures on the left have regions of alternating refined
and non-refined cells. The quadtree structures optimized by the continuous formulation consistently have smaller compliances. The difference in compliance values is especially significant
when the volume fraction is small (cf. the top row).
5. Conclusion
We have presented a topology optimization method to design adaptively refined structures.
The structures are refined within a closed-walled design domain. The adaptive structure offers
highly optimized stiffness for specified mechanical loads, and outperforms the commonly used
uniform infill pattern in 3D printing. Meanwhile the structure covers the interior of the mechanical part with controllable and smoothly varying void sizes, making it stable for uncertain and
small perturbation forces. It is expected that the design of optimized and adaptive infill structures
will be integrated in slicing software.
17

Core to our method is a continuous representation of the binary topology associated with
edges in adaptive meshing. Relaxing binary design variables to continuous values allows to
derive sensitivity for numerical optimization, and enables a novel refinement filter to encode the
dependency in recursive refinements. Numerical results suggest that the continuous formulation
finds solutions that are more optimal compared with a greedy approach for solving the discrete
formulation. The method is easy to implement. While the refinement filter is demonstrated on
2D Cartesian grids, its applicability to 3D is straightforward. We leave this as future work.
References
[1] M. P. Bendsøe, O. Sigmund, Topology optimization: theory, methods, and applications, Springer-Verlag Berlin
Heidelberg, 2004.
[2] O. Sigmund, K. Maute, Topology optimization approaches, Struct. Multidiscip. Optim. 48 (6) (2013) 1031–1055.
doi:10.1007/s00158-013-0978-6.
[3] J. D. Deaton, R. V. Grandhi, A survey of structural and multidisciplinary continuum topology optimization: post
2000, Struct. Multidiscip. Optim. 49 (1) (2014) 1–38. doi:10.1007/s00158-013-0956-z.
[4] D. Brackett, I. Ashcroft, R. Hague, Topology optimization for additive manufacturing, in: Proceedings of the solid
freeform fabrication symposium, Austin, TX, Vol. 1, S, 2011, pp. 348–362.
[5] T. Zegard, G. H. Paulino, Bridging topology optimization and additive manufacturing, Structural and Multidisciplinary Optimization 53 (1) (2016) 175–192. doi:10.1007/s00158-015-1274-4.
[6] M. Langelaar, An additive manufacturing filter for topology optimization of print-ready designs, Struct. Multidiscip. Optim. (2016) 1–13doi:10.1007/s00158-016-1522-2.
[7] X. Qian, Undercut and overhang angle control in topology optimization: A density gradient based integral
approach, International Journal for Numerical Methods in Engineering 111 (3) (2017) 247–272, nme.5461.
doi:10.1002/nme.5461.
[8] A. T. Gaynor, J. K. Guest, Topology optimization considering overhang constraints: Eliminating sacrificial support
material in additive manufacturing through design, Struct. Multidiscip. Optim. 54 (5) (2016) 1157–1172. doi:
10.1007/s00158-016-1551-x.
[9] A. M. Mirzendehdel, K. Suresh, Support structure constrained topology optimization for additive manufacturing,
Comput. Aided Des. 81 (C) (2016) 1–13. doi:10.1016/j.cad.2016.08.006.
[10] G. Allaire, C. Dapogny, R. Estevez, A. Faure, G. Michailidis, Structural optimization under overhang constraints
imposed by additive manufacturing technologies (2017) 1–41Working paper or preprint.
[11] X. Guo, J. Zhou, W. Zhang, Z. Du, C. Liu, Y. Liu, Self-supporting structure design in additive manufacturing
through explicit topology optimization, Computer Methods in Applied Mechanics and Engineering 323 (2017) 27
– 63. doi:10.1016/j.cma.2017.05.003.
[12] Y. Mass, O. Amir, Topology optimization for additive manufacturing: accounting for overhang limitations using a
virtual skeleton, in: Proceedings of ECCOMAS, 2016.
[13] J. K. Guest, J. Prévost, T. Belytschko, Achieving minimum length scale in topology optimization using nodal design
variables and projection functions, International Journal for Numerical Methods in Engineering 61 (2) (2004) 238–
254. doi:10.1002/nme.1064.
[14] J. Guest, Imposing maximum length scale in topology optimization, Struct. Multidiscip. Optim. 37 (5) (2009)
463–473. doi:10.1007/s00158-008-0250-7.
[15] B. S. Lazarov, F. Wang, O. Sigmund, Length scale and manufacturability in density-based topology optimization,
Archive of Applied Mechanics 86 (1) (2016) 189–218. doi:10.1007/s00419-015-1106-4.
[16] M. Zhou, B. S. Lazarov, F. Wang, O. Sigmund, Minimum length scale in topology optimization by geometric
constraints, Computer Methods in Applied Mechanics and Engineering 293 (2015) 266 – 282. doi:10.1016/j.
cma.2015.05.003.
[17] W. Zhang, D. Li, J. Zhang, X. Guo, Minimum length scale control in structural topology optimization based on the
moving morphable components (mmc) approach, Computer Methods in Applied Mechanics and Engineering 311
(2016) 327 – 355. doi:10.1016/j.cma.2016.08.022.
[18] B. S. Lazarov, F. Wang, Maximum length scale in density based topology optimization, Computer Methods in
Applied Mechanics and Engineering 318 (2017) 826 – 844. doi:10.1016/j.cma.2017.02.018.
[19] I. Gibson, D. W. Rosen, B. Stucker, Design for Additive Manufacturing, Springer US, Boston, MA, 2010, pp.
299–332. doi:10.1007/978-1-4419-1120-9_11.
[20] W. Gao, Y. Zhang, D. Ramanujan, K. Ramani, Y. Chen, C. B. Williams, C. C. Wang, Y. C. Shin, S. Zhang, P. D.
Zavattieri, The status, challenges, and future of additive manufacturing in engineering, Computer-Aided Design 69
(2015) 65–89. doi:10.1016/j.cad.2015.04.001.

18

[21] A. Clausen, N. Aage, O. Sigmund, Exploiting additive manufacturing infill in topology optimization for improved
buckling load, Engineering 2 (2) (2016) 250 – 257. doi:10.1016/J.ENG.2016.02.006.
[22] J. Wu, N. Aage, R. Westermann, O. Sigmund, Infill optimization for additive manufacturing – approaching bonelike porous structures, IEEE Transactions on Visualization and Computer Graphics PP (99) (2017) 1–14. doi:
10.1109/TVCG.2017.2655523.
[23] A. Clausen, N. Aage, O. Sigmund, Topology optimization of coated structures and material interface problems,
Computer Methods in Applied Mechanics and Engineering 290 (2015) 524 – 541. doi:10.1016/j.cma.2015.
02.011.
[24] J. Wu, A. Clausen, O. Sigmund, Minimum compliance topology optimization of shell-infill composites for additive
manufacturing, Computer Methods in Applied Mechanics and Engineering (2017) –doi:10.1016/j.cma.2017.
08.018.
[25] J. Wu, C. C. Wang, X. Zhang, R. Westermann, Self-supporting rhombic infill structures for additive manufacturing,
Computer-Aided Design 80 (2016) 32 – 42. doi:10.1016/j.cad.2016.07.006.
[26] F. Wang, B. S. Lazarov, O. Sigmund, On projection methods, convergence and robust formulations in topology optimization, Structural and Multidisciplinary Optimization 43 (6) (2010) 767–784. doi:10.1007/
s00158-010-0602-y.
[27] M. J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of computational Physics 53 (3) (1984) 484–512.
[28] K. Maute, E. Ramm, Adaptive topology optimization, Structural optimization 10 (2) (1995) 100–112. doi:10.
1007/BF01743537.
[29] H. Nguyen-Xuan, A polytree-based adaptive polygonal finite element method for topology optimization, International Journal for Numerical Methods in Engineering 110 (10) (2017) 972–1000, nme.5448. doi:10.1002/nme.
5448.
[30] Y. Wang, Z. Kang, Q. He, Adaptive topology optimization with independent error control for separated displacement and density fields, Computers & Structures 135 (2014) 50 – 61. doi:10.1016/j.compstruc.2014.01.
008.
[31] A. B. Lambe, A. Czekanski, Topology optimization using a continuous density field and adaptive mesh refinement,
International Journal for Numerical Methods in Engineering n/a–n/aNme.5617. doi:10.1002/nme.5617.
[32] O. Sigmund, A 99 line topology optimization code written in matlab, Struct. Multidiscip. Optim. 21 (2) (2001)
120–127. doi:10.1007/s001580050176.
[33] E. Andreassen, A. Clausen, M. Schevenels, B. S. Lazarov, O. Sigmund, Efficient topology optimization in matlab using 88 lines of code, Structural and Multidisciplinary Optimization 43 (1) (2010) 1–16. doi:10.1007/
s00158-010-0594-7.
[34] K. Svanberg, The method of moving asymptotes—a new method for structural optimization, International Journal
for Numerical Methods in Engineering 24 (2) (1987) 359–373. doi:10.1002/nme.1620240207.
[35] O. Sigmund, Morphology-based black and white filters for topology optimization, Structural and Multidisciplinary
Optimization 33 (4) (2007) 401–424. doi:10.1007/s00158-006-0087-x.

19

