using TopOpt, Zygote, TestMulti-Material Topology Optimization with Softmax Parametrization
Description
This tutorial demonstrates multi-material topology optimization, where the design distributes several candidate materials (plus void) over the domain. Each element holds a softmax over nmats - 1 decision variables, ensuring material fractions sum to one. Material interpolation maps these fractions to Young’s moduli (for compliance) and physical densities (for mass constraints).
We optimize a cantilever beam with 3 materials: void (E=1e-5), soft material (E=1.0), and stiff material (E=4.0), subject to a 40% mass constraint.
Setup
Define the Problem
Es = [1e-5, 1.0, 4.0] # Young's moduli: void, soft, stiff
densities = [0.0, 0.5, 1.0] # physical densities for mass constraint
nmats = 3 # number of materials
nu = 0.3 # Poisson's ratio
f = 1.0 # downward force
problem = PointLoadCantilever(
Val{:Linear}, # linear elements
(160, 40), # 160×40 mesh
(1.0, 1.0), # element size
1.0, # base Young's modulus
nu, # Poisson's ratio
f, # load
)
ncells = getncells(problem)FEA Solver and Filter
solver = FEASolver(DirectSolver, problem; xmin=0.0)
filter = DensityFilterFun(solver; rmin=4.0)
comp = ComplianceFun(solver)Material Interpolation
# Interpolation 1: Young's moduli for compliance (penalty=3.0)
penalty1 = PowerPenaltyFun(3.0)
interp1 = MaterialInterpolationFun(Es, penalty1)
# Interpolation 2: Physical densities for mass constraint (penalty=1.0)
penalty2 = PowerPenaltyFun(1.0)
interp2 = MaterialInterpolationFun(densities, penalty2)Two interpolations serve different purposes:
interp1maps softmax outputs to stiffness for FEAinterp2maps to physical densities for mass calculation
Objective and Constraint
obj = y -> begin
x = tounit(MultiMaterialVariablesFun(y, nmats))
_E = interp1(filter(x))
return comp(_E)
end
y0 = zeros(ncells * (nmats - 1))
# Sanity-check objective and gradient
println("Initial objective: $(obj(y0))")
grad_obj = Zygote.gradient(obj, y0)[1]
println("Gradient norm: $(norm(grad_obj))")
constr = y -> begin
_rhos = interp2(MultiMaterialVariablesFun(y, nmats))
return sum(_rhos.x) / ncells - 0.4 # 40% mass constraint
end
println("Initial constraint: $(constr(y0))")
grad_constr = Zygote.gradient(constr, y0)[1]
println("Constraint gradient norm: $(norm(grad_constr))")The softmax parametrization ensures material fractions sum to one at each element. Design variables y are unconstrained (−10 to 10), transformed to material fractions via softmax.
Optimization
model = Model(obj)
addvar!(model, fill(-10.0, length(y0)), fill(10.0, length(y0)))
add_ineq_constraint!(model, constr)
alg = MMA87()
options = MMAOptions(; s_init=0.1, tol=Tolerance(; kkt=1e-3))
y0 = zeros(ncells * (nmats - 1))
res = optimize(model, alg, y0; options=options)
y = res.minimizer[ Info: iter obj Δobj violation kkt_residual [ Info: 0 1.3e+03 Inf 1.0e-01 1.7e+00 [ Info: 1 1.5e+03 1.5e+02 0.0e+00 9.5e+02 [ Info: 2 5.6e+02 9.1e+02 0.0e+00 2.3e+02 [ Info: 3 3.4e+02 2.2e+02 0.0e+00 5.1e+01 [ Info: 4 2.6e+02 7.7e+01 0.0e+00 1.7e+01 [ Info: 5 2.3e+02 3.3e+01 0.0e+00 7.9e+00 [ Info: 6 2.1e+02 1.7e+01 0.0e+00 4.4e+00 [ Info: 7 2.0e+02 1.0e+01 0.0e+00 3.0e+00 [ Info: 8 1.9e+02 7.4e+00 0.0e+00 2.2e+00 [ Info: 9 1.9e+02 5.6e+00 0.0e+00 1.7e+00 [ Info: 10 1.8e+02 4.5e+00 0.0e+00 1.3e+00 [ Info: 11 1.8e+02 3.6e+00 0.0e+00 1.1e+00 [ Info: 12 1.8e+02 3.0e+00 0.0e+00 8.7e-01 [ Info: 13 1.8e+02 2.4e+00 0.0e+00 7.3e-01 [ Info: 14 1.7e+02 2.2e+00 0.0e+00 4.6e-01 [ Info: 15 1.7e+02 1.9e+00 0.0e+00 2.6e-01 [ Info: 16 1.7e+02 1.8e+00 4.3e-05 3.3e-02 [ Info: 17 1.7e+02 1.3e+00 5.2e-05 3.8e-02 [ Info: 18 1.7e+02 5.6e-01 0.0e+00 1.8e+00 [ Info: 19 1.7e+02 2.2e+00 0.0e+00 8.9e-01 [ Info: 20 1.7e+02 8.2e-01 0.0e+00 7.8e-01 [ Info: 21 1.6e+02 1.2e+00 0.0e+00 2.4e-01 [ Info: 22 1.6e+02 4.7e-01 0.0e+00 2.3e-01 [ Info: 23 1.6e+02 5.8e-01 0.0e+00 7.0e-02 [ Info: 24 1.6e+02 4.2e-01 0.0e+00 4.8e-02 [ Info: 25 1.6e+02 3.4e-01 0.0e+00 1.1e-01 [ Info: 26 1.6e+02 4.2e-01 0.0e+00 1.0e-01 [ Info: 27 1.6e+02 4.7e-01 0.0e+00 9.5e-02 [ Info: 28 1.6e+02 5.4e-01 0.0e+00 9.3e-02 [ Info: 29 1.6e+02 5.3e-01 0.0e+00 1.4e-01 [ Info: 30 1.6e+02 6.5e-01 0.0e+00 9.3e-02 [ Info: 31 1.6e+02 6.4e-01 0.0e+00 1.1e-01 [ Info: 32 1.6e+02 8.6e-01 0.0e+00 3.4e-02 [ Info: 33 1.6e+02 9.1e-01 5.3e-05 3.2e-02 [ Info: 34 1.6e+02 8.6e-01 4.6e-05 2.7e-02 [ Info: 35 1.6e+02 9.7e-01 0.0e+00 6.9e-02 [ Info: 36 1.5e+02 2.0e+00 8.4e-04 4.5e-01 [ Info: 37 1.5e+02 1.4e+00 1.4e-03 7.1e-01 [ Info: 38 1.5e+02 6.8e-03 0.0e+00 4.6e-02 [ Info: 39 1.5e+02 2.1e-01 0.0e+00 2.5e-01 [ Info: 40 1.5e+02 3.6e-01 0.0e+00 1.6e-01 [ Info: 41 1.5e+02 2.2e-01 0.0e+00 1.3e-01 [ Info: 42 1.5e+02 1.3e-01 0.0e+00 1.2e-01 [ Info: 43 1.5e+02 3.3e-01 2.2e-04 1.1e-01 [ Info: 44 1.5e+02 7.7e-02 0.0e+00 2.1e-02 [ Info: 45 1.5e+02 6.8e-02 1.9e-05 9.9e-03 [ Info: 46 1.5e+02 4.3e-02 4.3e-05 2.2e-02 [ Info: 47 1.5e+02 2.6e-02 0.0e+00 2.8e-02 [ Info: 48 1.5e+02 3.7e-02 0.0e+00 4.7e-03 [ Info: 49 1.5e+02 5.5e-02 6.3e-05 3.3e-02 [ Info: 50 1.5e+02 7.1e-02 0.0e+00 6.1e-02 [ Info: 51 1.5e+02 6.3e-02 0.0e+00 2.0e-02 [ Info: 52 1.5e+02 3.4e-02 1.6e-05 8.5e-03 [ Info: 53 1.5e+02 2.9e-03 1.5e-05 7.7e-03 [ Info: 54 1.5e+02 7.1e-03 0.0e+00 3.2e-03 [ Info: 55 1.5e+02 1.5e-02 1.4e-05 7.5e-03 [ Info: 56 1.5e+02 1.5e-02 3.5e-05 1.8e-02 [ Info: 57 1.5e+02 2.6e-03 2.2e-05 1.2e-02 [ Info: 58 1.5e+02 1.1e-02 0.0e+00 3.6e-03 [ Info: 59 1.5e+02 2.3e-02 0.0e+00 4.0e-02 [ Info: 60 1.5e+02 5.0e-02 0.0e+00 9.7e-02 [ Info: 61 1.5e+02 7.8e-02 0.0e+00 2.4e-02 [ Info: 62 1.5e+02 5.2e-03 0.0e+00 2.5e-02 [ Info: 63 1.5e+02 2.4e-02 0.0e+00 6.9e-03 [ Info: 64 1.5e+02 3.5e-03 0.0e+00 7.3e-03 [ Info: 65 1.5e+02 1.6e-03 0.0e+00 9.7e-03 [ Info: 66 1.5e+02 2.9e-02 0.0e+00 5.0e-02 [ Info: 67 1.5e+02 2.6e-02 0.0e+00 3.7e-02 [ Info: 68 1.5e+02 8.1e-02 7.1e-05 3.8e-02 [ Info: 69 1.5e+02 8.2e-03 7.4e-05 3.9e-02 [ Info: 70 1.5e+02 1.5e-02 4.0e-05 2.1e-02 [ Info: 71 1.5e+02 1.4e-02 1.0e-05 5.3e-03 [ Info: 72 1.5e+02 2.5e-03 3.5e-06 1.9e-03 [ Info: 73 1.5e+02 1.8e-03 5.4e-06 2.8e-03 [ Info: 74 1.5e+02 3.1e-03 9.2e-06 4.9e-03 [ Info: 75 1.5e+02 1.6e-02 3.4e-05 1.8e-02 [ Info: 76 1.5e+02 1.7e-02 0.0e+00 3.1e-04
12800-element Vector{Float64}:
-10.0
-10.0
-10.0
-10.0
-10.0
-10.0
-10.0
-10.0
-10.0
-10.0
⋮
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0Verify Results
println("Constraint value: $(constr(y))")
x = tounit(reshape(y, ncells, nmats - 1))
println("Non-void fraction: $(sum(x[:, 2:3]) / size(x, 1))")
# Verify constraints are satisfied
println("Mass constraint satisfied: $(constr(y) < 1e-6)")
println("Softmax sums to 1: $(all(x -> isapprox(x, 1), sum(x; dims=2)))")Constraint value: -1.0341408696046273e-7
Non-void fraction: 0.5683926946580995
Mass constraint satisfied: true
Softmax sums to 1: trueVisualization
The visualize function accepts a per-cell topology vector (length ncells) that controls each cell’s transparency, and an optional cell_colors vector (same length) mapped through a colormap. For multi-material designs we first convert the raw decision variables y into per-cell material fractions via tounit, yielding an ncells × nmats matrix whose rows sum to 1. Each column is one material’s fraction per element, which we can pass to visualize to inspect that material in isolation.
using CairoMakie
x = tounit(MultiMaterialVariablesFun(y, nmats)) # ncells × nmats, rows sum to 1
rhos = x.x * densities # physical density per cell, used for transparency
fig = visualize(problem; topology=rhos, cell_colors=interp1(x).x, draw_legend=true)Precompiling packages... 6245.0 ms ✓ QuartoNotebookWorkerMakieExt (serial) 1 dependency successfully precompiled in 6 seconds Precompiling packages... 5063.6 ms ✓ QuartoNotebookWorkerCairoMakieExt (serial) 1 dependency successfully precompiled in 5 seconds
Here topology=rhos sets each cell’s alpha to its physical density (0 for void, 0.5 for soft, 1.0 for stiff), and cell_colors is the per-cell Young’s modulus so the colormap distinguishes the three materials. To inspect a single material, pass that material’s fraction column as the topology and mask the others:
stiff_frac = x.x[:, 3] # fraction of stiff material per cell
fig_stiff = visualize(problem; topology=stiff_frac)
soft_frac = x.x[:, 2] # fraction of soft material per cell
fig_soft = visualize(problem; topology=soft_frac)
The combined view shows optimal material distribution — stiff material concentrates in high-stress regions, soft material in intermediate areas, and void where material is inefficient. The per-material figures isolate each phase, making it easier to verify the design.