From 5bb1c1fc9b19a619a6fa130c3dbb02404df09eea Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 27 Aug 2026 15:04:09 -0700 Subject: [PATCH 1/3] Add minimum-lengthscale constraints to the Python package Ports the geometric lengthscale constraints from the Julia package: constraint_solid and constraint_void, with the hyperparameter-free thresholds of Arrieta et al. (arXiv:2507.16108). Both are written in terms of rho_filtered and rho_projected only, so they compose with any projection order (ssp1_bilinear, ssp2, or a plain tanh projection). Includes a test suite covering the threshold functions, lengthscale detection on stripe/gap geometries, gradients against finite differences, composition with each projection, jit, and argument validation, plus a worked example that sweeps the constraint over feature widths and runs a two-stage constrained optimization. --- README.md | 1 + docs/features.md | 6 +- examples/python/ssp_constrained_example.py | 269 ++++++++++++++++ src/python/ssp_topopt/__init__.py | 3 + src/python/ssp_topopt/constraints.py | 271 ++++++++++++++++ tests/python/test_constraints.py | 352 +++++++++++++++++++++ 6 files changed, 899 insertions(+), 3 deletions(-) create mode 100644 examples/python/ssp_constrained_example.py create mode 100644 src/python/ssp_topopt/constraints.py create mode 100644 tests/python/test_constraints.py diff --git a/README.md b/README.md index f61f7af..f5dddb1 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,7 @@ The Python import package is `ssp_topopt`: ```python from ssp_topopt import conic_filter, get_conic_radius_from_eta_e, ssp1_bilinear,ssp2 +from ssp_topopt import constraint_solid, constraint_void ``` For local development: diff --git a/docs/features.md b/docs/features.md index a8c2266..c93a153 100644 --- a/docs/features.md +++ b/docs/features.md @@ -14,8 +14,8 @@ currently implements; **please keep it up to date when adding or removing functi | Second-order subpixel smoothing (SSP2), differentiable through topology changes | ✅ `ssp2` | ✅ `ssp2` | | Finite and infinite projection strength (0 ≤ β ≤ ∞) | ✅ | ✅ | | Dilation/erosion of the projected contour | ✅ `dilation_distance` argument | ❌ | -| Minimum-lengthscale constraints for solid and void | ✅ `constraint_solid`, `constraint_void` | ❌ | -| Lengthscale constraints compatible with any SSP order | ✅ (constraints act on `rho_filtered`/`rho_projected`) | ❌ | +| Minimum-lengthscale constraints for solid and void | ✅ `constraint_solid`, `constraint_void` | ✅ `constraint_solid`, `constraint_void` | +| Lengthscale constraints compatible with any SSP order | ✅ (constraints act on `rho_filtered`/`rho_projected`) | ✅ (constraints act on `rho_filtered`/`rho_projected`) | | Reverse-mode automatic differentiation | ✅ hand-written adjoints, exposed to Zygote.jl and friends through a ChainRulesCore.jl extension | ✅ through JAX (`grad`, `jit`, `vmap`) | | Dimensionality | N-dimensional code paths (only 2D is currently tested) | 2D only | | Periodic filter axes | ❌ | ✅ `periodic_axes` argument of `conic_filter` | @@ -37,7 +37,7 @@ currently implements; **please keep it up to date when adding or removing functi Contributions welcome — these are known gaps rather than fundamental limitations: -* Python: cubic-interpolation SSP1, dilation/erosion, and minimum-lengthscale constraints. +* Python: cubic-interpolation SSP1 and dilation/erosion. * Python: a low-level API with reusable workspaces. * Julia: `get_conic_radius_from_eta_e`-style helpers and periodic filter axes. * Both: validated 3D usage and anisotropic grid spacings in the projection. diff --git a/examples/python/ssp_constrained_example.py b/examples/python/ssp_constrained_example.py new file mode 100644 index 0000000..075a82a --- /dev/null +++ b/examples/python/ssp_constrained_example.py @@ -0,0 +1,269 @@ +"""Minimum-lengthscale (geometric) constraints combined with SSP. + +The constraints of `ssp_topopt.constraints` measure whether the solid or void +features of a design are thinner than a target lengthscale. They are formulated +in terms of the filtered and projected densities only, so they work with any +order of subpixel smoothing; this example uses `ssp2`. + +The example has two parts: + +1. A sweep over a simple stripe geometry, which shows that the solid constraint + changes sign at the target lengthscale (compare fig. 1c of Arrieta et al.). +2. A two-stage optimization that matches a target pattern containing bars that + are thinner than the target lengthscale. The first stage is unconstrained and + happily reproduces the sub-lengthscale bars; the second stage turns the + constraints on and repairs the design, at a small cost in the figure of merit. + +Ref: R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free +minimum-lengthscale constraints for topology optimization," arXiv.org e-Print +archive, 2507.16108, July 2025. +""" + +import time + +import nlopt +import numpy as np +from jax import jit, value_and_grad +from jax import numpy as jnp +from matplotlib import pyplot as plt + +from ssp_topopt import conic_filter, constraint_solid, constraint_void, ssp2 + +BETA = np.inf +ETA_I = 0.5 + + +def filter_and_project(rho, filter_radius, lx, ly, resolution): + """The usual filter-then-project pipeline.""" + rho_filtered = conic_filter(rho, filter_radius, lx, ly, resolution) + rho_projected = ssp2(rho_filtered, BETA, ETA_I, resolution) + return rho_filtered, rho_projected + + +def feature_width(rho_projected, resolution, solid=True): + """Width of the central feature of a stripe design, in physical units.""" + profile = np.asarray(rho_projected)[:, rho_projected.shape[1] // 2] + mask = profile > 0.5 if solid else profile < 0.5 + return float(np.count_nonzero(mask)) / resolution + + +def constraint_sweep(): + """Evaluate the constraints for stripes of varying width.""" + lx = ly = 0.5 + resolution = 100 + target_length = 0.1 + filter_radius = target_length + + nx = int(np.round(lx * resolution)) + 1 + ny = int(np.round(ly * resolution)) + 1 + coords = np.linspace(-lx / 2, lx / 2, nx) + x = coords[:, None] * np.ones((1, ny)) + + latent_widths = np.linspace(0.04, 0.24, 21) + widths = [] + solid_values = [] + void_values = [] + + # Note that the void constraint also fires for very thin stripes: once the + # stripe is thin enough that it barely projects to solid, the crest of the + # filtered density is a nearly-solid feature sitting inside the void region. + print("Constraint sweep over stripe widths " f"(target lengthscale {target_length})") + for latent_width in latent_widths: + # A binary stripe of width `latent_width`, which the filter and the + # projection turn into a stripe of some (smaller) physical width. + rho = jnp.asarray((np.abs(x) <= latent_width / 2).astype(float)) + rho_filtered, rho_projected = filter_and_project( + rho, filter_radius, lx, ly, resolution + ) + + width = feature_width(rho_projected, resolution) + solid_value = float( + constraint_solid(rho_filtered, rho_projected, resolution, target_length) + ) + void_value = float( + constraint_void(rho_filtered, rho_projected, resolution, target_length) + ) + + widths.append(width) + solid_values.append(solid_value) + void_values.append(void_value) + print( + f" latent width={latent_width:.3f} physical width={width:.3f} " + f"solid={solid_value:+.3e} void={void_value:+.3e}" + ) + + # Feasible designs have an exactly vanishing violation, so clip to a floor to + # keep them visible on a log scale. + floor = 1e-2 + plt.figure(figsize=(5, 3.5)) + plt.axhspan(floor / 2, 1.0, color="tab:green", alpha=0.1, label="feasible") + plt.plot( + widths, np.maximum(np.asarray(solid_values) + 1, floor), marker="o", label="solid" + ) + plt.plot( + widths, np.maximum(np.asarray(void_values) + 1, floor), marker="s", label="void" + ) + plt.axhline(1.0, color="k", linestyle=":", label="threshold") + plt.axvline(target_length, color="r", linestyle="--", label="target lengthscale") + plt.yscale("log") + plt.ylim(bottom=floor / 2) + plt.xlabel("physical width of the stripe") + plt.ylabel("constraint / threshold") + plt.legend(fontsize=8) + plt.title("Constraint value vs feature width") + plt.tight_layout() + plt.savefig("constraint_sweep.png") + + +def optimization_demo(): + """Repair a design with sub-lengthscale features using the constraints.""" + lx, ly = 1.2, 0.6 + resolution = 50 + target_length = 0.12 + filter_radius = target_length + + nx = int(np.round(lx * resolution)) + 1 + ny = int(np.round(ly * resolution)) + 1 + num_vars = nx * ny + coords = np.linspace(-lx / 2, lx / 2, nx) + + # A target pattern of vertical bars, two of which are thinner than the + # target lengthscale and therefore cannot be manufactured. + target = np.zeros((nx, ny)) + position = -lx / 2 + 0.08 + for width in (0.24, 0.16, 0.08, 0.04): + target[(coords >= position) & (coords <= position + width), :] = 1.0 + position += width + 0.12 + target_jnp = jnp.asarray(target) + + def pipeline(rho_flat): + return filter_and_project( + rho_flat.reshape((nx, ny)), filter_radius, lx, ly, resolution + ) + + def figure_of_merit(rho_flat): + _, rho_projected = pipeline(rho_flat) + return jnp.mean((rho_projected - target_jnp) ** 2) + + def solid_constraint(rho_flat): + rho_filtered, rho_projected = pipeline(rho_flat) + return constraint_solid( + rho_filtered, rho_projected, resolution, target_length + ) + + def void_constraint(rho_flat): + rho_filtered, rho_projected = pipeline(rho_flat) + return constraint_void(rho_filtered, rho_projected, resolution, target_length) + + objective_and_grad = jit(value_and_grad(figure_of_merit)) + constraints_and_grad = { + "solid": jit(value_and_grad(solid_constraint)), + "void": jit(value_and_grad(void_constraint)), + } + + fom_history = [] + constraint_history = {"solid": [], "void": []} + + def run_stage(x_init, constrained, maxeval): + def nlopt_objective(x, grad_out): + value, gradient = objective_and_grad(jnp.asarray(x)) + if grad_out.size > 0: + grad_out[:] = np.asarray(gradient, dtype=float) + fom_history.append(float(value)) + for name, value_and_grad_fn in constraints_and_grad.items(): + constraint_history[name].append( + float(value_and_grad_fn(jnp.asarray(x))[0]) + ) + return float(value) + + opt = nlopt.opt(nlopt.LD_CCSAQ, num_vars) + opt.set_lower_bounds(np.zeros(num_vars)) + opt.set_upper_bounds(np.ones(num_vars)) + opt.set_min_objective(nlopt_objective) + + if constrained: + for value_and_grad_fn in constraints_and_grad.values(): + + def nlopt_constraint(x, grad_out, fn=value_and_grad_fn): + value, gradient = fn(jnp.asarray(x)) + if grad_out.size > 0: + grad_out[:] = np.asarray(gradient, dtype=float) + return float(value) + + # The constraints are normalized, so the feasible region is + # exactly where they are nonpositive. + opt.add_inequality_constraint(nlopt_constraint, 0.0) + + opt.set_maxeval(maxeval) + start = time.perf_counter() + x_opt = opt.optimize(np.asarray(x_init, dtype=float)) + elapsed = time.perf_counter() - start + + label = "constrained" if constrained else "unconstrained" + print( + f" {label}: {maxeval} evaluations in {elapsed:6.1f}s " + f"FOM={fom_history[-1]:.4e} " + f"solid={constraint_history['solid'][-1]:+.3e} " + f"void={constraint_history['void'][-1]:+.3e}" + ) + return x_opt + + print(f"\nTwo-stage optimization (target lengthscale {target_length})") + rng = np.random.default_rng(0) + x_init = 0.5 * np.ones(num_vars) + 0.01 * rng.standard_normal(num_vars) + + stage1_evals = 120 + x_stage1 = run_stage(x_init, constrained=False, maxeval=stage1_evals) + # The constrained stage needs a few hundred iterations: CCSA has to first + # drag the design back into the feasible region and then re-minimize. + x_stage2 = run_stage(x_stage1, constrained=True, maxeval=400) + + _, projected_stage1 = pipeline(jnp.asarray(x_stage1)) + _, projected_stage2 = pipeline(jnp.asarray(x_stage2)) + + plt.figure(figsize=(9, 3)) + for index, (image, title) in enumerate( + ( + (target, "target"), + (np.asarray(projected_stage1), "stage 1: unconstrained"), + (np.asarray(projected_stage2), "stage 2: lengthscale constrained"), + ) + ): + plt.subplot(1, 3, index + 1) + plt.imshow(image.T, vmin=0, vmax=1, cmap="binary", origin="lower") + plt.title(title, fontsize=9) + plt.axis("off") + plt.tight_layout() + plt.savefig("constraint_designs.png") + + plt.figure(figsize=(8, 3)) + plt.subplot(1, 2, 1) + plt.semilogy(fom_history, linewidth=1.5) + plt.axvline(stage1_evals, color="k", linestyle="--", linewidth=1) + plt.xlabel("iteration") + plt.ylabel("FOM") + plt.title("Objective history", fontsize=9) + + plt.subplot(1, 2, 2) + for name, values in constraint_history.items(): + # Shift by one so that the (normalized) constraint can be shown on a log + # scale: values below one are feasible. + plt.semilogy(np.asarray(values) + 1, linewidth=1.5, label=name) + plt.axhline(1.0, color="k", linestyle=":", label="threshold") + plt.axvline(stage1_evals, color="k", linestyle="--", linewidth=1) + plt.xlabel("iteration") + plt.ylabel("constraint / threshold") + plt.title("Constraint history", fontsize=9) + plt.legend(fontsize=8) + plt.tight_layout() + plt.savefig("constraint_history.png") + plt.show() + + +def main(): + constraint_sweep() + optimization_demo() + + +if __name__ == "__main__": + main() diff --git a/src/python/ssp_topopt/__init__.py b/src/python/ssp_topopt/__init__.py index e58537f..d0f57f7 100644 --- a/src/python/ssp_topopt/__init__.py +++ b/src/python/ssp_topopt/__init__.py @@ -1,5 +1,6 @@ """Public Python API for smoothed subpixel projection (SSP) for topology optimization.""" +from .constraints import constraint_solid, constraint_void from .core import ssp1_bilinear,ssp2 from .utils import conic_filter, get_conic_radius_from_eta_e, tanh_projection @@ -7,6 +8,8 @@ "ssp1_bilinear", "ssp2", "conic_filter", + "constraint_solid", + "constraint_void", "get_conic_radius_from_eta_e", "tanh_projection", ] diff --git a/src/python/ssp_topopt/constraints.py b/src/python/ssp_topopt/constraints.py new file mode 100644 index 0000000..8a98757 --- /dev/null +++ b/src/python/ssp_topopt/constraints.py @@ -0,0 +1,271 @@ +"""Minimum-lengthscale (geometric) constraints for topology optimization. + +These are the geometric constraints of Zhou et al. (2015) combined with the +hyperparameter-free thresholds derived by Arrieta et al. (2025). The constraints +penalize solid (or void) features whose lengthscale falls below a target value, +and are formulated purely in terms of the filtered density and the projected +density. Consequently they work with *any* order of subpixel smoothing -- +`ssp1_bilinear`, `ssp2`, or even a plain `tanh_projection` -- since the +projection only enters through `rho_projected`. + +The solid constraint reads + + g_s = (1/N) Σ_i I_s,i [min(rho_filtered_i - eta_e, 0)]^2 , + I_s = rho_projected * exp(-c |∇ rho_filtered|^2) , + +and the void constraint is the complementary expression + + g_v = (1/N) Σ_i I_v,i [min(eta_d - rho_filtered_i, 0)]^2 , + I_v = (1 - rho_projected) * exp(-c |∇ rho_filtered|^2) , + +where the structural functions I_s and I_v single out the "inflection regions" +of the design (the interior of a feature, where the filtered density is +stationary), and eta_e, eta_d are the eroded/dilated threshold points of the +conic filter. Following Arrieta et al., the decay rate is c = 64 R^2 and the +constraint threshold is eps = 1e-8, where R is the conic filter radius; the +constraint is well behaved for target_length / R roughly in [0.25, 1.5]. + +Refs: + +R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free minimum-lengthscale +constraints for topology optimization," arXiv.org e-Print archive, 2507.16108, +July 2025. + +M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, "Minimum length scale in topology +optimization by geometric constraints," Computer Methods in Applied Mechanics and +Engineering, vol. 293, pp. 266-282, 2015. + +X. Qian and O. Sigmund, "Topological design of electromechanical actuators with +robustness toward over- and under-etching," Computer Methods in Applied Mechanics +and Engineering, vol. 253, pp. 237-251, 2013. +""" + +from typing import Optional + +from jax import numpy as jnp + +from .utils import ArrayLikeType, gradient + +# Hyperparameters derived in section 4.1 of Arrieta et al. (2025). The decay rate +# is expressed relative to the square of the conic filter radius, i.e. the actual +# decay rate is c = DEFAULT_CONSTRAINT_DECAYRATE * conic_radius**2. +DEFAULT_CONSTRAINT_THRESHOLD = 1e-8 +DEFAULT_CONSTRAINT_DECAYRATE = 64.0 + + +def solid_threshold(lengthscale_ratio: float): + """The eroded threshold point eta_e of a conic filter. + + Ref: Eq. (9) of Arrieta et al. (2025), originally from Qian and Sigmund (2013). + + Args: + lengthscale_ratio: the ratio of the target lengthscale to the conic + filter radius, which must be nonnegative. + + Returns: + The threshold point in the range [1/2, 1]. + """ + x = jnp.asarray(lengthscale_ratio) + return jnp.where( + x < 1, + x**2 / 4 + 1 / 2, + jnp.where(x < 2, -(x**2) / 4 + x, 1.0), + ) + + +def void_threshold(lengthscale_ratio: float): + """The dilated threshold point eta_d of a conic filter. + + Ref: Eq. (12) of Arrieta et al. (2025), originally from Qian and Sigmund (2013). + + Args: + lengthscale_ratio: the ratio of the target lengthscale to the conic + filter radius, which must be nonnegative. + + Returns: + The threshold point in the range [0, 1/2]. + """ + x = jnp.asarray(lengthscale_ratio) + return jnp.where( + x < 1, + 1 / 2 - x**2 / 4, + jnp.where(x < 2, 1 + x**2 / 4 - x, 0.0), + ) + + +def _geometric_constraint( + rho_filtered: ArrayLikeType, + rho_projected: ArrayLikeType, + resolution: float, + target_length: float, + conic_radius: Optional[float], + constraint_threshold: float, + constraint_decayrate: float, + solid: bool, +): + """Shared implementation of the solid and void lengthscale constraints.""" + if target_length < 0: + raise ValueError("The target lengthscale must be nonnegative.") + + if conic_radius is None: + conic_radius = target_length + if conic_radius <= 0: + raise ValueError("The conic filter radius must be positive.") + + rho_filtered = jnp.asarray(rho_filtered) + rho_projected = jnp.asarray(rho_projected) + if rho_filtered.ndim != 2: + raise ValueError( + f"Only 2D designs are supported, got {rho_filtered.ndim} dimensions." + ) + if rho_filtered.shape != rho_projected.shape: + raise ValueError( + "rho_filtered and rho_projected must have the same shape, got " + f"{rho_filtered.shape} and {rho_projected.shape}." + ) + + decayrate = constraint_decayrate * conic_radius**2 + + # The gradient of the filtered density vanishes in the interior of a feature, + # so this term restricts the constraint to those "inflection regions". + rho_filtered_grad = gradient(rho_filtered, resolution) + rho_filtered_grad_normsq = jnp.sum(rho_filtered_grad**2, axis=-1) + extremal_region = jnp.exp(-decayrate * rho_filtered_grad_normsq) + + if solid: + eta_m = solid_threshold(target_length / conic_radius) + inflection_region = rho_projected * extremal_region + # Only densities below the eroded threshold, i.e. features that are too + # thin to survive an erosion by the target lengthscale, are penalized. + beyond_threshold = jnp.minimum(rho_filtered - eta_m, 0.0) + else: + eta_m = void_threshold(target_length / conic_radius) + inflection_region = (1 - rho_projected) * extremal_region + beyond_threshold = jnp.minimum(eta_m - rho_filtered, 0.0) + + violation = jnp.mean(inflection_region * beyond_threshold**2) + + # Normalize so that the constraint is satisfied when it is nonpositive. + return violation / constraint_threshold - 1 + + +def constraint_solid( + rho_filtered: ArrayLikeType, + rho_projected: ArrayLikeType, + resolution: float, + target_length: float, + conic_radius: Optional[float] = None, + constraint_threshold: float = DEFAULT_CONSTRAINT_THRESHOLD, + constraint_decayrate: float = DEFAULT_CONSTRAINT_DECAYRATE, +): + """Calculate a solid minimum-lengthscale constraint function. + + This technique takes smoothed data, e.g. from filtering, `rho_filtered` and + binary data, e.g. from projection, `rho_projected` both defined on the same + grid, and measures whether features in the solid region, i.e. where + `rho_projected` takes values of 1, violate the minimum `target_length`. + + The returned value is normalized by the constraint threshold, so the + constraint is satisfied when the value is nonpositive and violated when it is + positive. It can therefore be handed directly to a nonlinear optimizer such + as `nlopt` as an inequality constraint. The unnormalized constraint value of + Eq. (7) of Arrieta et al. (2025) is `(value + 1) * constraint_threshold`. + + Any projection may be used to produce `rho_projected`, including + `ssp1_bilinear`, `ssp2`, and `tanh_projection`. + + Args: + rho_filtered: the (2D) filtered design parameters, e.g. from + `conic_filter`. + rho_projected: the (2D) projected design parameters, e.g. from `ssp2`. + resolution: resolution of the design grid. + target_length: the minimum lengthscale to impose on the solid region. + conic_radius: the radius of the conic filter used to obtain + `rho_filtered`. Defaults to `target_length`, which is the + recommended choice. + constraint_threshold: the threshold that separates feasible from + infeasible designs. May be tuned if feasible designs still don't + meet the target lengthscale. + constraint_decayrate: the decay rate of the structural function outside + of the inflection region, relative to `conic_radius**2`. + + Returns: + The normalized constraint value, which is nonpositive when the design + satisfies the minimum lengthscale. + + Example: + >>> rho_filtered = conic_filter(rho, filter_radius, lx, ly, resolution) + >>> rho_projected = ssp2(rho_filtered, beta, eta_i, resolution) + >>> constraint_solid(rho_filtered, rho_projected, resolution, filter_radius) + """ + return _geometric_constraint( + rho_filtered, + rho_projected, + resolution, + target_length, + conic_radius, + constraint_threshold, + constraint_decayrate, + solid=True, + ) + + +def constraint_void( + rho_filtered: ArrayLikeType, + rho_projected: ArrayLikeType, + resolution: float, + target_length: float, + conic_radius: Optional[float] = None, + constraint_threshold: float = DEFAULT_CONSTRAINT_THRESHOLD, + constraint_decayrate: float = DEFAULT_CONSTRAINT_DECAYRATE, +): + """Calculate a void minimum-lengthscale constraint function. + + This technique takes smoothed data, e.g. from filtering, `rho_filtered` and + binary data, e.g. from projection, `rho_projected` both defined on the same + grid, and measures whether features in the void region, i.e. where + `rho_projected` takes values of 0, violate the minimum `target_length`. + + The returned value is normalized by the constraint threshold, so the + constraint is satisfied when the value is nonpositive and violated when it is + positive. It can therefore be handed directly to a nonlinear optimizer such + as `nlopt` as an inequality constraint. The unnormalized constraint value of + Eq. (10) of Arrieta et al. (2025) is `(value + 1) * constraint_threshold`. + + Any projection may be used to produce `rho_projected`, including + `ssp1_bilinear`, `ssp2`, and `tanh_projection`. + + Args: + rho_filtered: the (2D) filtered design parameters, e.g. from + `conic_filter`. + rho_projected: the (2D) projected design parameters, e.g. from `ssp2`. + resolution: resolution of the design grid. + target_length: the minimum lengthscale to impose on the void region. + conic_radius: the radius of the conic filter used to obtain + `rho_filtered`. Defaults to `target_length`, which is the + recommended choice. + constraint_threshold: the threshold that separates feasible from + infeasible designs. May be tuned if feasible designs still don't + meet the target lengthscale. + constraint_decayrate: the decay rate of the structural function outside + of the inflection region, relative to `conic_radius**2`. + + Returns: + The normalized constraint value, which is nonpositive when the design + satisfies the minimum lengthscale. + + Example: + >>> rho_filtered = conic_filter(rho, filter_radius, lx, ly, resolution) + >>> rho_projected = ssp2(rho_filtered, beta, eta_i, resolution) + >>> constraint_void(rho_filtered, rho_projected, resolution, filter_radius) + """ + return _geometric_constraint( + rho_filtered, + rho_projected, + resolution, + target_length, + conic_radius, + constraint_threshold, + constraint_decayrate, + solid=False, + ) diff --git a/tests/python/test_constraints.py b/tests/python/test_constraints.py new file mode 100644 index 0000000..db5cc6d --- /dev/null +++ b/tests/python/test_constraints.py @@ -0,0 +1,352 @@ +"""Tests for the minimum-lengthscale (geometric) constraints. + +The reference behavior is the one described in R. Arrieta, G. Romano, and +S. G. Johnson, "Hyperparameter-free minimum-lengthscale constraints for topology +optimization," arXiv:2507.16108 (2025): with the derived hyperparameters, the +constraint is violated (positive) when the physical lengthscale of a feature +falls below the target lengthscale and satisfied (nonpositive) otherwise. +""" + +import unittest + +import numpy as np +from jax import grad, jit, value_and_grad +from jax import numpy as jnp +from jax.experimental import enable_x64 + +from ssp_topopt import ( + conic_filter, + constraint_solid, + constraint_void, + ssp1_bilinear, + ssp2, + tanh_projection, +) +from ssp_topopt.constraints import solid_threshold, void_threshold + + +class TestThresholdFunctions(unittest.TestCase): + """The conic-filter threshold points of Qian and Sigmund (2013).""" + + def test_known_values(self): + # Eqs. (9) and (12) of Arrieta et al. (2025). + self.assertAlmostEqual(float(solid_threshold(0.0)), 0.5) + self.assertAlmostEqual(float(solid_threshold(1.0)), 0.75) + self.assertAlmostEqual(float(solid_threshold(2.0)), 1.0) + self.assertAlmostEqual(float(solid_threshold(3.0)), 1.0) + + self.assertAlmostEqual(float(void_threshold(0.0)), 0.5) + self.assertAlmostEqual(float(void_threshold(1.0)), 0.25) + self.assertAlmostEqual(float(void_threshold(2.0)), 0.0) + self.assertAlmostEqual(float(void_threshold(3.0)), 0.0) + + def test_thresholds_are_complementary_and_monotonic(self): + ratios = np.linspace(0.0, 3.0, 61) + eta_e = np.asarray(solid_threshold(ratios)) + eta_d = np.asarray(void_threshold(ratios)) + + np.testing.assert_allclose(eta_e + eta_d, 1.0, atol=1e-6) + self.assertTrue(np.all(np.diff(eta_e) >= -1e-7)) + self.assertTrue(np.all(np.diff(eta_d) <= 1e-7)) + self.assertTrue(np.all((eta_e >= 0.5 - 1e-7) & (eta_e <= 1.0 + 1e-7))) + self.assertTrue(np.all((eta_d >= -1e-7) & (eta_d <= 0.5 + 1e-7))) + + +class LengthscaleFixture(unittest.TestCase): + """A stripe geometry whose physical lengthscale can be measured directly.""" + + def setUp(self): + self.lx = 0.5 + self.ly = 0.5 + self.resolution = 100 + self.target_length = 0.1 + self.filter_radius = self.target_length + self.beta = np.inf + self.eta_i = 0.5 + self.nx = int(np.round(self.lx * self.resolution)) + 1 + self.ny = int(np.round(self.ly * self.resolution)) + 1 + + coords = np.linspace(-self.lx / 2, self.lx / 2, self.nx) + self.x = coords[:, None] * np.ones((1, self.ny)) + + def _filter_and_project(self, rho, projection=ssp2): + rho_filtered = conic_filter( + jnp.asarray(rho), self.filter_radius, self.lx, self.ly, self.resolution + ) + rho_projected = projection( + rho_filtered, self.beta, self.eta_i, self.resolution + ) + return rho_filtered, rho_projected + + def _feature_width(self, rho_projected, solid=True): + """Width of the central feature of the projected design, in physical units.""" + profile = np.asarray(rho_projected)[:, self.ny // 2] + mask = profile > 0.5 if solid else profile < 0.5 + return float(np.count_nonzero(mask)) / self.resolution + + +class TestSolidConstraint(LengthscaleFixture): + def test_detects_thin_solid_features(self): + """A stripe thinner than the target violates the solid constraint.""" + for latent_width in (0.06, 0.08): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) <= latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected) + self.assertGreater(width, 0.0) + self.assertLess(width, self.target_length) + + value = float( + constraint_solid( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertGreater(value, 0.0) + + def test_accepts_thick_solid_features(self): + """A stripe at or above the target satisfies the solid constraint.""" + for latent_width in (0.10, 0.14, 0.20): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) <= latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected) + self.assertGreaterEqual(width, self.target_length) + + value = float( + constraint_solid( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertLessEqual(value, 0.0) + + def test_uniform_solid_design_is_feasible(self): + """A fully solid design has no interfaces and no lengthscale violation.""" + rho_filtered = jnp.ones((self.nx, self.ny)) + rho_projected = jnp.ones((self.nx, self.ny)) + + value = float( + constraint_solid( + rho_filtered, rho_projected, self.resolution, self.target_length + ) + ) + self.assertAlmostEqual(value, -1.0, places=6) + + +class TestVoidConstraint(LengthscaleFixture): + def test_detects_thin_void_features(self): + """A gap thinner than the target violates the void constraint.""" + for latent_width in (0.06, 0.08): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) > latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected, solid=False) + self.assertGreater(width, 0.0) + self.assertLess(width, self.target_length) + + value = float( + constraint_void( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertGreater(value, 0.0) + + def test_accepts_thick_void_features(self): + """A gap at or above the target satisfies the void constraint.""" + for latent_width in (0.10, 0.14, 0.20): + with self.subTest(latent_width=latent_width): + rho = (np.abs(self.x) > latent_width / 2).astype(float) + rho_filtered, rho_projected = self._filter_and_project(rho) + + width = self._feature_width(rho_projected, solid=False) + self.assertGreaterEqual(width, self.target_length) + + value = float( + constraint_void( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + ) + self.assertLessEqual(value, 0.0) + + def test_uniform_void_design_is_feasible(self): + rho_filtered = jnp.zeros((self.nx, self.ny)) + rho_projected = jnp.zeros((self.nx, self.ny)) + + value = float( + constraint_void( + rho_filtered, rho_projected, self.resolution, self.target_length + ) + ) + self.assertAlmostEqual(value, -1.0, places=6) + + +class TestAnyProjectionOrder(LengthscaleFixture): + """The constraints only see `rho_projected`, so any SSP order works.""" + + def _design(self): + rng = np.random.default_rng(42) + return rng.random((self.nx, self.ny)) + + def test_gradient_flows_through_each_projection(self): + projections = { + "tanh_projection": tanh_projection, + "ssp1_bilinear": ssp1_bilinear, + "ssp2": ssp2, + } + rho = self._design() + + for name, projection in projections.items(): + for constraint in (constraint_solid, constraint_void): + with self.subTest(projection=name, constraint=constraint.__name__): + + def objective(rho_flat, projection=projection, constraint=constraint): + rho_design = rho_flat.reshape((self.nx, self.ny)) + rho_filtered = conic_filter( + rho_design, + self.filter_radius, + self.lx, + self.ly, + self.resolution, + ) + if projection is tanh_projection: + # A finite beta keeps the plain tanh projection smooth. + rho_projected = projection(rho_filtered, 8.0, self.eta_i) + else: + rho_projected = projection( + rho_filtered, self.beta, self.eta_i, self.resolution + ) + return constraint( + rho_filtered, + rho_projected, + self.resolution, + self.target_length, + ) + + value, gradient = value_and_grad(objective)( + jnp.asarray(rho.ravel()) + ) + + self.assertTrue(np.isfinite(float(value))) + gradient = np.asarray(gradient) + self.assertTrue(np.isfinite(gradient).all()) + self.assertGreater(np.linalg.norm(gradient), 0.0) + + def test_jit_matches_eager(self): + rho = self._design() + rho_filtered, rho_projected = self._filter_and_project(rho) + + def constraints(rho_filtered, rho_projected): + return ( + constraint_solid( + rho_filtered, rho_projected, self.resolution, self.target_length + ), + constraint_void( + rho_filtered, rho_projected, self.resolution, self.target_length + ), + ) + + eager = [float(v) for v in constraints(rho_filtered, rho_projected)] + compiled = [ + float(v) for v in jit(constraints)(rho_filtered, rho_projected) + ] + + np.testing.assert_allclose(compiled, eager, rtol=1e-5) + + +class TestGradientsAgainstFiniteDifferences(unittest.TestCase): + """Compare reverse-mode gradients with central finite differences. + + Double precision is required here: the normalized constraint is O(1/epsilon), + so differencing it in single precision is dominated by roundoff. + """ + + def test_adjoints_match_finite_differences(self): + with enable_x64(): + lx = ly = 0.4 + resolution = 50 + target_length = 0.15 + nx = int(np.round(lx * resolution)) + 1 + ny = int(np.round(ly * resolution)) + 1 + + rng = np.random.default_rng(0) + rho = rng.random((nx, ny)) + rho_filtered = conic_filter( + jnp.asarray(rho), target_length, lx, ly, resolution + ) + rho_projected = ssp2(rho_filtered, np.inf, 0.5, resolution) + perturbation = rng.standard_normal((nx, ny)) + step = 1e-6 + + for constraint in (constraint_solid, constraint_void): + for argument in ("rho_filtered", "rho_projected"): + with self.subTest( + constraint=constraint.__name__, argument=argument + ): + + def scalar(value, argument=argument, constraint=constraint): + if argument == "rho_filtered": + return constraint( + value, rho_projected, resolution, target_length + ) + return constraint( + rho_filtered, value, resolution, target_length + ) + + base = ( + rho_filtered + if argument == "rho_filtered" + else rho_projected + ) + adjoint = float( + np.sum(np.asarray(grad(scalar)(base)) * perturbation) + ) + finite_difference = ( + float(scalar(base + step * perturbation)) + - float(scalar(base - step * perturbation)) + ) / (2 * step) + + self.assertAlmostEqual( + adjoint / finite_difference, 1.0, places=5 + ) + + +class TestArgumentValidation(unittest.TestCase): + def setUp(self): + self.rho_filtered = jnp.full((8, 8), 0.5) + self.rho_projected = jnp.zeros((8, 8)) + + def test_negative_target_length(self): + with self.assertRaises(ValueError): + constraint_solid(self.rho_filtered, self.rho_projected, 10, -0.1) + + def test_nonpositive_conic_radius(self): + with self.assertRaises(ValueError): + constraint_void( + self.rho_filtered, self.rho_projected, 10, 0.1, conic_radius=0.0 + ) + + def test_mismatched_shapes(self): + with self.assertRaises(ValueError): + constraint_solid(self.rho_filtered, jnp.zeros((8, 4)), 10, 0.1) + + def test_non_2d_input(self): + with self.assertRaises(ValueError): + constraint_solid(jnp.full((4, 4, 4), 0.5), jnp.zeros((4, 4, 4)), 10, 0.1) + + +if __name__ == "__main__": + unittest.main() From 75edbcb26c1e015020807d2f6eecfa223c5600a0 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 18 Sep 2026 14:33:22 -0700 Subject: [PATCH 2/3] Address review comments on the lengthscale constraints MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cite the published version of Arrieta et al. (Struct. Multidisc. Optim. 69, 210, 2026) in place of the arXiv preprint, use the paper's ρ̃ / ρ̄ notation for the filtered and projected densities in the documentation, and report the figure of merit and constraints recomputed at the design nlopt returns rather than at the last design it happened to evaluate. Co-Authored-By: Claude Opus 5 --- examples/python/ssp_constrained_example.py | 17 ++++-- src/python/ssp_topopt/constraints.py | 62 +++++++++++----------- tests/python/test_constraints.py | 9 ++-- 3 files changed, 48 insertions(+), 40 deletions(-) diff --git a/examples/python/ssp_constrained_example.py b/examples/python/ssp_constrained_example.py index 075a82a..b49b1df 100644 --- a/examples/python/ssp_constrained_example.py +++ b/examples/python/ssp_constrained_example.py @@ -15,8 +15,8 @@ constraints on and repairs the design, at a small cost in the figure of merit. Ref: R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free -minimum-lengthscale constraints for topology optimization," arXiv.org e-Print -archive, 2507.16108, July 2025. +minimum-lengthscale constraints for topology optimization," Structural and +Multidisciplinary Optimization, vol. 69, p. 210, September 2026. """ import time @@ -199,12 +199,19 @@ def nlopt_constraint(x, grad_out, fn=value_and_grad_fn): x_opt = opt.optimize(np.asarray(x_init, dtype=float)) elapsed = time.perf_counter() - start + # The last entry of the histories is the last design nlopt *evaluated*, + # which is not necessarily the one it returns, so recompute at x_opt. + x_opt_jnp = jnp.asarray(x_opt) + fom_opt = float(objective_and_grad(x_opt_jnp)[0]) + solid_opt = float(constraints_and_grad["solid"](x_opt_jnp)[0]) + void_opt = float(constraints_and_grad["void"](x_opt_jnp)[0]) + label = "constrained" if constrained else "unconstrained" print( f" {label}: {maxeval} evaluations in {elapsed:6.1f}s " - f"FOM={fom_history[-1]:.4e} " - f"solid={constraint_history['solid'][-1]:+.3e} " - f"void={constraint_history['void'][-1]:+.3e}" + f"FOM={fom_opt:.4e} " + f"solid={solid_opt:+.3e} " + f"void={void_opt:+.3e}" ) return x_opt diff --git a/src/python/ssp_topopt/constraints.py b/src/python/ssp_topopt/constraints.py index 8a98757..38764a7 100644 --- a/src/python/ssp_topopt/constraints.py +++ b/src/python/ssp_topopt/constraints.py @@ -1,35 +1,35 @@ """Minimum-lengthscale (geometric) constraints for topology optimization. These are the geometric constraints of Zhou et al. (2015) combined with the -hyperparameter-free thresholds derived by Arrieta et al. (2025). The constraints +hyperparameter-free thresholds derived by Arrieta et al. (2026). The constraints penalize solid (or void) features whose lengthscale falls below a target value, -and are formulated purely in terms of the filtered density and the projected -density. Consequently they work with *any* order of subpixel smoothing -- -`ssp1_bilinear`, `ssp2`, or even a plain `tanh_projection` -- since the -projection only enters through `rho_projected`. +and are formulated purely in terms of the filtered density ρ̃ (`rho_filtered`) +and the projected density ρ̄ (`rho_projected`). Consequently they work with +*any* order of subpixel smoothing -- `ssp1_bilinear`, `ssp2`, or even a plain +`tanh_projection` -- since the projection only enters through ρ̄. The solid constraint reads - g_s = (1/N) Σ_i I_s,i [min(rho_filtered_i - eta_e, 0)]^2 , - I_s = rho_projected * exp(-c |∇ rho_filtered|^2) , + g_s = (1/N) Σ_i I_s,i [min(ρ̃_i - η_e, 0)]^2 , + I_s = ρ̄ * exp(-c |∇ρ̃|^2) , and the void constraint is the complementary expression - g_v = (1/N) Σ_i I_v,i [min(eta_d - rho_filtered_i, 0)]^2 , - I_v = (1 - rho_projected) * exp(-c |∇ rho_filtered|^2) , + g_v = (1/N) Σ_i I_v,i [min(η_d - ρ̃_i, 0)]^2 , + I_v = (1 - ρ̄) * exp(-c |∇ρ̃|^2) , where the structural functions I_s and I_v single out the "inflection regions" -of the design (the interior of a feature, where the filtered density is -stationary), and eta_e, eta_d are the eroded/dilated threshold points of the -conic filter. Following Arrieta et al., the decay rate is c = 64 R^2 and the -constraint threshold is eps = 1e-8, where R is the conic filter radius; the -constraint is well behaved for target_length / R roughly in [0.25, 1.5]. +of the design (the interior of a feature, where ρ̃ is stationary), and η_e, η_d +are the eroded/dilated threshold points of the conic filter. Following Arrieta +et al., the decay rate is c = 64 R^2 and the constraint threshold is eps = 1e-8, +where R is the conic filter radius; the constraint is well behaved for +target_length / R roughly in [0.25, 1.5]. Refs: R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free minimum-lengthscale -constraints for topology optimization," arXiv.org e-Print archive, 2507.16108, -July 2025. +constraints for topology optimization," Structural and Multidisciplinary +Optimization, vol. 69, p. 210, September 2026. M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, "Minimum length scale in topology optimization by geometric constraints," Computer Methods in Applied Mechanics and @@ -46,7 +46,7 @@ from .utils import ArrayLikeType, gradient -# Hyperparameters derived in section 4.1 of Arrieta et al. (2025). The decay rate +# Hyperparameters derived in section 4.1 of Arrieta et al. (2026). The decay rate # is expressed relative to the square of the conic filter radius, i.e. the actual # decay rate is c = DEFAULT_CONSTRAINT_DECAYRATE * conic_radius**2. DEFAULT_CONSTRAINT_THRESHOLD = 1e-8 @@ -54,9 +54,9 @@ def solid_threshold(lengthscale_ratio: float): - """The eroded threshold point eta_e of a conic filter. + """The eroded threshold point η_e of a conic filter. - Ref: Eq. (9) of Arrieta et al. (2025), originally from Qian and Sigmund (2013). + Ref: Eq. (9) of Arrieta et al. (2026), originally from Qian and Sigmund (2013). Args: lengthscale_ratio: the ratio of the target lengthscale to the conic @@ -74,9 +74,9 @@ def solid_threshold(lengthscale_ratio: float): def void_threshold(lengthscale_ratio: float): - """The dilated threshold point eta_d of a conic filter. + """The dilated threshold point η_d of a conic filter. - Ref: Eq. (12) of Arrieta et al. (2025), originally from Qian and Sigmund (2013). + Ref: Eq. (12) of Arrieta et al. (2026), originally from Qian and Sigmund (2013). Args: lengthscale_ratio: the ratio of the target lengthscale to the conic @@ -160,16 +160,16 @@ def constraint_solid( ): """Calculate a solid minimum-lengthscale constraint function. - This technique takes smoothed data, e.g. from filtering, `rho_filtered` and - binary data, e.g. from projection, `rho_projected` both defined on the same - grid, and measures whether features in the solid region, i.e. where - `rho_projected` takes values of 1, violate the minimum `target_length`. + This technique takes smoothed data ρ̃, e.g. from filtering, `rho_filtered` + and binary data ρ̄, e.g. from projection, `rho_projected` both defined on the + same grid, and measures whether features in the solid region, i.e. where ρ̄ + takes values of 1, violate the minimum `target_length`. The returned value is normalized by the constraint threshold, so the constraint is satisfied when the value is nonpositive and violated when it is positive. It can therefore be handed directly to a nonlinear optimizer such as `nlopt` as an inequality constraint. The unnormalized constraint value of - Eq. (7) of Arrieta et al. (2025) is `(value + 1) * constraint_threshold`. + Eq. (7) of Arrieta et al. (2026) is `(value + 1) * constraint_threshold`. Any projection may be used to produce `rho_projected`, including `ssp1_bilinear`, `ssp2`, and `tanh_projection`. @@ -221,16 +221,16 @@ def constraint_void( ): """Calculate a void minimum-lengthscale constraint function. - This technique takes smoothed data, e.g. from filtering, `rho_filtered` and - binary data, e.g. from projection, `rho_projected` both defined on the same - grid, and measures whether features in the void region, i.e. where - `rho_projected` takes values of 0, violate the minimum `target_length`. + This technique takes smoothed data ρ̃, e.g. from filtering, `rho_filtered` + and binary data ρ̄, e.g. from projection, `rho_projected` both defined on the + same grid, and measures whether features in the void region, i.e. where ρ̄ + takes values of 0, violate the minimum `target_length`. The returned value is normalized by the constraint threshold, so the constraint is satisfied when the value is nonpositive and violated when it is positive. It can therefore be handed directly to a nonlinear optimizer such as `nlopt` as an inequality constraint. The unnormalized constraint value of - Eq. (10) of Arrieta et al. (2025) is `(value + 1) * constraint_threshold`. + Eq. (10) of Arrieta et al. (2026) is `(value + 1) * constraint_threshold`. Any projection may be used to produce `rho_projected`, including `ssp1_bilinear`, `ssp2`, and `tanh_projection`. diff --git a/tests/python/test_constraints.py b/tests/python/test_constraints.py index db5cc6d..b1fbd0c 100644 --- a/tests/python/test_constraints.py +++ b/tests/python/test_constraints.py @@ -2,9 +2,10 @@ The reference behavior is the one described in R. Arrieta, G. Romano, and S. G. Johnson, "Hyperparameter-free minimum-lengthscale constraints for topology -optimization," arXiv:2507.16108 (2025): with the derived hyperparameters, the -constraint is violated (positive) when the physical lengthscale of a feature -falls below the target lengthscale and satisfied (nonpositive) otherwise. +optimization," Struct. Multidisc. Optim. 69, 210 (2026): with the derived +hyperparameters, the constraint is violated (positive) when the physical +lengthscale of a feature falls below the target lengthscale and satisfied +(nonpositive) otherwise. """ import unittest @@ -29,7 +30,7 @@ class TestThresholdFunctions(unittest.TestCase): """The conic-filter threshold points of Qian and Sigmund (2013).""" def test_known_values(self): - # Eqs. (9) and (12) of Arrieta et al. (2025). + # Eqs. (9) and (12) of Arrieta et al. (2026). self.assertAlmostEqual(float(solid_threshold(0.0)), 0.5) self.assertAlmostEqual(float(solid_threshold(1.0)), 0.75) self.assertAlmostEqual(float(solid_threshold(2.0)), 1.0) From 7ed53ab363bd73af504e2d871abc8ddaff1de05a Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 18 Sep 2026 14:46:58 -0700 Subject: [PATCH 3/3] Cite the published lengthscale-constraints paper in the README The top-level reference list still pointed at the arXiv preprint for the paper the constraints are based on, which appeared as Struct. Multidisc. Optim. 69, 210 (September 2026). Link the journal DOI instead, matching the style of the Optics Express entry above it. Co-Authored-By: Claude Opus 5 --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f5dddb1..63c671b 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ This is a repository for code implementing the **smoothed subpixel projection (S * A. M. Hammond, A. Oskooi, I. M. Hammond, M. Chen, S. E. Ralph, and S. G. Johnson, [“Unifying and accelerating level-set and density-based topology optimization by subpixel-smoothed projection,”](http://doi.org/10.1364/OE.563512) Optics Express, vol. 33, pp. 33620–33642, July 2025. Editor's Pick. * G. Romano, R. Arrieta, and S. G. Johnson, [“Differentiating through binarized topology changes: Second-order subpixel-smoothed projection,”](http://arxiv.org/abs/2601.10737) arXiv.org e-Print archive, 2601.10737, January 2026. -* R. Arrieta, G. Romano, and S. G. Johnson, [“Hyperparameter-free minimum-lengthscale constraints for topology optimization,”](http://arxiv.org/abs/2507.16108) arXiv.org e-Print archive, 2507.16108, July 2025. +* R. Arrieta, G. Romano, and S. G. Johnson, [“Hyperparameter-free minimum-lengthscale constraints for topology optimization,”](http://doi.org/10.1007/s00158-026-04388-6) Structural and Multidisciplinary Optimization, vol. 69, p. 210, September 2026. ## Documentation