Density-Based Topology Optimization

Overview

The density method keeps the mesh fixed and gives each element a material fraction φᵉ ∈ [0, 1]; solid elements (φᵉ ≈ 1) carry full stiffness, void elements (φᵉ ≈ 0) carry none. This tutorial drives that with Rasmah’s native SIMP (Solid Isotropic Material with Penalization) topology-optimization engine, built on the same weak-form FEM and the same region => value boundary-condition syntax as the forward solves. Three optimizers are available: MMA (default), BESO, and TOBS.

The worked part is a design job an engineer would run: a steel engine-mount bracket, a 120 × 60 × 8 mm plate bolted to the block on its left face and pulling 4800 N (a 10 MPa traction) on its right face. The full plate weighs 452 g; the optimizer must find the lightest layout that keeps 50% of the material while minimizing the compliance — the measure of how “soft” the mount has become.

using Rasmah
using LinearAlgebra

g(v) = round(ustrip(uconvert(u"kg", v)) * 1000; digits = 0)

g (generic function with 1 method)

Geometry, face selection, and meshing

The workflow mirrors the forward FEM: build the geometry, take its exact BRep, select the faces before meshing, then mesh. The bracket is a plain box in millimetres; the two faces_in_box calls grab the bolt faces — the left face inside the clamped end, the right face inside the loaded end:

# 1. the geometry: 120 × 60 × 8 mm bracket plate
part = box(0.120, 0.060, 0.008)

# 2. exact boundary representation + face selection
b = brep(part)
left  = faces_in_box(b, [-0.062, -1.0, -1.0], [-0.058, 1.0, 1.0])
right = faces_in_box(b, [ 0.058, -1.0, -1.0], [ 0.062, 1.0, 1.0])

# 3. mesh the geometry — the face sets carry over
mesh = tetrahedralize(b, 20)
println(size(mesh.elements, 2), " tetrahedra")

1200 tetrahedra

The BRep renders directly:

vtk(b)

Zoom Opacity

The bill of materials is the closed-form volume — the starting weight the optimizer has to spend:

V_full = 0.120 * 0.060 * 0.008
println("full-plate mass (steel) = ", g(V_full * 7850u"kg/m^3"), " g")

full-plate mass (steel) = 452.0 g

Elasticity compliance minimization (SIMP + MMA)

solve(ComplianceTopOpt(), …) minimizes the mechanical compliance fᵀu subject to a volume-fraction constraint, using the same fixed/loads pair syntax as solve(Elasticity(), …). The 10 MPa traction on the 60 × 8 mm face is a 4800 N pull — a realistic service load for a bracket of this size. It returns the per-element density field and the compliance history:

φ, hist = solve(ComplianceTopOpt(), mesh, steel;
                fixed = left => [0.0, 0.0, 0.0],
                loads = right => [1e7, 0.0, 0.0],
                volfrac = 0.5, max_iter = 40)

println("initial compliance    = ", round(hist[1]; digits = 3))
println("final compliance      = ", round(hist[end]; digits = 3))
println("final volume fraction = ", round(sum(φ) / length(φ); digits = 4))

initial compliance    = 0.228 J
final compliance      = 0.065 J
final volume fraction = 0.4999

The density concentrates into load-bearing members while the compliance drops — the layout got stiffer at half the weight. The mass budget closes exactly: 50% of 452 g:

println("mass saved ≈ ", g(0.5 * V_full * 7850u"kg/m^3"), " g  (full plate ", g(V_full * 7850u"kg/m^3"), " g)")

mass saved ≈ 226.0 g  (full plate 452.0 g)

Plot the final density field on the mesh:

vtk(mesh; field = φ, fieldname = "material density")

Zoom Opacity

Alternative optimizers

algorithm = :beso does bi-directional evolutionary add/remove (soft-kill {φmin, 1} densities), and algorithm = :tobs produces binary designs by sequential linearization with integer cuts:

φ_beso, _ = solve(ComplianceTopOpt(), mesh, steel;
                  fixed = left => [0.0, 0.0, 0.0], loads = right => [1e7, 0.0, 0.0],
                  volfrac = 0.5, max_iter = 40, algorithm = :beso)
vtk(mesh; field = φ_beso, fieldname = "BESO density")

Zoom Opacity

φ_tobs, _ = solve(ComplianceTopOpt(), mesh, steel;
                  fixed = left => [0.0, 0.0, 0.0], loads = right => [1e7, 0.0, 0.0],
                  volfrac = 0.5, max_iter = 40, algorithm = :tobs)
vtk(mesh; field = φ_tobs, fieldname = "TOBS density")

Zoom Opacity

Thermal compliance minimization

solve(ThermalComplianceTopOpt(), …) is the scalar (heat-conduction) twin: it minimizes the thermal compliance qᵀT over the same density field. The bracket now acts as a copper heat spreader: a chip-scale 1 MW/m² flux enters the left face, the right is held at T = 0:

# heat: a 1 MW/m² flux is injected on the left, the right is held at T = 0
φT, histT = solve(ThermalComplianceTopOpt(), mesh, copper;
                 fixed = right => 0.0, loads = left => 1e6,
                 volfrac = 0.5, max_iter = 40)

println("final thermal compliance = ", round(histT[end]; digits = 2))

final thermal compliance = 1.14912716e6 m^2·kg·K/s^3

vtk(mesh; field = φT, fieldname = "thermal density")

Zoom Opacity

The underlying framework

ComplianceTopOpt/ThermalComplianceTopOpt are thin wrappers over a general epigraph problem: a TopOptProblem bundles a state map (AffineStateMap) with an objective and constraints, and optimize runs the chosen updater. Everything is expressed as BilinearForm/FEOperator weak forms, so the same machinery serves any single-field linear physics.

Penalties, filters, and projections plug into the same problem, and the typed Densities{F,P,I} pipeline exposes compliance at element/block/global granularity:

# penalize (SIMP), filter, and project are independent, composable stages
ρ = Densities(φ)                       # raw per-element densities
penalized = penalize(ρ, PowerPenalty(3.0))
println("penalized[1..3] = ", penalized[1:3])

filter = SensitivityFilter(mesh, 1.5)  # 99-line `check` filter radius
println("filter ready: ", filter isa SensitivityFilter)

proj = HeavisideProjection(8.0)        # smooth 0/1 projection
println("projection ready: ", proj isa HeavisideProjection)

penalized[1..3] = [1.0, 1.0, 1.0]
filter ready: true
projection ready: true

The optimizer is selected by algorithm = :mma | :beso | :tobs; the objective and constraints are evaluated together in one state solve plus one adjoint bundle (the epigraph form min c s.t. f(x) ≤ c).

Next steps

Heat Conduction and SIMP Topology Optimization — the thermal solve and its adjoint gradient in detail.
Shape optimization — the boundary-moving counterpart.
Linear elasticity — the forward FEM behind ComplianceTopOpt.

References

Bendsøe, M. P., & Sigmund, O. (2003). Topology Optimization: Theory, Methods, and Applications (2nd ed.). Springer, Berlin. https://doi.org/10.1007/978-3-662-05086-6
Svanberg, K. (1987). The method of moving asymptotes. Int. J. Numer. Methods Eng. 24(2): 359–373.
Sivapuram, R., & Picelli, R. (2018). Topology optimization of binary structures using integer linear programming. Finite Elem. Anal. Des. 139: 49–61.