A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D

Abstract

Compact and efficient Matlab implementations of compliance topology optimization (TO) for 2D and 3D continua are given, consisting of 99 and 125 lines respectively. On discretizations ranging from 3 ⋅ 10⁴ to 4.8 ⋅ 10⁵ elements, the 2D version, named top99neo, shows speedups from 2.55 to 5.5 times compared to the well-known top88 code of Andreassen et al. (Struct Multidiscip Optim 43(1):1–16, 2011). The 3D version, named top3D125, is the most compact and efficient Matlab implementation for 3D TO to date, showing a speedup of 1.9 times compared to the code of Amir et al. (Struct Multidiscip Optim 49(5):815–829, 2014), on a discretization with 2.2 ⋅ 10⁵ elements. For both codes, improvements are due to much more efficient procedures for the assembly and implementation of filters and shortcuts in the design update step. The use of an acceleration strategy, yielding major cuts in the overall computational time, is also discussed, stressing its easy integration within the basic codes.

Appendices

Appendix A: Elaboration on the OC update

Let us consider (3) at a given design point x_k assuming the reciprocal and linear approximation for the compliance and volume functions, respectively (Christensen and Klarbring 2008)

$$ \begin{cases} & \min\limits_{\mathbf{x}\in [\delta_{-}, \delta_{+}]^{m}} c\left( \mathbf{x} \right) \simeq c_{k} + {\sum}^{m}_{e=1} (-x^{2}_{k, e}\partial_{e} c(\mathbf{x}_{k})) x^{-1}_{e} \\ \text{s.t.} & {\sum}^{m}_{e=1} \partial_{e} V(\mathbf{x}_{k}) x_{e} - f|{\Omega}_{h}| \leq 0 \end{cases} $$

(24)

We set up the Lagrangian associated with (24)

$$ L(\mathbf{x}, \lambda ) = c(\mathbf{x}) + \lambda \left( \sum\limits^{m}_{e=1} \partial_{e} V(\mathbf{x}_{k}) x_{e} - f|{\Omega}_{h}| \right) $$

and seek the pair $(\mathbf {x}_{k+1},\lambda ^{\ast }_{k})\in \mathbb {R}^{m}\times \mathbb {R}_{+}$ solving the subproblem

$$ \max_{\lambda > 0} \left\{ \psi(\lambda) := \min\limits_{\mathbf{x}\in \mathcal{C}} L(\mathbf{x}, \lambda )\right\} $$

(25)

where $\mathcal {C} = \{ \mathbf {x} \in \mathbb {R}^{m} \mid \delta _{-}\leq x_{e} \leq \delta _{+}, e = 1 \ldots , m \}$ and ψ(λ) is the dual function. Equation (25) is solved by primal-dual (PD) iterations, as x and λ are interlaced. Replacing ξ = x_k and using subscripts (j) to denote inner PD iterations, we have

• 1. Fixed λ = λ_(j), the inner minimization in (25) gives

  $$ {\xi^{2}_{e}}\partial_{e}c(\boldsymbol{\xi})x^{-2}_{e} + \lambda \partial_{e}V(\boldsymbol{\xi}) = 0 \Longrightarrow x_{e} = \xi_{e} \left( -\frac{\partial_{e}c(\boldsymbol{\xi})} {\lambda\partial_{e}V(\boldsymbol{\xi})} \right) ^{\frac{1}{2}} $$

  due to separability of the approximation. Let us denote the rightmost expression $x_{e} = \mathcal {F}_{(j)e}(\lambda )$, and taking into account the box constraints in $\mathcal {C}$, we have

  $$ \mathcal{U}(x_{e}) = \begin{cases} x_{(j+1),e} = \delta_{-} & \text{if} e \in \mathcal{L} = \{e\mid x_{(j+1),e} \leq \delta_{-} \} \\ x_{(j+1),e} = \delta_{+} & \text{if} e \in \mathcal{U} = \{e\mid x_{(j+1),e} \geq \delta_{+} \} \\ x_{(j+1),e} = \mathcal{F}_{(j),e} & \text{if} e \in \mathcal{M} = \{e\mid\delta_{-} < x_{(j+1),e} < \delta_{+} \} \\ \end{cases} $$

  (26)

  where $\mathcal {C} = {\mathscr{L}} + \mathcal {U} + {\mathscr{M}}$. The above is equivalent to (10).

• 2. We then evaluate the dual function for x_(j+ 1) given by (26), and the stationarity (∂_λψ = 0) gives

  $$ \sum\limits_{e=1}^{m} \partial_{e} V(\boldsymbol{\xi}) (\chi_{\mathcal{U}}\delta_{+} + \chi_{\mathcal{L}}\delta_{-} + \mathcal{F}_{(j),e}(\lambda)\chi_{\mathcal{M}} ) - f|{\Omega}_{h}| = 0 $$

  where χ_[⋅] is the characteristic function of a set. In this simple case, the above can be solved for λ_(j+ 1), the Lagrange multiplier enforcing the volume constraint for the updated density x_(j+ 1), and after some simplifications, we obtain

  $$ \lambda_{(j+1)} = \left( \frac{{\sum}_{e\in\mathcal{M}}x_{(j+1)e}(\partial_{e}c(\boldsymbol{\xi})/\partial_{e}V(\boldsymbol{\xi}))^{1/2}} {f|{\Omega}_{h}|/\partial_{e}V(\boldsymbol{\xi}) - |\mathcal{L}|\delta_{-} - |\mathcal{U}|\delta_{+}}\right)^{2} $$

  (27)

  where |⋅| denotes the number of elements in a set.

Equations (26) and (27) can be iteratively used to compute the new solution $(\mathbf {x}_{k+1},\lambda ^{\ast }_{k})$, as implemented in the code here below (again, note that lm here represents $\sqrt {\lambda }$)

and, for the MBB beam example, this performs as shown by the green curves in Fig. 4b.

However, a closed form expression such as (27) cannot be obtained for more involved constraint expressions and therefore a root finding strategy must be employed to approximate the Lagrange multiplier. The application of (27) to the current, feasible design point (x_(j+ 1) = x_k) reduces to

$$ \lambda^{\#} = \left[ \frac{1}{mf} \sum\limits^{m}_{e=1} x_{k,e} \left( -\frac{\partial_{e}c(\boldsymbol{\xi})} {\partial_{e}V(\boldsymbol{\xi})} \right)^{1/2} \right]^{2} $$

(28)

since $|{\mathscr{M}}| = |{\Omega }_{h}| = m$, $|{\mathscr{L}}| = |\mathcal {U}| = 0$ and we made use of (7). We immediately verify that (28) is identical to (19).

Appendix B: The 2D code for compliance minimization

Appendix C: 3D code for compliance minimization