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 ≤ σ_lim has 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

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)

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
end

stress 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.3e-01   2.8e+01   3.0e-01   1.0e+03   3.0e+01     295         max_eval         update_μ
[ Info:      4  -1.3e-01   1.4e+02   3.0e-01   1.0e+03   2.7e+02     299         max_eval         update_y
[ Info:      5  -1.3e-01   1.4e+02   3.0e-01   1.0e+04   2.7e+02     303         max_eval         update_μ
[ Info:      6  -1.3e-01   1.7e+03   3.0e-01   1.0e+04   3.3e+03     307         max_eval         update_y
[ Info:      7  -1.3e-01   1.7e+03   3.0e-01   1.0e+05   3.3e+03     311         max_eval         update_μ
[ Info:      8  -1.3e-01   1.8e+04   3.0e-01   1.0e+05   3.3e+04     315         max_eval         update_y
[ Info:      9  -1.3e-01   1.8e+04   3.0e-01   1.0e+06   3.3e+04     319         max_eval         update_μ
[ Info:     10  -1.3e-01   1.8e+04   3.0e-01   1.0e+07   3.3e+04     323         max_eval         update_μ
[ Info:     11  -1.3e-01   1.8e+04   3.0e-01   1.0e+08   3.3e+04     327         max_eval         update_μ
[ Info:     12  -1.3e-01   1.8e+04   3.0e-01   1.0e+09   3.3e+04     331         max_eval         update_μ
[ Info:     13  -1.3e-01   1.8e+04   3.0e-01   1.0e+10   3.3e+04     335         max_eval         update_μ
[ Info:     14  -1.3e-01   1.8e+04   3.0e-01   1.0e+11   3.3e+04     339         max_eval         update_μ
[ Info:     15  -1.3e-01   1.8e+04   3.0e-01   1.0e+12   3.3e+04     343         max_eval         update_μ
[ Info:     16  -1.3e-01   1.8e+04   3.0e-01   1.0e+13   3.3e+04     347         max_eval         update_μ
[ Info:     17  -1.3e-01   1.8e+04   3.0e-01   1.0e+14   3.3e+04     351         max_eval         update_μ
[ Info:     18  -1.3e-01   1.8e+04   3.0e-01   1.0e+15   3.3e+04     355         max_eval         update_μ
[ Info:     19  -1.4e-01   1.8e+04   3.0e-01   1.0e+16   3.3e+04     359         max_eval         update_μ
[ Info:     20  -1.9e-01   1.8e+04   3.0e-01   1.0e+17   3.3e+04     363         max_eval         update_μ
Penalty 1.0: obj=-0.1342641641983185, max_stress=1.2237054052900755
[ Info:   iter        fx    normgp    normcx         μ     normy    sumc     inner_status        iter_type  
[ Info:      0   4.9e-02   4.2e+01   4.3e+01   1.0e+01   4.0e+01       5
[ Info:      1   4.9e-02   4.2e+01   4.3e+01   1.0e+02   4.0e+01       9      first_order         update_μ
[ Info:      2   1.9e-01   5.5e+00   4.1e-01   1.0e+02   7.9e+00      58      first_order         update_y
[ Info:      3   1.7e-01   5.5e+00   7.8e-02   1.0e+03   7.9e+00     299         max_eval         update_μ
[ Info:      4   1.7e-01   6.2e+01   7.8e-02   1.0e+03   7.1e+01     303         max_eval         update_y
[ Info:      5   1.7e-01   6.2e+01   7.8e-02   1.0e+04   7.1e+01     307         max_eval         update_μ
[ Info:      6   1.7e-01   7.2e+02   7.8e-02   1.0e+04   8.5e+02     311         max_eval         update_y
[ Info:      7   1.7e-01   7.2e+02   7.8e-02   1.0e+05   8.5e+02     315         max_eval         update_μ
[ Info:      8   1.7e-01   7.3e+03   7.8e-02   1.0e+05   8.7e+03     319         max_eval         update_y
[ Info:      9   1.7e-01   7.3e+03   7.8e-02   1.0e+06   8.7e+03     323         max_eval         update_μ
[ Info:     10   1.7e-01   7.3e+04   7.8e-02   1.0e+06   8.7e+04     327         max_eval         update_y
[ Info:     11   1.7e-01   7.3e+04   7.8e-02   1.0e+07   8.7e+04     331         max_eval         update_μ
[ Info:     12   1.7e-01   7.3e+05   7.8e-02   1.0e+07   8.7e+05     335         max_eval         update_y
[ Info:     13   1.7e-01   7.3e+05   7.8e-02   1.0e+08   8.7e+05     339         max_eval         update_μ
[ Info:     14   1.7e-01   7.3e+06   7.8e-02   1.0e+08   8.7e+06     343         max_eval         update_y
[ Info:     15   1.7e-01   7.3e+06   7.8e-02   1.0e+09   8.7e+06     347         max_eval         update_μ
[ Info:     16   1.7e-01   7.3e+07   7.8e-02   1.0e+09   8.7e+07     351         max_eval         update_y
[ Info:     17   1.7e-01   7.3e+07   7.8e-02   1.0e+10   8.7e+07     355         max_eval         update_μ
[ Info:     18   1.7e-01   7.3e+08   7.8e-02   1.0e+10   8.7e+08     359         max_eval         update_y
[ Info:     19   1.7e-01   7.3e+08   7.8e-02   1.0e+11   8.7e+08     363         max_eval         update_μ
[ Info:     20   1.7e-01   7.3e+09   7.8e-02   1.0e+11   8.7e+09     367         max_eval         update_y
[ Info:     21   1.7e-01   7.3e+09   7.8e-02   1.0e+12   8.7e+09     371         max_eval         update_μ
Penalty 2.0: obj=0.17060264354513585, max_stress=1.2276419806046244
[ Info:   iter        fx    normgp    normcx         μ     normy    sumc     inner_status        iter_type  
[ Info:      0   2.0e-01   4.7e+01   1.6e+01   1.0e+01   4.0e+01       5
[ Info:      1   2.0e-01   4.7e+01   1.6e+01   1.0e+02   4.0e+01       9      first_order         update_μ
[ Info:      2   1.9e-01   1.2e+01   4.5e-01   1.0e+02   1.9e+01     111      first_order         update_y
[ Info:      3   2.0e-01   1.2e+01   1.9e-01   1.0e+03   1.9e+01     293         max_eval         update_μ
[ Info:      4   2.0e-01   1.5e+02   1.9e-01   1.0e+03   1.7e+02     297         max_eval         update_y
[ Info:      5   2.0e-01   1.5e+02   1.9e-01   1.0e+04   1.7e+02     301         max_eval         update_μ
[ Info:      6   2.0e-01   1.8e+03   1.9e-01   1.0e+04   2.0e+03     305         max_eval         update_y
[ Info:      7   2.0e-01   1.8e+03   1.9e-01   1.0e+05   2.0e+03     309         max_eval         update_μ
[ Info:      8   2.0e-01   1.8e+04   1.9e-01   1.0e+05   2.1e+04     313         max_eval         update_y
[ Info:      9   2.0e-01   1.8e+04   1.9e-01   1.0e+06   2.1e+04     317         max_eval         update_μ
[ Info:     10   2.0e-01   1.8e+05   1.9e-01   1.0e+06   2.1e+05     321         max_eval         update_y
[ Info:     11   2.0e-01   1.8e+05   1.9e-01   1.0e+07   2.1e+05     325         max_eval         update_μ
[ Info:     12   2.0e-01   1.8e+06   1.9e-01   1.0e+07   2.1e+06     329         max_eval         update_y
[ Info:     13   2.0e-01   1.8e+06   1.9e-01   1.0e+08   2.1e+06     333         max_eval         update_μ
[ Info:     14   2.0e-01   1.8e+06   1.9e-01   1.0e+09   2.1e+06     337         max_eval         update_μ
[ Info:     15   2.0e-01   1.8e+06   1.9e-01   1.0e+10   2.1e+06     341         max_eval         update_μ
[ Info:     16   2.0e-01   1.8e+06   1.9e-01   1.0e+11   2.1e+06     345         max_eval         update_μ
[ Info:     17   2.0e-01   1.8e+06   1.9e-01   1.0e+12   2.1e+06     349         max_eval         update_μ
[ Info:     18   2.0e-01   1.8e+06   1.9e-01   1.0e+13   2.1e+06     353         max_eval         update_μ
[ Info:     19   2.0e-01   1.8e+06   1.9e-01   1.0e+14   2.1e+06     357         max_eval         update_μ
[ Info:     20   2.0e-01   1.8e+06   1.9e-01   1.0e+15   2.1e+06     361         max_eval         update_μ
[ Info:     21   2.0e-01   1.8e+06   1.9e-01   1.0e+16   2.1e+06     365         max_eval         update_μ
Penalty 3.0: obj=0.19768495582484197, max_stress=1.2254852344858034

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.2254852344858034
max stress on solid elements: 1.2254852344858034 (threshold: 1.223183190956383)
stress constraint satisfied: true

Visualization

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

Figure 1: Local stress-constrained design with element-wise relaxed stress limits

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.