From 355b5986c7a8b04b5fca44ea1f776b4c0fe74787 Mon Sep 17 00:00:00 2001 From: Magnus Fjeld Olsen <77971426+magnusfjeldolsen@users.noreply.github.com> Date: Fri, 29 May 2026 09:11:23 +0200 Subject: [PATCH 1/2] perf: avoid np.isclose in preprocess_strains_with_limits preprocess_strains_with_limits runs at the top of every constitutive law's get_stress/get_tangent, so it is called for every strain evaluation in the section integrators. It called np.isclose twice and allocated two zeros_like temporaries on every call, even though it only compares each strain against the constant ultimate-strain limit. Replace this with precomputed scalar tolerances and plain abs(...) <= tol comparisons that exactly reproduce np.isclose(atol=1e-6) with the default rtol=1e-5, so results are bit-identical. The fiber integrator benefits most, since every fiber goes through get_stress: ~1.8x faster calculations across randomly generated RC sections (median over 40 sections). Marin gains little, as it does far fewer strain evaluations. --- structuralcodes/core/base/_constitutive_law.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/structuralcodes/core/base/_constitutive_law.py b/structuralcodes/core/base/_constitutive_law.py index badca37e..74d57efc 100644 --- a/structuralcodes/core/base/_constitutive_law.py +++ b/structuralcodes/core/base/_constitutive_law.py @@ -61,16 +61,21 @@ def preprocess_strains_with_limits( eps = eps if np.isscalar(eps) else np.atleast_1d(eps) eps_min, eps_max = self.get_ultimate_strain() + # Absolute thresholds reproducing np.isclose(atol=1e-6) with the + # default rtol=1e-5, i.e. atol + rtol * abs(limit). Precomputing + # these scalars avoids the per-call np.isclose machinery and the + # zeros_like allocations, which are costly in the integration loop. + tol_max = 1e-6 + 1e-5 * abs(eps_max) + tol_min = 1e-6 + 1e-5 * abs(eps_min) + if np.isscalar(eps): - if np.isclose(eps, eps_max, atol=1e-6): + if abs(eps - eps_max) <= tol_max: return eps_max - if np.isclose(eps, eps_min, atol=1e-6): + if abs(eps - eps_min) <= tol_min: return eps_min return eps - idxs = np.isclose(eps, np.zeros_like(eps) + eps_max, atol=1e-6) - eps[idxs] = eps_max - idxs = np.isclose(eps, np.zeros_like(eps) + eps_min, atol=1e-6) - eps[idxs] = eps_min + eps[np.abs(eps - eps_max) <= tol_max] = eps_max + eps[np.abs(eps - eps_min) <= tol_min] = eps_min return eps From bbc0b656937c0f236fd431c77b96a0da29f41d5b Mon Sep 17 00:00:00 2001 From: Magnus Fjeld Olsen <77971426+magnusfjeldolsen@users.noreply.github.com> Date: Fri, 5 Jun 2026 12:43:40 +0200 Subject: [PATCH 2/2] test: cover preprocess_strains_with_limits across all constitutive laws Assert that the precomputed-threshold implementation reproduces the documented np.isclose(atol=1e-6) tolerance bit-for-bit, and that strains within the band snap to the exact limit while strains outside it are returned unchanged. Parametrized over all nine constitutive laws. --- .../test_materials/test_constitutive_laws.py | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/tests/test_materials/test_constitutive_laws.py b/tests/test_materials/test_constitutive_laws.py index c339a0e9..f92e4723 100644 --- a/tests/test_materials/test_constitutive_laws.py +++ b/tests/test_materials/test_constitutive_laws.py @@ -742,3 +742,92 @@ def test_parallel_get_ultimate_strain(E1, fy1, eps_su1, E2, fy2, eps_su2): expected_ult_strain = max(eps_su1, eps_su2) assert math.isclose(ult_strain[1], expected_ult_strain) + + +def _snap_with_isclose(law, eps): + """Reference snapping using the original ``np.isclose`` formulation. + + Documents the tolerance that ``preprocess_strains_with_limits`` must + reproduce: ``np.isclose`` with ``atol=1e-6`` and the default + ``rtol=1e-5``. + """ + eps = eps if np.isscalar(eps) else np.atleast_1d(eps) + eps_min, eps_max = law.get_ultimate_strain() + if np.isscalar(eps): + if np.isclose(eps, eps_max, atol=1e-6): + return eps_max + if np.isclose(eps, eps_min, atol=1e-6): + return eps_min + return eps + idxs = np.isclose(eps, np.zeros_like(eps) + eps_max, atol=1e-6) + eps[idxs] = eps_max + idxs = np.isclose(eps, np.zeros_like(eps) + eps_min, atol=1e-6) + eps[idxs] = eps_min + return eps + + +_PREPROCESS_LAWS = [ + Elastic(E=30000), + ElasticPlastic(E=200000, fy=500, eps_su=0.07), + ParabolaRectangle(fc=35), + ParabolaRectangle(fc=35, n=2.3), + BilinearCompression(fc=35, eps_c=-0.002, eps_cu=-0.0035), + Sargin(fc=38), + Popovics(fc=35), + UserDefined(x=[-0.0035, -0.002, 0.0, 0.001], y=[0.0, -35.0, 0.0, 3.0]), + InitialStrain( + constitutive_law=ElasticPlastic(E=200000, fy=500, eps_su=0.07), + initial_strain=1e-3, + ), + Parallel( + [ + ElasticPlastic(E=200000, fy=450, eps_su=0.06), + ElasticPlastic(E=150000, fy=800, eps_su=0.01), + ] + ), +] + + +@pytest.mark.parametrize('law', _PREPROCESS_LAWS) +def test_preprocess_strains_matches_isclose_tolerance(law): + """preprocess_strains_with_limits reproduces the np.isclose tolerance. + + The optimized scalar-threshold implementation must give results that are + bit-identical to the documented ``np.isclose(atol=1e-6)`` behaviour for + every constitutive law. + """ + eps_min, eps_max = law.get_ultimate_strain() + rng = np.random.default_rng(0) + eps = np.concatenate( + [ + rng.uniform(eps_min, eps_max, 5000), + np.linspace(eps_min - 3e-6, eps_min + 3e-6, 2000), + np.linspace(eps_max - 3e-6, eps_max + 3e-6, 2000), + np.array([eps_min, eps_max, 0.0]), + ] + ) + expected = _snap_with_isclose(law, eps.copy()) + result = law.preprocess_strains_with_limits(eps.copy()) + np.testing.assert_array_equal(result, expected) + + +@pytest.mark.parametrize('law', _PREPROCESS_LAWS) +def test_preprocess_strains_snaps_within_band(law): + """Strains within the tolerance band snap to the exact limit and + strains outside it are returned unchanged. + """ + eps_min, eps_max = law.get_ultimate_strain() + tol_max = 1e-6 + 1e-5 * abs(eps_max) + tol_min = 1e-6 + 1e-5 * abs(eps_min) + + # Just inside the band -> snapped to the exact limit. + assert law.preprocess_strains_with_limits(eps_max - 0.4 * tol_max) == ( + eps_max + ) + assert law.preprocess_strains_with_limits(eps_min + 0.4 * tol_min) == ( + eps_min + ) + + # Clearly outside the band -> returned unchanged. + inside = eps_max - 10 * tol_max + assert law.preprocess_strains_with_limits(inside) == inside