Skip to main content

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):

  1. 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.
  2. 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.2231831909518107
#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
Penalty 1.0: obj=-0.07082011911183865, max_stress=1.1772363812975966
Penalty 2.0: obj=0.057144999120655834, max_stress=1.236514617484557
Penalty 3.0: obj=0.23101868555045507, max_stress=1.2232352113815335

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.762395482626007
maximum(stress(filter(PseudoDensities(x)))) = 1.2232352113815335
max stress on solid elements: 1.2232352113815335 (threshold: 1.2231831909518107)
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.