using TopOpt, LinearAlgebraLocal Stress: Element-Wise Constraints with Percival Solver
Description
This tutorial demonstrates local stress constraints in topology optimization. Unlike global stress aggregation (p-norm or KS; see the global_stress tutorial), local constraints enforce stress limits at every element individually. This provides stricter stress control but requires handling thousands of constraints efficiently.
Two ingredients from the stress-constrained literature are essential even in the local setting (Duysinx & Bendsøe, 1998; Le et al., 2010):
- Relaxed stress. The raw microscopic stress
σ = C_0 : εstays finite as an element’s density vanishes, so the local constraintσ_e ≤ σ_limhas singular optima: gradient methods cannot remove a stressed member. We therefore constrain the relaxed stressσ̃_e = ρ_e^q σ_e(stress_exponent = q = 0.5), which vanishes at zero density and makes the true optima reachable. - A calibrated threshold. The stress limit must be on the scale of the actual stresses in the structure — an arbitrary value is either vacuous (never active) or makes the problem infeasible. We calibrate it from the solid design below.
We use Percival.jl, a specialized solver for large-scale constrained optimization, with continuation SIMP (penalty 1.0 → 3.0) to converge to stress-constrained designs.
Setup
Define the Problem
E = 1.0 # Young's modulus
v = 0.3 # Poisson's ratio
f = 1.0 # downward force
rmin = 3.0 # filter radius
problem = PointLoadCantilever(Val{:Linear}, (80, 20), (1.0, 1.0), E, v, f)A 80×20 mesh (1,600 elements) keeps 1,600 local stress constraints affordable for repeated subproblem solves in CI.
Parameter Settings
V = 0.5 # volume fraction
xmin = 0.0001 # minimum density
x0 = fill(0.5, 80 * 20) # initial design
N = length(x0)
penalty = PowerPenaltyFun(1.0) # start with no penalization
solver = FEASolver(DirectSolver, problem; xmin=xmin, penalty=penalty)
# Relaxed stress (q = 0.5): removes the stress singularity at low densities
stress = von_mises_stress_function(solver; stress_exponent=0.5)
filter = DensityFilterFun(solver; rmin=rmin)
volfrac = VolumeFun(solver)Local Stress Constraints
comp = ComplianceFun(solver)
λ = 1e-4
# Minimize volume with a compliance regularization term. The regularizer
# steers the optimizer toward stiff layouts, preventing the trivial
# all-void solution that pure volume minimization would produce.
obj = x -> volfrac(filter(PseudoDensities(x))) - V + λ * comp(filter(PseudoDensities(x)))
# Calibrate the stress limit from the solid design's peak relaxed stress:
# slightly above it, so the constraint is active but the problem feasible
thr = 1.2 * maximum(stress(filter(PseudoDensities(ones(N)))))
println("stress threshold: $thr")
# Local constraints: σ̃ᵢ/thr - 1 ≤ 0 for each element (scaled to O(1)
# magnitudes; scaling by N instead would dilute violations below the
# solver's tolerances)
constr = x -> begin
s = stress(filter(PseudoDensities(x)))
return s ./ thr .- 1
endstress threshold: 1.223183190956383#5 (generic function with 1 method)Each element’s relaxed von Mises stress must satisfy the threshold. The constraint function returns a vector (one constraint per element), which Percival handles efficiently using interior-point methods. The compliance regularizer (λ·uᵀKu) penalizes overly compliant (near-void) designs, keeping the optimization away from the degenerate empty-domain solution.
Optimization with Continuation
alg = PercivalAlg()
options = PercivalOptions(; max_iter=20, subsolver_max_eval=100)
model = Model(obj)
addvar!(model, zeros(N), ones(N))
add_ineq_constraint!(model, constr)
x = copy(x0)
for p in [1.0, 2.0, 3.0]
setpenalty!(solver, p)
global r = optimize(model, alg, x; options=options)
global x = r.minimizer
println("Penalty $p: obj=$(obj(x)), max_stress=$(maximum(stress(filter(PseudoDensities(x)))))")
end[ Info: iter fx normgp normcx μ normy sumc inner_status iter_type [ Info: 0 4.9e-02 4.2e+01 2.8e+01 1.0e+01 4.0e+01 5 [ Info: 1 4.9e-02 4.2e+01 2.8e+01 1.0e+02 4.0e+01 9 first_order update_μ [ Info: 2 -2.2e-02 2.8e+01 4.5e-01 1.0e+02 3.0e+01 36 first_order update_y [ Info: 3 -1.2e-01 2.8e+01 3.0e-01 1.0e+03 3.0e+01 285 max_eval update_μ [ Info: 4 -1.2e-01 1.4e+02 3.0e-01 1.0e+03 2.7e+02 289 max_eval update_y [ Info: 5 -1.2e-01 1.4e+02 3.0e-01 1.0e+04 2.7e+02 293 max_eval update_μ [ Info: 6 -1.2e-01 1.7e+03 3.0e-01 1.0e+04 3.3e+03 297 max_eval update_y [ Info: 7 -1.2e-01 1.7e+03 3.0e-01 1.0e+05 3.3e+03 301 max_eval update_μ [ Info: 8 -1.2e-01 1.7e+04 3.0e-01 1.0e+05 3.3e+04 305 max_eval update_y [ Info: 9 -1.2e-01 1.7e+04 3.0e-01 1.0e+06 3.3e+04 309 max_eval update_μ [ Info: 10 -1.2e-01 1.7e+04 3.0e-01 1.0e+07 3.3e+04 313 max_eval update_μ [ Info: 11 -1.2e-01 1.7e+04 3.0e-01 1.0e+08 3.3e+04 317 max_eval update_μ [ Info: 12 -1.2e-01 1.7e+04 3.0e-01 1.0e+09 3.3e+04 321 max_eval update_μ [ Info: 13 -1.2e-01 1.7e+04 3.0e-01 1.0e+10 3.3e+04 325 max_eval update_μ [ Info: 14 -1.2e-01 1.7e+04 3.0e-01 1.0e+11 3.3e+04 329 max_eval update_μ [ Info: 15 -1.2e-01 1.7e+04 3.0e-01 1.0e+12 3.3e+04 333 max_eval update_μ [ Info: 16 -1.2e-01 1.7e+04 3.0e-01 1.0e+13 3.3e+04 337 max_eval update_μ [ Info: 17 -1.2e-01 1.7e+04 3.0e-01 1.0e+14 3.3e+04 341 max_eval update_μ [ Info: 18 -1.2e-01 1.7e+04 3.0e-01 1.0e+15 3.3e+04 345 max_eval update_μ [ Info: 19 -1.2e-01 1.7e+04 3.0e-01 1.0e+16 3.3e+04 349 max_eval update_μ [ Info: 20 -1.2e-01 1.7e+04 3.0e-01 1.0e+17 3.3e+04 353 max_eval update_μ Penalty 1.0: obj=-0.12354066763399325, max_stress=1.2329834248008038 [ Info: iter fx normgp normcx μ normy sumc inner_status iter_type [ Info: 0 4.6e-02 4.2e+01 4.0e+01 1.0e+01 4.0e+01 5 [ Info: 1 4.6e-02 4.2e+01 4.0e+01 1.0e+02 4.0e+01 9 first_order update_μ [ Info: 2 7.7e-02 6.4e+00 4.4e-01 1.0e+02 1.6e+01 132 first_order update_y [ Info: 3 8.3e-02 6.3e+00 1.6e-01 1.0e+03 1.6e+01 295 max_eval update_μ [ Info: 4 8.3e-02 1.4e+02 1.6e-01 1.0e+03 1.4e+02 299 max_eval update_y [ Info: 5 8.3e-02 1.4e+02 1.6e-01 1.0e+04 1.4e+02 303 max_eval update_μ [ Info: 6 8.3e-02 1.6e+03 1.6e-01 1.0e+04 1.7e+03 307 max_eval update_y [ Info: 7 8.3e-02 1.6e+03 1.6e-01 1.0e+05 1.7e+03 311 max_eval update_μ [ Info: 8 8.3e-02 1.7e+04 1.6e-01 1.0e+05 1.7e+04 315 max_eval update_y [ Info: 9 8.3e-02 1.7e+04 1.6e-01 1.0e+06 1.7e+04 319 max_eval update_μ [ Info: 10 8.3e-02 1.7e+05 1.6e-01 1.0e+06 1.7e+05 323 max_eval update_y [ Info: 11 8.3e-02 1.7e+05 1.6e-01 1.0e+07 1.7e+05 327 max_eval update_μ [ Info: 12 8.3e-02 1.7e+06 1.6e-01 1.0e+07 1.7e+06 331 max_eval update_y [ Info: 13 8.3e-02 1.7e+06 1.6e-01 1.0e+08 1.7e+06 335 max_eval update_μ [ Info: 14 8.3e-02 1.7e+07 1.6e-01 1.0e+08 1.7e+07 339 max_eval update_y [ Info: 15 8.3e-02 1.7e+07 1.6e-01 1.0e+09 1.7e+07 343 max_eval update_μ [ Info: 16 8.3e-02 1.7e+07 1.6e-01 1.0e+10 1.7e+07 347 max_eval update_μ [ Info: 17 8.3e-02 1.7e+07 1.6e-01 1.0e+11 1.7e+07 351 max_eval update_μ [ Info: 18 8.3e-02 1.7e+07 1.6e-01 1.0e+12 1.7e+07 355 max_eval update_μ [ Info: 19 8.3e-02 1.7e+07 1.6e-01 1.0e+13 1.7e+07 359 max_eval update_μ [ Info: 20 8.3e-02 1.7e+07 1.6e-01 1.0e+14 1.7e+07 363 max_eval update_μ [ Info: 21 8.3e-02 1.7e+07 1.6e-01 1.0e+15 1.7e+07 367 max_eval update_μ Penalty 2.0: obj=0.08343657183799422, max_stress=1.2197165158650343 [ Info: iter fx normgp normcx μ normy sumc inner_status iter_type [ Info: 0 1.4e-01 4.6e+01 5.3e+01 1.0e+01 4.0e+01 5 [ Info: 1 1.4e-01 4.6e+01 5.3e+01 1.0e+02 4.0e+01 9 first_order update_μ [ Info: 2 3.4e-01 1.1e+01 4.6e-01 1.0e+02 2.2e+01 80 first_order update_y [ Info: 3 3.5e-01 1.0e+01 2.2e-01 1.0e+03 2.2e+01 295 max_eval update_μ [ Info: 4 3.5e-01 1.8e+02 2.2e-01 1.0e+03 1.9e+02 299 max_eval update_y [ Info: 5 3.5e-01 1.8e+02 2.2e-01 1.0e+04 1.9e+02 303 max_eval update_μ [ Info: 6 3.5e-01 2.2e+03 2.2e-01 1.0e+04 2.4e+03 307 max_eval update_y [ Info: 7 3.5e-01 2.2e+03 2.2e-01 1.0e+05 2.4e+03 311 max_eval update_μ [ Info: 8 3.5e-01 2.2e+04 2.2e-01 1.0e+05 2.4e+04 315 max_eval update_y [ Info: 9 3.5e-01 2.2e+04 2.2e-01 1.0e+06 2.4e+04 319 max_eval update_μ [ Info: 10 3.5e-01 2.2e+05 2.2e-01 1.0e+06 2.4e+05 323 max_eval update_y [ Info: 11 3.5e-01 2.2e+05 2.2e-01 1.0e+07 2.4e+05 327 max_eval update_μ [ Info: 12 3.5e-01 2.2e+05 2.2e-01 1.0e+08 2.4e+05 331 max_eval update_μ [ Info: 13 3.5e-01 2.2e+05 2.2e-01 1.0e+09 2.4e+05 335 max_eval update_μ [ Info: 14 3.5e-01 2.2e+05 2.2e-01 1.0e+10 2.4e+05 339 max_eval update_μ [ Info: 15 3.5e-01 2.2e+05 2.2e-01 1.0e+11 2.4e+05 343 max_eval update_μ [ Info: 16 3.5e-01 2.2e+05 2.2e-01 1.0e+12 2.4e+05 347 max_eval update_μ [ Info: 17 3.5e-01 2.2e+05 2.2e-01 1.0e+13 2.4e+05 351 max_eval update_μ [ Info: 18 3.5e-01 2.2e+05 2.2e-01 1.0e+14 2.4e+05 355 max_eval update_μ [ Info: 19 3.5e-01 2.2e+05 2.2e-01 1.0e+15 2.4e+05 359 max_eval update_μ [ Info: 20 3.5e-01 2.2e+05 2.2e-01 1.0e+16 2.4e+05 363 max_eval update_μ [ Info: 21 3.4e-01 2.2e+05 2.2e-01 1.0e+17 2.4e+05 367 max_eval update_μ Penalty 3.0: obj=0.3500752007510398, max_stress=1.2245689727078921
The continuation strategy (penalty 1.0 → 2.0 → 3.0) helps converge to feasible designs by first solving a relaxed problem, then progressively penalizing intermediate densities.
Verify Results
@show maximum(stress(filter(PseudoDensities(x0)))) # initial max stress
@show maximum(stress(filter(PseudoDensities(x)))) # final max stress
# Verify peak stress on solid elements doesn't exceed the threshold.
# Gray boundary-layer elements are penalized material whose true stress
# exceeds the relaxed measure by ρ^(-q) — a known artifact of filtered SIMP
# stress evaluation (da Silva et al., 2019) — so they are excluded here.
topo = filter(PseudoDensities(x)).x
s = stress(filter(PseudoDensities(x)))
max_solid = maximum(s[topo .> 0.5])
println("max stress on solid elements: $max_solid (threshold: $thr)")
println("stress constraint satisfied: $((max_solid - thr) / thr < 0.01)")maximum(stress(filter(PseudoDensities(x0)))) = 5.762395482668618
maximum(stress(filter(PseudoDensities(x)))) = 1.2245689727078921
max stress on solid elements: 1.2245689727078921 (threshold: 1.223183190956383)
stress constraint satisfied: trueVisualization
using CairoMakie
fig = visualize(
problem;
topology=topo,
default_exagg_scale=0.07,
scale_range=10.0,
vector_arrowsize=0.5,
)Precompiling packages... 6756.7 ms ✓ QuartoNotebookWorkerMakieExt (serial) 1 dependency successfully precompiled in 7 seconds Precompiling packages... 5531.1 ms ✓ QuartoNotebookWorkerCairoMakieExt (serial) 1 dependency successfully precompiled in 6 seconds
The optimized design shows smooth stress distributions without concentration points, demonstrating effective local stress control. See STRESS_CONSTRAINED_TO.md in the repository root for the theory and the comparison with global aggregation approaches.