using TopOpt, LinearAlgebra
using WGLMakie
WGLMakie.activate!(; resize_to=:parent)
using Bonito
if haskey(ENV, "QUARTO_PROJECT_DIR")
Bonito.Page(exportable=true, offline=true)
else
Bonito.browser_display()
end
display_app(app) = display(app)Local 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
WGLMakie.activate!(; resize_to=:parent) selects the browser renderer and fills the Quarto output column. Bonito.Page(exportable=true, offline=true) embeds the assets needed by visualize(...; static=true) in the Quarto output, so the visualization does not require a running Julia process.
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((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.2231831909556796#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.4e-01 2.8e+01 3.0e-01 1.0e+03 3.0e+01 295 max_eval update_μ [ Info: 4 -1.4e-01 1.4e+02 3.0e-01 1.0e+03 2.7e+02 299 max_eval update_y [ Info: 5 -1.4e-01 1.4e+02 3.0e-01 1.0e+04 2.7e+02 303 max_eval update_μ [ Info: 6 -1.4e-01 1.7e+03 3.0e-01 1.0e+04 3.3e+03 307 max_eval update_y [ Info: 7 -1.4e-01 1.7e+03 3.0e-01 1.0e+05 3.3e+03 311 max_eval update_μ [ Info: 8 -1.4e-01 1.7e+04 3.0e-01 1.0e+05 3.3e+04 315 max_eval update_y [ Info: 9 -1.4e-01 1.7e+04 3.0e-01 1.0e+06 3.3e+04 319 max_eval update_μ [ Info: 10 -1.4e-01 1.7e+04 3.0e-01 1.0e+07 3.3e+04 323 max_eval update_μ [ Info: 11 -1.4e-01 1.7e+04 3.0e-01 1.0e+08 3.3e+04 327 max_eval update_μ [ Info: 12 -1.4e-01 1.7e+04 3.0e-01 1.0e+09 3.3e+04 331 max_eval update_μ [ Info: 13 -1.4e-01 1.7e+04 3.0e-01 1.0e+10 3.3e+04 335 max_eval update_μ [ Info: 14 -1.4e-01 1.7e+04 3.0e-01 1.0e+11 3.3e+04 339 max_eval update_μ [ Info: 15 -1.4e-01 1.7e+04 3.0e-01 1.0e+12 3.3e+04 343 max_eval update_μ [ Info: 16 -1.4e-01 1.7e+04 3.0e-01 1.0e+13 3.3e+04 347 max_eval update_μ [ Info: 17 -1.4e-01 1.7e+04 3.0e-01 1.0e+14 3.3e+04 351 max_eval update_μ [ Info: 18 -1.4e-01 1.7e+04 3.0e-01 1.0e+15 3.3e+04 355 max_eval update_μ [ Info: 19 -1.2e-01 1.7e+04 3.0e-01 1.0e+16 3.3e+04 359 max_eval update_μ [ Info: 20 -1.2e-01 1.7e+04 3.0e-01 1.0e+17 3.3e+04 363 max_eval update_μ Penalty 1.0: obj=-0.13913102745464836, max_stress=1.21891626182133 [ Info: iter fx normgp normcx μ normy sumc inner_status iter_type [ Info: 0 5.1e-02 4.2e+01 4.5e+01 1.0e+01 4.0e+01 5 [ Info: 1 5.1e-02 4.2e+01 4.5e+01 1.0e+02 4.0e+01 9 first_order update_μ [ Info: 2 2.1e-01 5.2e+00 4.0e-01 1.0e+02 7.4e+00 56 first_order update_y [ Info: 3 1.9e-01 5.1e+00 7.4e-02 1.0e+03 7.4e+00 303 max_eval update_μ [ Info: 4 1.9e-01 6.0e+01 7.4e-02 1.0e+03 6.6e+01 307 max_eval update_y [ Info: 5 1.9e-01 6.0e+01 7.4e-02 1.0e+04 6.6e+01 311 max_eval update_μ [ Info: 6 1.9e-01 6.9e+02 7.4e-02 1.0e+04 8.0e+02 315 max_eval update_y [ Info: 7 1.9e-01 6.9e+02 7.4e-02 1.0e+05 8.0e+02 319 max_eval update_μ [ Info: 8 1.9e-01 7.0e+03 7.4e-02 1.0e+05 8.2e+03 323 max_eval update_y [ Info: 9 1.9e-01 7.0e+03 7.4e-02 1.0e+06 8.2e+03 327 max_eval update_μ [ Info: 10 1.9e-01 7.0e+04 7.4e-02 1.0e+06 8.2e+04 331 max_eval update_y [ Info: 11 1.9e-01 7.0e+04 7.4e-02 1.0e+07 8.2e+04 335 max_eval update_μ [ Info: 12 1.9e-01 7.0e+05 7.4e-02 1.0e+07 8.2e+05 339 max_eval update_y [ Info: 13 1.9e-01 7.0e+05 7.4e-02 1.0e+08 8.2e+05 343 max_eval update_μ [ Info: 14 1.9e-01 7.0e+06 7.4e-02 1.0e+08 8.2e+06 347 max_eval update_y [ Info: 15 1.9e-01 7.0e+06 7.4e-02 1.0e+09 8.2e+06 351 max_eval update_μ [ Info: 16 1.9e-01 7.0e+07 7.4e-02 1.0e+09 8.2e+07 355 max_eval update_y [ Info: 17 1.9e-01 7.0e+07 7.4e-02 1.0e+10 8.2e+07 359 max_eval update_μ [ Info: 18 1.9e-01 7.0e+08 7.4e-02 1.0e+10 8.2e+08 363 max_eval update_y [ Info: 19 1.9e-01 7.0e+08 7.4e-02 1.0e+11 8.2e+08 367 max_eval update_μ [ Info: 20 1.9e-01 7.0e+09 7.4e-02 1.0e+11 8.2e+09 371 max_eval update_y [ Info: 21 1.9e-01 7.0e+09 7.4e-02 1.0e+12 8.2e+09 375 max_eval update_μ Penalty 2.0: obj=0.19444441497492118, max_stress=1.2185104537699378 [ Info: iter fx normgp normcx μ normy sumc inner_status iter_type [ Info: 0 2.2e-01 4.8e+01 1.7e+01 1.0e+01 4.0e+01 5 [ Info: 1 2.2e-01 4.8e+01 1.7e+01 1.0e+02 4.0e+01 9 first_order update_μ [ Info: 2 1.9e-01 1.4e+01 4.6e-01 1.0e+02 2.1e+01 116 first_order update_y [ Info: 3 2.0e-01 1.4e+01 2.1e-01 1.0e+03 2.1e+01 291 max_eval update_μ [ Info: 4 2.0e-01 1.7e+02 2.1e-01 1.0e+03 1.9e+02 295 max_eval update_y [ Info: 5 2.0e-01 1.7e+02 2.1e-01 1.0e+04 1.9e+02 299 max_eval update_μ [ Info: 6 2.0e-01 2.1e+03 2.1e-01 1.0e+04 2.3e+03 303 max_eval update_y [ Info: 7 2.0e-01 2.1e+03 2.1e-01 1.0e+05 2.3e+03 307 max_eval update_μ [ Info: 8 2.0e-01 2.1e+04 2.1e-01 1.0e+05 2.3e+04 311 max_eval update_y [ Info: 9 2.0e-01 2.1e+04 2.1e-01 1.0e+06 2.3e+04 315 max_eval update_μ [ Info: 10 2.0e-01 2.1e+05 2.1e-01 1.0e+06 2.3e+05 319 max_eval update_y [ Info: 11 2.0e-01 2.1e+05 2.1e-01 1.0e+07 2.3e+05 323 max_eval update_μ [ Info: 12 2.0e-01 2.1e+05 2.1e-01 1.0e+08 2.3e+05 327 max_eval update_μ [ Info: 13 2.0e-01 2.1e+05 2.1e-01 1.0e+09 2.3e+05 331 max_eval update_μ [ Info: 14 2.0e-01 2.1e+05 2.1e-01 1.0e+10 2.3e+05 335 max_eval update_μ [ Info: 15 2.0e-01 2.1e+05 2.1e-01 1.0e+11 2.3e+05 339 max_eval update_μ [ Info: 16 2.0e-01 2.1e+05 2.1e-01 1.0e+12 2.3e+05 343 max_eval update_μ [ Info: 17 2.0e-01 2.1e+05 2.1e-01 1.0e+13 2.3e+05 347 max_eval update_μ [ Info: 18 2.0e-01 2.1e+05 2.1e-01 1.0e+14 2.3e+05 351 max_eval update_μ [ Info: 19 2.0e-01 2.1e+05 2.1e-01 1.0e+15 2.3e+05 355 max_eval update_μ [ Info: 20 2.0e-01 2.1e+05 2.1e-01 1.0e+16 2.3e+05 359 max_eval update_μ [ Info: 21 1.9e-01 2.1e+05 2.1e-01 1.0e+17 2.3e+05 363 max_eval update_μ Penalty 3.0: obj=0.1962924096903936, max_stress=1.2258160297395517
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.76239548267132
maximum(stress(filter(PseudoDensities(x)))) = 1.2258160297395517
max stress on solid elements: 1.2258160297395517 (threshold: 1.2231831909556796)
stress constraint satisfied: trueVisualization
fig = visualize(
problem;
static=true,
topology=topo,
default_exagg_scale=0.07,
scale_range=10.0,
vector_arrowsize=0.5,
)
display_app(fig)undeformed mesh
load arrows
support arrows
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.
Hyper-parameter guide
The reliable defaults and their ranges:
| Parameter | Default here | Range | Effect of increasing it |
|---|---|---|---|
stress_exponent (q) |
0.5 | 0.3–0.5 | weaker relaxation → singular optima |
SIMP p |
continuation 1→3 | 1→3 | skip the p=1 warm start if the design is already feasible |
stress threshold thr |
1.2× solid peak | 1.1–1.5× solid peak | more feasible, higher volume |
compliance weight λ |
1e-4 | 1e-5–1e-3 | stronger regularizer → stiffer but less volume-reduced designs |
filter rmin |
3.0 | 1.5–3 elements | smoother designs |
Percival max_iter |
20 | 10–40 | more outer AL iterations |
Percival subsolver_max_eval |
100 | 50–200 | more accurate subproblem solves |
As in the global case, the threshold must be calibrated from the solid design’s peak stress — an arbitrary value is either vacuous or infeasible, and the latter stalls the optimizer at the solid design.