Topology optimization

Overview

Topology optimization answers a different question than sizing a part: not how thick should this beam be, but where should the material go at all. Given a design space, the supports, and the loads, the optimizer carves holes into the space until what remains is the lightest structure that still carries the load.

The method Rasmah uses is the density method (SIMP). It works on a fixed mesh — the mesh never changes during the optimization. Instead, every element carries a density $x_e \in [0,1]$ that means "how much material is here": $x_e = 1$ is solid, $x_e = 0$ is void. The element's stiffness is then scaled by that density, so a void element contributes almost nothing to the structure. The optimizer slowly pushes every $x_e$ toward $0$ or $1$ until the design is a crisp black-and-white layout.

Because the mesh is fixed, the whole $\text{CAD} \to \text{density} \to \text{FEM} \to \text{objective}$ chain is differentiable end to end — which is exactly what a gradient-based optimizer needs.

The density method (SIMP)

The idea is one line: make the material stiffness a function of a per-element density. In SIMP (Solid Isotropic Material with Penalization) the Young's modulus of element $e$ is

\[E(x_e) = x_e^{p} E_0 ,\]

with the penalty $p$ (usually $p = 3$). $p$ is what makes the result black-and-white: because the stiffness cost of a gray $x_e = 0.5$ element is $0.5^3 \approx 0.125$ of a solid element while its material cost is still $0.5$, gray material is inefficient, and the optimizer is pushed toward the extremes.

The optimization problem itself is to make the structure as stiff as possible under a material budget. Stiffness is measured by the compliance — the strain energy $C = \mathbf{f}^{\top}\mathbf{u}$, where $\mathbf{u}$ solves the equilibrium $\mathbf{K}(x)\,\mathbf{u} = \mathbf{f}$. A stiffer structure stores less energy under the same load, so we minimize compliance:

\[\min_{x} \; C(x) = \mathbf{f}^{\top}\mathbf{u}(x) \quad\text{subject to}\quad \frac{\sum_e x_e}{V_0} \le v_f , \quad 0 \le x_e \le 1 .\]

$V_0$ is the full design volume and $v_f$ the allowed volume fraction. The solver iterates: solve for $\mathbf{u}$, compute the density gradient of the compliance (one adjoint solve, $\partial C/\partial x_e = -p\,x_e^{p-1} \mathbf{u}_e^{\top}\mathbf{K}_e \mathbf{u}_e$), update $x$ with a bound-aware updater (MMA / BESO / TOBS), and repeat.

The design variable and the penalty

A design vector is a Densities, whose type parameters record where it is in the filter → project → penalize pipeline (so downstream code can chain-rule correctly):

Densities

The material interpolation is a penalty object. SIMP is PowerPenalty; RAMP, sinh, and the fluid Darcy interpolations are the others:

AbstractPenalty

Applying the penalty elementwise is penalize, with its derivative penalize_derivative:

penalize(d::Densities, p::AbstractPenalty) -> Densities

penalize_derivative(d::Densities, p::AbstractPenalty) -> Vector

d = Densities([0.0, 0.5, 1.0])
penalize(d, PowerPenalty(3.0))

3-element Densities{false, false, true, Float64, Vector{Float64}}:
 0.0
 0.125
 1.0

The mid-gray $x = 0.5$ becomes $0.125$ under $p = 3$, which is the whole mechanism that drives the design to black and white.

The fixed-mesh density method

The bridge between geometry and density is the signed distance field. Each element's density is read from an SDF at its centroid — a smooth sigmoid of the signed distance, or a hard inside/outside cutoff:

element_density(g::ImplicitGeometry, m::TetMesh; β = 0.1) -> vector

compliance_gradient_density(gfun, m::TetMesh, material, fixed_dofs, forces, θ; β = 0.1, p = 3) -> number

Because the mesh is fixed, the derivative of the compliance with respect to a design parameter of the SDF is computed analytically with a single adjoint solve — this is what keeps the CAD-to-objective chain differentiable.

Objectives and constraints

An objective or a constraint is an AbstractFunctional — a quantity that depends on the design through the state. Each exposes a value and a gradient, so one forward solve plus one adjoint yields every value and gradient at once:

AbstractFunctional

The two you use most are the compliance objective and the volume constraint:

The state map

The forward problem — "given densities, solve the physics" — is bundled into a state map. The linear (elasticity/heat) one is AffineStateMap:

AffineStateMap(form, load, fixed; penalty, method, reuse, solid, void, void_density)

The state map owns the assembly and the solve, and reuses its factorization (or multigrid hierarchy) across iterations, since SIMP scaling changes values but not the sparsity pattern.

Assembling a problem and running it

A complete problem couples a state map to an objective and constraints:

TopOptProblem(state, objective, constraints...; xmin=1e-3, xmax=1.0)

The problem is run by optimize(prob; algorithm=…), which picks the updater — :mma (the method of moving asymptotes, the default), :beso (bi-directional evolutionary), or :tobs (topology optimization of binary structures) — and returns (; x, history) with the per-element densities and the objective trace.

The high-level driver

In practice you rarely assemble the state map by hand. Name the physics and pass a material, exactly as you would for a forward solve:

ComplianceTopOpt <: AbstractPhysics

ThermalComplianceTopOpt <: AbstractPhysics

m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4)
n = size(m.nodes, 2)
fixed = Int[]
for i in 1:n
    m.nodes[1, i] < 1e-6 && append!(fixed, 3i - 2, 3i - 1, 3i)
end
loads = zeros(3n)
for i in 1:n
    m.nodes[1, i] > 0.999 && (loads[3i - 2] = 1e6)
end
x, history = solve(ComplianceTopOpt(), m, steel; fixed=fixed, loads=loads, volfrac=0.5, p=3.0, max_iter=40)

([0.9999998030733432, 0.9999995273760242, 0.001014250805391356, 0.001014250805391356, 0.9999995273760242, 0.9999998030733435, 0.9999997444254357, 0.9999920225237264, 0.0010714750410144973, 0.0010714750410144973  …  0.9999995273760245, 0.9999995273760245, 0.9999967860574556, 0.9977344920225029, 0.999880839989524, 0.9999952584935178, 0.9999983397707033, 0.9999983397707034, 0.9999952584935178, 0.9998808399895237], [25866.825770885473 J, 17972.32404280367 J, 28724.23298207425 J, 24814.1894542239 J, 23940.241963934193 J, 23637.48329407373 J, 23344.049662433878 J, 22989.325490353247 J, 22546.133193639715 J, 21975.539359017883 J, …, 6497.457785910194 J, 6487.546035703661 J, 6474.965555175202 J, 6474.569316606495 J, 6525.154466195304 J, 6498.546244502275 J, 6502.798467943279 J, 6488.725932649259 J, 6499.0584114162975 J, 6491.860990051545 J])

history[1], history[end]

(25866.825770885473 J, 6491.860990051545 J)

The compliance falls by about $4\times$ as the optimizer carves a stiff structure out of half the box. x is the per-element density ($\approx 0.5$ on average, saturating the volume budget), ready to be visualized with vtk(m; field=x).

Filters and projection

Raw density optimization produces checkerboards — alternating solid/void elements that are a numerical artifact, not a real structure. Two tools clean them up: a density (or sensitivity) filter that averages each element's neighborhood, and a smooth Heaviside projection that sharpens the result back to 0/1:

filter_sensitivity(f::SensitivityFilter, dc::AbstractVector, x::AbstractVector) -> Vector

filter_density(f::DensityFilter, x::AbstractVector) -> Vector

project_derivative(f::HeavisideProjection, x::AbstractVector) -> Vector

The wider family

Everything above is the single-material, linear-elasticity baseline, but the same solve(Physics(), mesh, material; …) entry point exposes dozens of specialized methods: thermal and eigenfrequency design, stress-constrained and buckling design, multi-material (DMO) and moving-morphable-component (MMC/MMV) methods, fluid (Stokes/Navier–Stokes/Darcy) topology, level-set and phase-field methods, robust and overhang-constrained design for additive manufacturing, and functional- device inverse design for acoustics, electromagnetics, and photonics (AcousticCloakTopOpt, PhotonicBandGapTopOpt, ElectromagneticInverseDesignTopOpt, and many more). Each follows the same pattern — a physics marker, a material, and the boundary conditions — with a density (or level-set) field as the design variable.

Next steps

Topology optimization sits on top of the forward FEM & simulation chapter and consumes the Materials library it interpolates. The updaters that drive the densities — MMA, BESO, TOBS — are general optimizers, covered in the Optimization & nonlinear solvers chapter.