Source-linked AI summary
A new generation 99 line Matlab code for compliance Topology Optimization and its extension to 3D
Federico Ferrari, Ole Sigmund
TL;DR
The paper addresses the need for more efficient Matlab implementations of compliance topology optimization as existing codes became outdated and computational bottlenecks shifted to assembly and filtering. It presents compact 2D and 3D codes with more efficient operations and discusses Anderson acceleration; the resulting implementations show remarkable speedups and allocate most runtime to the equilibrium solve.
Problem
Existing top88-based implementations may be outdated for increasingly challenging topology optimization problems, while matrix assembly and filter application remain computational bottlenecks.
Method
The paper develops compact Matlab codes using efficient assembly, filtering, projection, design-update shortcuts, and an optionally integrated Anderson acceleration strategy.
Results
85 to 90% of 2D runtime and 96% of 3D runtime is spent solving the equilibrium system, while the new codes show remarkable speedups over prior implementations.
Takeaways & Limitations
The implementations are computationally well balanced, leaving the dominant equilibrium solve to preconditioned iterative solvers and enabling acceleration of the overall optimization loop.
Takeaways & Limitations
The influence of acceleration parameters is not discussed in depth, and the authors note that small OC-update differences can produce local optima or bad designs in some settings.
Abstract
from arXiv · showhide
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\cdot 10^{4}$ to $4.8\cdot10^{5}$ elements, the 2D version, named top99neo, shows speedups from 2.55 to 5.5 times compared to the well-known top88 code (Andreassen-etal 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-etal 2014, on a discretization with $2.2\cdot10^{5}$ 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.
1 Introduction
The paper introduces compact Matlab implementations for efficient minimum-compliance topology optimization in 2D and 3D. They target medium- to large-scale problems by combining implementation shortcuts, speedups, and acceleration strategies.
- Motivation: The implementation targets medium- and large-scale topology optimization problems that can be solved efficiently on a laptop.The authors present the code as an exemplary collection of shortcuts and speedups rather than as new optimization theory.
- Motivation: Preconditioned iterative solvers have shifted the computational bottleneck toward matrix assembly and repeated filter application.
- New implementations: top99neo and top3d125 provide compact Matlab implementations for minimum-compliance topology optimization in 2D and 3D continua.The 2D code has 99 lines, while the 3D extension requires remarkably few changes.
- New implementations: The codes include passive domains, volume-preserving density projection, and continuation strategies for penalization and projection parameters.
- Acceleration: Anderson acceleration can be integrated with a few extra code lines and may produce major speedups in the overall optimization loop.
2 Problem formulation and solution scheme
The paper formulates compliance topology optimization on 2D/3D discretizations with active and passive element variables, filtered and projected physical densities, and a volume constraint. A nested iterative scheme solves equilibrium, evaluates sensitivities, and updates active variables using an Optimality Criterion rule with a bisection-determined Lagrange multiplier.
- Problem formulation: The discretization contains m equal-sized quadrilateral elements, with n global and d local degrees of freedom.
- Problem formulation: Design variables are partitioned into active variables and passive solid or void elements, while physical active variables are obtained through relaxed Heaviside projection of a filtered field.The filtered field is produced by a linear operator using neighboring elements within radius rmin and weights he,i.
- Problem formulation: The optimization minimizes compliance c(ˆx) = uT f subject to an upper bound on overall volume for a prescribed volume fraction.
- Solution scheme: At each nested-loop iteration, the displacement u is computed from the equilibrium problem using a stiffness matrix interpolated from physical variables through SIMP.
- Solution scheme: Active variables are updated by the Optimality Criterion rule using compliance and volume sensitivities, move limits, and an approximation to the volume-constraint Lagrange multiplier.The multiplier is obtained by approximately enforcing the next iterate’s volume constraint, for example by bisection.
3 Matlab implementation and speedups
The implementation reorganizes Matlab topology optimization around compact initialization, efficient filtering, sparse assembly, and accelerated design updates. These changes reduce memory and computation while preserving extensions such as projection, continuation, and passive domains.
- Implementation structure: top99neo organizes the routine into one-time initialization operations and an iterative topology-optimization redesign loop.Initialization includes material, discretization, loads, supports, passive domains, implicit functions, and other parameters.
- Implementation structure: The filter uses Matlab’s imfilter function instead of explicitly constructing a neighboring array.This is presented as a more efficient alternative to explicit filter construction.
- Assembly speedups: The H assembly strategy reduces indexing-array entries and memory by approximately 45%, while fsparse further improves sparse assembly efficiency.The H strategy exploits stiffness-matrix symmetry; fsparse improves operation sorting and supports int32 indexing.
- Assembly speedups: 10–15 times speedups are achieved in assembly by combining symmetry-based H assembly with fsparse, with linear CPU-time scaling in the number of DOFs.The reported speedups are 10 for 2D and 15 for 3D; H assembly alone cuts about half the CPU time.
- Acceleration: Anderson acceleration substantially reduces iterations to the stopping criterion while increasing compliance by only 0.2–3%.Starting acceleration later can yield lower compliance, and differences are negligible for T3 and T4 after design stabilization.
- Performance comparison: 2.66, 3.85 and 5.5 times speedups over top88 are obtained for the three tested discretizations.Relative to top88U, the corresponding improvements are 1.55, 1.57 and 1.78 times.
4 Extension to 3D
The 2D implementation extends to 3D with few code changes, while top3D125 substantially reduces assembly and update costs and achieves faster performance than the compared 3D code.
- Implementation: The 2D implementation extends to 3D problems with only 12 minor line modifications, mainly accounting for the extra spatial dimension.The 3D changes define the hexahedral elemental stiffness matrix and modify equilibrium-solver instructions.
- 3D example: For the cantilever test, zero-Neumann filtering reproduces the reference design, whereas zero-Dirichlet filtering produces a clearly worse local minimum.The designs use 100 redesign steps and include volume-preserving Heaviside projection.
- 3D example: 96% of the 1741 s total runtime is spent solving the state equation, while assembly takes 1.2% and filtering and design updates take 0.4%.This runtime profile reflects the 3D cantilever experiment over 100 iterations.
- Performance comparison: top3D125 achieves speedups of about 1.8 and 1.9 on the 48 × 24 × 24 and 96 × 48 × 48 discretizations, respectively.The comparison uses the same multigrid-preconditioned CG solver as Amir et al. (2014).
- Performance comparison: Matrix assembly accounts for 7–10% of top3D125 runtime, compared with about 50% in Amir et al. (2014).Most of the speedup comes from reducing assembly time; the OC update is also faster.
5 Concluding remarks
The paper presents compact 2D and 3D Matlab implementations that achieve substantial speedups through efficient assembly, filtering, and design updates, with iterative solvers addressing the dominant equilibrium-solve cost.
- Contributions: The new compliance Topology Optimization codes for 2D and 3D show remarkable speedups over top88 and available 3D implementations.The 3D comparison includes codes by Liu and Tovar (2014) and Amir et al. (2014).
- Contributions: Efficient assembly uses integer mesh quantities and assembles only one half of the stiffness matrix.The approach exploits matrix symmetry to reduce assembly work.
- Contributions: The explicit Lagrange-multiplier expression drastically cuts the number of OC iterations, while filtering and volume-preserving projection speed up redesign.These changes target the optimization update step rather than the equilibrium solve.
- Computational balance: As problem size increases, equilibrium-system solution consumes 85 to 90% of 2D runtime and 96% of 3D runtime.The computational balance makes preconditioned iterative solvers particularly relevant to the implementation.
- Scope: The methods can also apply, with stated extensions or qualifications, to level-set and other Topology Optimization approaches.The authors specifically note broader applicability of the assembly, update, filtering, and Anderson-acceleration ideas.
- Practical relevance: The codes are intended to help researchers and practitioners tackle medium/large-scale Topology Optimization problems on laptops.Matlab codes are listed in the Appendix and made available through the TopOpt website.
A Elaboration on the OC update
The OC update is derived through a primal-dual formulation in which separability yields elementwise updates subject to box constraints. For the volume constraint, the Lagrange multiplier can be computed directly in the simple case, while more involved constraints require root finding.
- Primal-dual formulation: The constrained approximation is formulated with a Lagrangian and solved by interlaced primal-dual iterations over densities and the Lagrange multiplier.The feasible set imposes lower and upper bounds on each design variable.
- Elementwise update: Separability reduces the density update to elementwise expressions that are then clipped according to lower-bound, upper-bound, and interior sets.The sets are denoted L, U, and M, with the interior update given by F^(j).
- Multiplier computation: The dual stationarity condition determines the Lagrange multiplier enforcing the volume constraint for the updated density.The derivation uses the characteristic function of the active sets and the number of elements in each set.
- General constraints: More involved constraint expressions do not generally yield a closed-form multiplier, so a root-finding strategy is required to approximate it.The paper states that bisection can be used to impose approximate volume feasibility.
- Connection to OC: For the simple volume constraint, the resulting expression is identical to the earlier OC update when applied at the current feasible design.The equivalence follows from the active-set cardinalities and equation (7).
B The 2D code for compliance minimization
The 2D top99neo implementation organizes compliance minimization into preprocessing, filtering and projection, equilibrium solving, sensitivity evaluation, OC design updates, and continuation. Its compact code uses reduced assembly indexing, sparse matrix construction, and convolution-based filters.
- Discretization and assembly: The code builds a 2D element connectivity structure and uses a symmetric elemental stiffness matrix with reduced assembly indexing.The elemental matrix is first represented through its lower triangular part and then completed.
- Filtering and projection: Physical densities are obtained by convolution filtering, optional projection, and projection-parameter continuation.The filter uses a radius-based kernel and selectable boundary conditions; Newton iterations determine eta when projection is active.
- Equilibrium solve: Each iteration interpolates element stiffness, assembles the global sparse stiffness matrix, and solves the free-degree-of-freedom equilibrium system by Cholesky decomposition.The interpolated stiffness is sK = Emin + xPhys^penal (E0 - Emin).
- Sensitivity and update: Compliance sensitivities and volume sensitivities are filtered before the OC resizing rule updates active design variables within move and box bounds.The Lagrange-multiplier estimate is refined by bisection until the volume condition is met, followed by continuation on penal and beta.
C The 3D code for compliance minimization
The 3D top3D125 implementation extends the same compliance-minimization workflow to hexahedral continua in a 125-line Matlab function. It uses 3D connectivity, sparse assembly, volumetric filtering, equilibrium solution, sensitivity filtering, and OC updates with continuation.
- Filtering and projection: Physical densities are produced by a radius-based 3D convolution filter, optional projection, Newton updates for eta, and sensitivity modification.The filter kernel depends on dx, dy, and dz, and the filtered-weight matrix Hs accounts for boundary conditions.
- Equilibrium solve: The equilibrium equations use interpolated elemental stiffnesses, sparse global assembly, Cholesky factorization, and forward/back substitution.The stiffness interpolation uses the physical density raised to the penalization power.
- Sensitivity, update, and output: The design update filters compliance and volume sensitivities, applies OC resizing under move and box bounds, and continues penalization and projection parameters.The 3D implementation also prints compliance, volume, change, continuation parameters, and the Lagrange multiplier, then plots an isosurface.