From e2c56e41b4fb0518df97f6264bd97cb34801e233 Mon Sep 17 00:00:00 2001
From: Morten Engen <morten.engen@multiconsult.no>
Date: Mon, 29 Sep 2025 22:10:53 +0200
Subject: [PATCH 1/8] Implement custom warning classes

---
 structuralcodes/core/errors.py | 11 +++++++++++
 1 file changed, 11 insertions(+)
 create mode 100644 structuralcodes/core/errors.py

diff --git a/structuralcodes/core/errors.py b/structuralcodes/core/errors.py
new file mode 100644
index 00000000..059dd6af
--- /dev/null
+++ b/structuralcodes/core/errors.py
@@ -0,0 +1,11 @@
+"""General exception and warning classes."""
+
+
+class StructuralCodesWarning(Warning):
+    """Base class for StructurlCodes warnings."""
+
+
+class NoConvergenceWarning(StructuralCodesWarning):
+    """A warning that indicates that no convergence was reached for an
+    iterative solver.
+    """

From 9db64167fbf4fef4acafc5a462678de565cec1da Mon Sep 17 00:00:00 2001
From: Morten Engen <morten.engen@multiconsult.no>
Date: Mon, 29 Sep 2025 22:11:42 +0200
Subject: [PATCH 2/8] Use NoConvergenceWarning to indicate no convergence in
 iterative solver

---
 structuralcodes/sections/_generic.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index 75a0440f..dfaee213 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -14,6 +14,7 @@
 
 import structuralcodes.core._section_results as s_res
 from structuralcodes.core.base import Section, SectionCalculator
+from structuralcodes.core.errors import NoConvergenceWarning
 from structuralcodes.geometry import (
     CompoundGeometry,
     PointGeometry,
@@ -1474,6 +1475,12 @@ def calculate_strain_profile(
                 break
 
         if num_iter >= max_iter:
-            raise StopIteration('Maximum number of iterations reached.')
+            warnings.warn(
+                message=(
+                    'Maximum number of iterations reached. '
+                    'Please review the strain profile carefully.'
+                ),
+                category=NoConvergenceWarning,
+            )
 
         return strain.tolist()

From 6bcb522757889f8a10949c7408ff31231f6a89d2 Mon Sep 17 00:00:00 2001
From: Morten Engen <morten.engen@multiconsult.no>
Date: Mon, 29 Sep 2025 22:12:05 +0200
Subject: [PATCH 3/8] Add filter to make sure custom warnings are always shown

---
 structuralcodes/__init__.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/structuralcodes/__init__.py b/structuralcodes/__init__.py
index 46182d71..aeb8a49e 100644
--- a/structuralcodes/__init__.py
+++ b/structuralcodes/__init__.py
@@ -1,7 +1,10 @@
 """A Python package that contains models from structural design codes."""
 
+import warnings
+
 from . import codes, core, geometry, materials, sections
 from .codes import get_design_codes, set_design_code, set_national_annex
+from .core.errors import StructuralCodesWarning
 
 __version__ = '0.6.0'
 
@@ -15,3 +18,5 @@
     'geometry',
     'sections',
 ]
+
+warnings.filterwarnings(action='always', category=StructuralCodesWarning)

From 70dad1e8ce901e01cefc37367de55ddb16d4c0ec Mon Sep 17 00:00:00 2001
From: Morten Engen <morten.engen@multiconsult.no>
Date: Mon, 29 Sep 2025 22:12:46 +0200
Subject: [PATCH 4/8] Update test

---
 tests/test_sections/test_generic_section.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py
index b811fde0..83ab4d3a 100644
--- a/tests/test_sections/test_generic_section.py
+++ b/tests/test_sections/test_generic_section.py
@@ -7,6 +7,7 @@
 from shapely import Polygon
 
 from structuralcodes.codes.ec2_2004 import reinforcement_duct_props
+from structuralcodes.core.errors import NoConvergenceWarning
 from structuralcodes.geometry import (
     CircularGeometry,
     RectangularGeometry,
@@ -1109,8 +1110,8 @@ def test_strain_plane_calculation_rectangular_rc_high_load(
 
     # Check that with given loads we don't reach convergence
     section = GenericSection(geo)
-    with pytest.raises(
-        StopIteration, match='Maximum number of iterations reached'
+    with pytest.warns(
+        NoConvergenceWarning, match='Maximum number of iterations reached'
     ):
         section.section_calculator.calculate_strain_profile(
             n, my, mz, tol=1e-7

From 73e2e47539914989bfd0f90a587f2102e76f5a27 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Fri, 17 Oct 2025 14:52:05 +0200
Subject: [PATCH 5/8] max_iter as optional input arguments on section
 calculator

---
 structuralcodes/sections/_generic.py | 78 +++++++++++++++++++++-------
 1 file changed, 58 insertions(+), 20 deletions(-)

diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index dfaee213..e35c5f4a 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -370,9 +370,14 @@ def get_balanced_failure_strain(
         return (y_n, y_p, strain)
 
     def find_equilibrium_fixed_pivot(
-        self, geom: CompoundGeometry, n: float, yielding: bool = False
+        self,
+        geom: CompoundGeometry,
+        n: float,
+        yielding: bool = False,
+        max_iter: int = 100,
+        tol: float = 1e-2,
     ) -> t.List[float]:
-        """Find the equilibrium changing curvature fixed a pivot.
+        r"""Find the equilibrium changing curvature fixed a pivot.
         The algorithm uses bisection algorithm between curvature
         of balanced failure and 0. Selected the pivot point as
         the top or the bottom one, the neutral axis is lowered or
@@ -382,15 +387,19 @@ def find_equilibrium_fixed_pivot(
             geom (CompoundGeometry): A geometry in the rotated reference
                 system.
             n (float): Value of external axial force needed to be equilibrated.
-            yielding (bool): ...
+            yielding (bool): If true, the yielding strain is used as the
+                ultimate one, therefore finding yielding strength and not
+                ultimate strength
+            max_iter (int): the maximum number of iterations in the iterative
+                process (default = 100).
+            tol (float): the tolerance for convergence test in terms of
+                $\deltaN_a - \deltaN_b$ (default = 1e-2).
 
         Returns:
             List(float): 3 floats: Axial strain at (0,0), and curvatures of y*
             and z* axes. Note that being uniaxial bending,
             curvature along z* is 0.0.
         """
-        # Number of maximum iteration for the bisection algorithm
-        ITMAX = 100
         # 1. Start with a balanced failure: this is found from all ultimate
         # strains for all materials, checking the minimum curvature value
         y_n, y_p, strain = self.get_balanced_failure_strain(geom, yielding)
@@ -431,7 +440,7 @@ def find_equilibrium_fixed_pivot(
         )
         dn_b = n_int - n
         it = 0
-        while (abs(dn_a - dn_b) > 1e-2) and (it < ITMAX):
+        while (abs(dn_a - dn_b) > tol) and (it < max_iter):
             chi_c = (chi_a + chi_b) / 2.0
             eps_0 = strain_pivot - chi_c * pivot
             (
@@ -450,9 +459,14 @@ def find_equilibrium_fixed_pivot(
                 chi_a = chi_c
                 dn_a = dn_c
             it += 1
-        if it >= ITMAX:
-            s = f'Last iteration reached a unbalance of {dn_c}'
-            raise ValueError(f'Maximum number of iterations reached.\n{s}')
+        if it >= max_iter:
+            msg = 'GenericSectionCalculator::find_equilibrium_fixed_pivot\n\t'
+            msg += 'Maximum number of iterations reached.\n\t'
+            msg += f'Last iteration reached a unbalance of {dn_c}'
+            warnings.warn(
+                message=msg,
+                category=NoConvergenceWarning,
+            )
         # Found equilibrium
         # Return the strain distribution
         return [eps_0, chi_c, 0]
@@ -464,6 +478,8 @@ def _prefind_range_curvature_equilibrium(
         curv: float,
         eps_0_a: float,
         dn_a: float,
+        max_iter: int = 20,
+        max_restart_attemps: int = 20,
     ):
         """Perfind range where the curvature equilibrium is located.
 
@@ -471,8 +487,6 @@ def _prefind_range_curvature_equilibrium(
         existence of at least one zero in the function dn vs. curv in order to
         apply the bisection algorithm.
         """
-        ITMAX = 20
-        MAXRESTATTEMPTS = 20
         sign = -1 if dn_a > 0 else 1
         found = False
         it = 0
@@ -482,7 +496,7 @@ def _prefind_range_curvature_equilibrium(
         r = 2.0
         diverging = False
         diverging_steps = 0
-        while not found and it < ITMAX and restarts < MAXRESTATTEMPTS:
+        while not found and it < max_iter and restarts < max_restart_attemps:
             eps_0_b = eps_0_a + sign * delta * r ** (it)
             (
                 n_int,
@@ -511,16 +525,24 @@ def _prefind_range_curvature_equilibrium(
                     diverging = False
                     diverging_steps = 0
             it += 1
-        if it >= ITMAX and not found:
+        if it >= max_iter and not found:
             s = f'Last iteration reached a unbalance of: \
                 dn_a = {dn_a} dn_b = {dn_b})'
+            # This should be kept as an exception otherwise afterwards
+            # the bisection will fail
             raise ValueError(f'Maximum number of iterations reached.\n{s}')
         return (eps_0_b, dn_b)
 
     def find_equilibrium_fixed_curvature(
-        self, geom: CompoundGeometry, n: float, curv: float, eps_0: float
+        self,
+        geom: CompoundGeometry,
+        n: float,
+        curv: float,
+        eps_0: float,
+        max_iter: int = 100,
+        tol: float = 1e-2,
     ) -> t.Tuple[float, float, float]:
-        """Find strain profile with equilibrium with fixed curvature.
+        r"""Find strain profile with equilibrium with fixed curvature.
 
         Given curvature and external axial force, find the strain profile that
         makes internal and external axial force in equilibrium.
@@ -530,14 +552,16 @@ def find_equilibrium_fixed_curvature(
             n (float): The external axial load.
             curv (float): The value of curvature.
             eps_0 (float): A first attempt for neutral axis position.
+            max_iter (int): the maximum number of iterations in the iterative
+                process (default = 100).
+            tol (float): the tolerance for convergence test in terms of
+                $\deltaN_a - \deltaN_b$. (default = 1e-2)
 
         Returns:
             Tuple(float, float, float): The axial strain and the two
             curvatures.
         """
         # Useful for Moment Curvature Analysis
-        # Number of maximum iteration for the bisection algorithm
-        ITMAX = 100
         # Start from previous position of N.A.
         eps_0_a = eps_0
         # find internal axial force by integration
@@ -561,7 +585,7 @@ def find_equilibrium_fixed_curvature(
         )
         # Found a range within there is the solution, apply bisection
         it = 0
-        while (abs(dn_a - dn_b) > 1e-2) and (it < ITMAX):
+        while (abs(dn_a - dn_b) > tol) and (it < max_iter):
             eps_0_c = (eps_0_a + eps_0_b) / 2
             (
                 n_int,
@@ -579,9 +603,11 @@ def find_equilibrium_fixed_curvature(
                 dn_a = dn_c
                 eps_0_a = eps_0_c
             it += 1
-        if it >= ITMAX:
+        if it >= max_iter:
             s = f'Last iteration reached a unbalance of: \
                 dn_c = {dn_c}'
+            # This is used from moment-curvature anaylsis, so we call an
+            # exception otherwise the next part of the curve could be wrong
             raise ValueError(f'Maximum number of iterations reached.\n{s}')
         return eps_0_c, curv, 0
 
@@ -1371,7 +1397,19 @@ def calculate_mm_interaction_domain(
         res.strains = np.zeros((num_theta, 3))
         # Compute strength for given angle of NA
         for i, th in enumerate(res.theta):
-            res_bend_strength = self.calculate_bending_strength(theta=th, n=n)
+            # If for some reason there is the warning of non convergence,
+            # catch it and raise an Exception because the domain is not
+            # correctly computed
+            try:
+                with warnings.catch_warnings():
+                    warnings.simplefilter('error', NoConvergenceWarning)
+                    res_bend_strength = self.calculate_bending_strength(
+                        theta=th, n=n
+                    )
+            except NoConvergenceWarning as e:
+                raise RuntimeError(
+                    'Convergence cannot be found, the algorithm is stopped.'
+                ) from e
             # Save forces
             res.forces[i, 0] = n
             res.forces[i, 1] = res_bend_strength.m_y

From 788739489efb703b8f3bffa375999d036511ad4f Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Mon, 20 Oct 2025 12:29:30 +0200
Subject: [PATCH 6/8] add warnings and tests

---
 structuralcodes/sections/_generic.py        | 266 +++++++++++++-------
 tests/test_sections/test_generic_section.py |  88 +++++++
 2 files changed, 270 insertions(+), 84 deletions(-)

diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index e35c5f4a..24d4fc8d 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -461,8 +461,8 @@ def find_equilibrium_fixed_pivot(
             it += 1
         if it >= max_iter:
             msg = 'GenericSectionCalculator::find_equilibrium_fixed_pivot\n\t'
-            msg += 'Maximum number of iterations reached.\n\t'
-            msg += f'Last iteration reached a unbalance of {dn_c}'
+            msg += 'Maximum number of iterations reached.'
+            msg += f' Last iteration reached a unbalance of {dn_c:.3e}'
             warnings.warn(
                 message=msg,
                 category=NoConvergenceWarning,
@@ -526,11 +526,16 @@ def _prefind_range_curvature_equilibrium(
                     diverging_steps = 0
             it += 1
         if it >= max_iter and not found:
-            s = f'Last iteration reached a unbalance of: \
-                dn_a = {dn_a} dn_b = {dn_b})'
-            # This should be kept as an exception otherwise afterwards
-            # the bisection will fail
-            raise ValueError(f'Maximum number of iterations reached.\n{s}')
+            msg = (
+                'GenericSectionCalculator::'
+                '_prefind_range_curvature_equilibrium\n\t'
+                'Maximum number of iterations reached.'
+                f' Last iteration reached a unbalance of {dn_b:.3e}'
+            )
+            warnings.warn(
+                message=msg,
+                category=NoConvergenceWarning,
+            )
         return (eps_0_b, dn_b)
 
     def find_equilibrium_fixed_curvature(
@@ -604,11 +609,13 @@ def find_equilibrium_fixed_curvature(
                 eps_0_a = eps_0_c
             it += 1
         if it >= max_iter:
-            s = f'Last iteration reached a unbalance of: \
-                dn_c = {dn_c}'
-            # This is used from moment-curvature anaylsis, so we call an
-            # exception otherwise the next part of the curve could be wrong
-            raise ValueError(f'Maximum number of iterations reached.\n{s}')
+            msg = 'GenericSectionCalculator::find_equilibrium_fixed_curvature'
+            msg += '\n\tMaximum number of iterations reached. '
+            msg += f' Last iteration reached a unbalance of {dn_c:.3e}'
+            warnings.warn(
+                message=msg,
+                category=NoConvergenceWarning,
+            )
         return eps_0_c, curv, 0
 
     def calculate_limit_axial_load(self):
@@ -745,9 +752,9 @@ def integrate_strain_profile(
         return result[:-1]
 
     def calculate_bending_strength(
-        self, theta=0, n=0
+        self, theta=0, n=0, max_iter: int = 100, tol: float = 1e-2
     ) -> s_res.UltimateBendingMomentResults:
-        """Calculates the bending strength for given inclination of n.a. and
+        r"""Calculates the bending strength for given inclination of n.a. and
         axial load.
 
         Arguments:
@@ -755,6 +762,10 @@ def calculate_bending_strength(
                 radians, default = 0.
             n (float): Axial load applied to the section (+: tension, -:
                 compression), default = 0.
+            max_iter (int): the maximum number of iterations in the iterative
+                bisection process (default = 100).
+            tol (float): the tolerance for convergence test in terms of
+                $\deltaN_a - \deltaN_b$ (default = 1e-2
 
         Returns:
             UltimateBendingMomentResults: The results from the calculation.
@@ -771,7 +782,9 @@ def calculate_bending_strength(
 
         # Find the strain distribution corresponding to failure and equilibrium
         # with external axial force
-        strain = self.find_equilibrium_fixed_pivot(rotated_geom, n)
+        strain = self.find_equilibrium_fixed_pivot(
+            rotated_geom, n, max_iter=max_iter, tol=tol
+        )
         # Compute the internal forces with this strain distribution
         N, My, Mz, _ = self.integrator.integrate_strain_response_on_geometry(
             geo=rotated_geom, strain=strain, tri=self.triangulated_data
@@ -800,6 +813,92 @@ def calculate_bending_strength(
 
         return res
 
+    def _prepare_chi_array(
+        self,
+        rotated_geom,
+        n,
+        num_pre_yield,
+        num_post_yield,
+        chi_first,
+        max_iter,
+        tol,
+    ):
+        # Find ultimate curvature from the strain distribution
+        # corresponding to failure and equilibrium with external axial
+        # force
+        # This coul not converge, let's catch the warning to customize
+        # the message
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always', NoConvergenceWarning)
+            strain = self.find_equilibrium_fixed_pivot(
+                rotated_geom, n, max_iter=max_iter, tol=tol
+            )
+            chi_ultimate = strain[1]
+        if w:
+            # I should expect only one warning of this category
+            warn = w[0]
+            if issubclass(warn.category, NoConvergenceWarning):
+                new_msg = (
+                    f'\nGenericSectionCalculator::calculate_moment_curvature'
+                    f'\n\tNo convergence during computation of ultimate'
+                    f' curvature. Please check results properly.\n'
+                    f'{warn.message}'
+                )
+                warnings.warn(new_msg, NoConvergenceWarning)
+        # Find the yielding curvature
+        # This could not converge, let0's catch the warning to customize
+        # the message
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always', NoConvergenceWarning)
+            strain = self.find_equilibrium_fixed_pivot(
+                rotated_geom, n, yielding=True, max_iter=max_iter, tol=tol
+            )
+            chi_yield = strain[1]
+        if w:
+            # I should expect only one warning of this category
+            warn = w[0]
+            if issubclass(warn.category, NoConvergenceWarning):
+                new_msg = (
+                    f'\nGenericSectionCalculator::calculate_moment_curvature'
+                    f'\n\tNo convergence during computation of yielding'
+                    f' curvature. Please check results properly.\n'
+                    f'{warn.message}'
+                )
+                warnings.warn(new_msg, NoConvergenceWarning)
+
+        if chi_ultimate * chi_yield < 0:
+            # They cannot have opposite signs!
+            raise ValueError(
+                'Curvature at yield and ultimate cannot have opposite signs!'
+            )
+
+        # Make sure the sign of the first curvature matches the sign of the
+        # yield curvature
+        chi_first *= -1.0 if chi_first * chi_yield < 0 else 1.0
+
+        # The first curvature should be less than the yield curvature
+        if abs(chi_first) >= abs(chi_yield):
+            chi_first = chi_yield / num_pre_yield
+
+        # Define the array of curvatures
+        if abs(chi_ultimate) <= abs(chi_yield) + 1e-8:
+            # We don't want a plastic branch in the analysis
+            # this is done to speed up analysis
+            chi = np.linspace(chi_first, chi_yield, num_pre_yield)
+        else:
+            chi = np.concatenate(
+                (
+                    np.linspace(
+                        chi_first,
+                        chi_yield,
+                        num_pre_yield - 1,
+                        endpoint=False,
+                    ),
+                    np.linspace(chi_yield, chi_ultimate, num_post_yield + 1),
+                )
+            )
+        return chi
+
     def calculate_moment_curvature(
         self,
         theta: float = 0.0,
@@ -808,8 +907,10 @@ def calculate_moment_curvature(
         num_pre_yield: int = 10,
         num_post_yield: int = 10,
         chi: t.Optional[ArrayLike] = None,
+        max_iter: int = 100,
+        tol: float = 1e-2,
     ) -> s_res.MomentCurvatureResults:
-        """Calculates the moment-curvature relation for given inclination of
+        r"""Calculates the moment-curvature relation for given inclination of
         n.a. and axial load.
 
         Arguments:
@@ -830,6 +931,10 @@ def calculate_moment_curvature(
                 If chi is not None, chi_first, num_pre_yield and num_post_yield
                 are disregarded, and the provided chi is used directly in the
                 calculations.
+            max_iter (int): the maximum number of iterations in the iterative
+                process (default = 10).
+            tol (float): the tolerance for convergence test in terms of
+                $\deltaN_a - \deltaN_b$ (default = 1e-2).
 
         Returns:
             MomentCurvatureResults: The calculation results.
@@ -847,50 +952,15 @@ def calculate_moment_curvature(
             self._rotate_triangulated_data(-theta)
 
         if chi is None:
-            # Find ultimate curvature from the strain distribution
-            # corresponding to failure and equilibrium with external axial
-            # force
-            strain = self.find_equilibrium_fixed_pivot(rotated_geom, n)
-            chi_ultimate = strain[1]
-            # Find the yielding curvature
-            strain = self.find_equilibrium_fixed_pivot(
-                rotated_geom, n, yielding=True
+            chi = self._prepare_chi_array(
+                rotated_geom,
+                n,
+                num_pre_yield,
+                num_post_yield,
+                chi_first,
+                max_iter,
+                tol,
             )
-            chi_yield = strain[1]
-            if chi_ultimate * chi_yield < 0:
-                # They cannot have opposite signs!
-                raise ValueError(
-                    'curvature at yield and ultimate cannot have opposite '
-                    'signs!'
-                )
-
-            # Make sure the sign of the first curvature matches the sign of the
-            # yield curvature
-            chi_first *= -1.0 if chi_first * chi_yield < 0 else 1.0
-
-            # The first curvature should be less than the yield curvature
-            if abs(chi_first) >= abs(chi_yield):
-                chi_first = chi_yield / num_pre_yield
-
-            # Define the array of curvatures
-            if abs(chi_ultimate) <= abs(chi_yield) + 1e-8:
-                # We don't want a plastic branch in the analysis
-                # this is done to speed up analysis
-                chi = np.linspace(chi_first, chi_yield, num_pre_yield)
-            else:
-                chi = np.concatenate(
-                    (
-                        np.linspace(
-                            chi_first,
-                            chi_yield,
-                            num_pre_yield - 1,
-                            endpoint=False,
-                        ),
-                        np.linspace(
-                            chi_yield, chi_ultimate, num_post_yield + 1
-                        ),
-                    )
-                )
 
         # prepare results
         eps_a = np.zeros_like(chi)
@@ -905,31 +975,59 @@ def calculate_moment_curvature(
             # find the new position of neutral axis for mantaining equilibrium
             # store the information in the results object for the current
             # value of curvature
-            strain = self.find_equilibrium_fixed_curvature(
-                rotated_geom, n, curv, strain[0]
-            )
-            (
-                _,
-                My,
-                Mz,
-                _,
-            ) = self.integrator.integrate_strain_response_on_geometry(
-                geo=rotated_geom, strain=strain, tri=self.triangulated_data
-            )
-            # Rotate back to section CRS
-            T = np.array(
-                [
-                    [np.cos(theta), -np.sin(theta)],
-                    [np.sin(theta), np.cos(theta)],
-                ]
-            )
-            M = T @ np.array([[My], [Mz]])
-            eps_a[i] = strain[0]
-            my[i] = M[0, 0]
-            mz[i] = M[1, 0]
-            chi_mat = T @ np.array([[curv], [0]])
-            chi_y[i] = chi_mat[0, 0]
-            chi_z[i] = chi_mat[1, 0]
+            # DEBUG to remove
+            fact = 100 if i < 10 else 1
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter('always', NoConvergenceWarning)
+                strain = self.find_equilibrium_fixed_curvature(
+                    rotated_geom,
+                    n,
+                    curv,
+                    strain[0],
+                    max_iter=max_iter * fact,
+                    tol=tol,
+                )
+                (
+                    _,
+                    My,
+                    Mz,
+                    _,
+                ) = self.integrator.integrate_strain_response_on_geometry(
+                    geo=rotated_geom,
+                    strain=strain,
+                    tri=self.triangulated_data,
+                )
+                # Rotate back to section CRS
+                T = np.array(
+                    [
+                        [np.cos(theta), -np.sin(theta)],
+                        [np.sin(theta), np.cos(theta)],
+                    ]
+                )
+                M = T @ np.array([[My], [Mz]])
+                eps_a[i] = strain[0]
+                my[i] = M[0, 0]
+                mz[i] = M[1, 0]
+                chi_mat = T @ np.array([[curv], [0]])
+                chi_y[i] = chi_mat[0, 0]
+                chi_z[i] = chi_mat[1, 0]
+            if w:
+                # I should expect only one warning of this category
+                warn = w[0]
+                if issubclass(warn.category, NoConvergenceWarning):
+                    new_msg = (
+                        f'\nGenericSectionCalculator::calculate_moment_curvature\n'
+                        f'\tNo convergence during computation of step {i}.'
+                        f' The computation is stopped at the current step\n'
+                        f'{warn.message}'
+                    )
+                    warnings.warn(new_msg, NoConvergenceWarning)
+                    chi_y = chi_y[:i]
+                    chi_z = chi_z[:i]
+                    eps_a = eps_a[:i]
+                    my = my[:i]
+                    mz = mz[:i]
+                    break
 
         if self.triangulated_data is not None:
             # Rotate back also triangulated data!
diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py
index 83ab4d3a..742ff248 100644
--- a/tests/test_sections/test_generic_section.py
+++ b/tests/test_sections/test_generic_section.py
@@ -1486,3 +1486,91 @@ def test_rectangular_section_init_strain(fck, fyk, ductility_class):
     # Check they are all the same
     assert math.isclose(res_fiber.m_y, res_fiber_i.m_y, rel_tol=1e-3)
     assert math.isclose(res_marin.m_y, res_fiber_i.m_y, rel_tol=1e-3)
+
+
+@pytest.mark.parametrize(
+    'b, h, n_bars, diameter, fck, fyk', [(200, 400, 4, 16, 30, 500)]
+)
+def test_section_strength_warning(b, h, n_bars, diameter, fck, fyk):
+    """Test that a warning for no convergence is raised."""
+    # Create materials to use
+    concrete = ConcreteMC2010(fck)
+    steel = ReinforcementMC2010(fyk=fyk, Es=210000, ftk=fyk, epsuk=0.0675)
+
+    # Create the section
+    geo = RectangularGeometry(width=b, height=h, material=concrete)
+    geo = add_reinforcement_line(
+        geo,
+        coords_i=(-b / 2 + 40, -h / 2 + 40),
+        coords_j=(b / 2 - 40, -h / 2 + 40),
+        diameter=diameter,
+        material=steel,
+        n=n_bars,
+    )
+    geo = add_reinforcement_line(
+        geo,
+        coords_i=(-b / 2 + 40, h / 2 - 40),
+        coords_j=(b / 2 - 40, h / 2 - 40),
+        diameter=diameter,
+        material=steel,
+        n=n_bars,
+    )
+    section = GenericSection(geo)
+
+    # Compute bending strength without reaching equilibrium
+    # This is fictitiously tested forcing a low max_iter
+    with pytest.warns(NoConvergenceWarning):
+        section.section_calculator.calculate_bending_strength(max_iter=5)
+
+
+@pytest.mark.parametrize(
+    'b, h, n_bars, diameter, fck, fyk', [(200, 400, 4, 16, 30, 500)]
+)
+def test_section_moment_curvature_warning(b, h, n_bars, diameter, fck, fyk):
+    """Test that a warning for no convergence is raised."""
+    # Create materials to use
+    concrete = ConcreteMC2010(fck)
+    steel = ReinforcementMC2010(fyk=fyk, Es=210000, ftk=fyk, epsuk=0.0675)
+
+    # Create the section
+    geo = RectangularGeometry(width=b, height=h, material=concrete)
+    geo = add_reinforcement_line(
+        geo,
+        coords_i=(-b / 2 + 40, -h / 2 + 40),
+        coords_j=(b / 2 - 40, -h / 2 + 40),
+        diameter=diameter,
+        material=steel,
+        n=n_bars,
+    )
+    geo = add_reinforcement_line(
+        geo,
+        coords_i=(-b / 2 + 40, h / 2 - 40),
+        coords_j=(b / 2 - 40, h / 2 - 40),
+        diameter=diameter,
+        material=steel,
+        n=n_bars,
+    )
+    section = GenericSection(geo)
+    n_u = section.section_calculator.calculate_limit_axial_load()[0]
+    # Compute moment curvature with no warning
+    res_good = section.section_calculator.calculate_moment_curvature(
+        n=n_u * 0.5, max_iter=100
+    )
+    chi = res_good.chi_y.copy()
+
+    # Compute moment curvature without reaching equilibrium
+    # This is fictitiously tested forcing a low max_iter
+    with pytest.warns(NoConvergenceWarning):
+        res = section.section_calculator.calculate_moment_curvature(
+            n=n_u * 0.5, max_iter=8, chi=chi
+        )
+
+    # Check that the curvature arrays are the same up to the last
+    # convergence point
+    assert len(res.chi_y) < len(res_good.chi_y)
+    assert np.allclose(
+        res.chi_y,
+        res_good.chi_y[: len(res.chi_y)],
+        rtol=1e-5,
+        atol=1e-8,
+    )

From 4b39268046850ced1a1e865c642c92b2679c255f Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Mon, 20 Oct 2025 19:44:31 +0200
Subject: [PATCH 7/8] More robust version of _prefind_range and updated tests

---
 structuralcodes/sections/_generic.py        | 176 +++++++++++++++++---
 tests/test_sections/test_generic_section.py |   8 -
 2 files changed, 151 insertions(+), 33 deletions(-)

diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index 24d4fc8d..52b3ef03 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -471,33 +471,34 @@ def find_equilibrium_fixed_pivot(
         # Return the strain distribution
         return [eps_0, chi_c, 0]
 
-    def _prefind_range_curvature_equilibrium(
+    def _quick_exponential_find(
         self,
-        geom: CompoundGeometry,
-        n: float,
-        curv: float,
-        eps_0_a: float,
-        dn_a: float,
-        max_iter: int = 20,
-        max_restart_attemps: int = 20,
+        geom,
+        n,
+        curv,
+        max_iter,
+        max_restart_attempts,
+        eps_0_a,
+        dn_a,
+        sign,
     ):
-        """Perfind range where the curvature equilibrium is located.
-
-        This algorithms quickly finds a position of NA that guaranteed the
-        existence of at least one zero in the function dn vs. curv in order to
-        apply the bisection algorithm.
-        """
-        sign = -1 if dn_a > 0 else 1
+        restarts = 0
         found = False
         it = 0
-        restarts = 0
         delta = 1e-3
         # Use a growth factor for an exponential finding
         r = 2.0
         diverging = False
         diverging_steps = 0
-        while not found and it < max_iter and restarts < max_restart_attemps:
-            eps_0_b = eps_0_a + sign * delta * r ** (it)
+        while not found and it < max_iter and restarts < max_restart_attempts:
+            # If there has been at least 20% of maximum restart attempts,
+            # slow down the exponent growth to linear. This is helpful to
+            # solve cases where the solution is close to the limit axial
+            # load
+            exponent = (
+                it if restarts <= int(0.2 * max_restart_attempts) else 1.0
+            )
+            eps_0_b = eps_0_a + sign * delta * r ** (exponent)
             (
                 n_int,
                 _,
@@ -516,21 +517,148 @@ def _prefind_range_curvature_equilibrium(
                 if diverging:
                     # Count for how many steps we are diverging
                     diverging_steps += 1
-                # If we are consistently diverging for more than 10 steps,
-                # Restart the process with a small delta
-                if diverging_steps > 10:
+                # If we are consistently diverging for more than 5 steps,
+                # Restart the process with a smaller delta
+                if diverging_steps > 5:
                     delta /= 2
                     it = 0
                     restarts += 1
                     diverging = False
                     diverging_steps = 0
             it += 1
+        return eps_0_b, dn_b, it, found
+
+    def _slow_linear_find(
+        self,
+        geom,
+        n,
+        curv,
+        max_iter,
+        eps_0_a,
+        dn_a,
+        sign,
+    ):
+        found = False
+        it = 0
+        delta = abs(eps_0_a) / (max_iter * 10)
+        while True:
+            eps_0_b = eps_0_a + sign * delta * (it + 1)
+            (
+                n_int,
+                _,
+                _,
+                _,
+            ) = self.integrator.integrate_strain_response_on_geometry(
+                geom, [eps_0_b, curv, 0], tri=self.triangulated_data
+            )
+            dn_b = n_int - n
+            if dn_a * dn_b < 0:
+                found = True
+                break
+            it += 1
+            if it >= max_iter * 10:
+                break
+        return eps_0_b, dn_b, it, found
+
+    def _brute_force_find_minimum(
+        self,
+        geom,
+        n,
+        curv,
+        eps_0_a,
+        dn_a,
+        sign,
+    ):
+        found = False
+        suffix = ''
+        eps_0_a_attempt = eps_0_a
+        eps_0_b_attempt = eps_0_a + abs(eps_0_a) * sign
+        for i in range(10):
+            eps_0_attempts = np.linspace(eps_0_a_attempt, eps_0_b_attempt, 100)
+            dn_attempts = np.zeros_like(eps_0_attempts)
+            for j in range(len(eps_0_attempts)):
+                (
+                    n_int,
+                    _,
+                    _,
+                    _,
+                ) = self.integrator.integrate_strain_response_on_geometry(
+                    geom,
+                    [eps_0_attempts[j], curv, 0],
+                    tri=self.triangulated_data,
+                )
+                dn_attempts[j] = n_int - n
+            if dn_a > 0:
+                idx_min = np.argmin(dn_attempts)
+            else:
+                idx_min = np.argmax(dn_attempts)
+            eps_0_b = eps_0_attempts[idx_min]
+            dn_b = dn_attempts[idx_min]
+            eps_0_a_attempt = eps_0_attempts[max(idx_min - 2, 0)]
+            eps_0_b_attempt = eps_0_attempts[
+                min(idx_min + 2, len(eps_0_attempts))
+            ]
+            # Check if it is on the opposite side
+            if dn_a * dn_b < 0:
+                found = True
+                break
+
+        if not found:
+            # It is impossibile to find a solution for the given case
+            # The best solution we have is dn_b
+            suffix = (
+                '(No solution can be found because it seems dN never '
+                ' becomes of opposite sign).'
+                '\n\tThe best that we have found is with an unbalance of'
+                f' {dn_b:.3e} for eps_0_b = {eps_0_b:e}.'
+            )
+        return eps_0_b, dn_b, found, suffix
+
+    def _prefind_range_curvature_equilibrium(
+        self,
+        geom: CompoundGeometry,
+        n: float,
+        curv: float,
+        eps_0_a: float,
+        dn_a: float,
+        max_iter: int = 20,
+        max_restart_attemps: int = 20,
+    ):
+        """Perfind range where the curvature equilibrium is located.
+
+        This algorithms quickly finds a position of NA that guaranteed the
+        existence of at least one zero in the function dn vs. curv in order to
+        apply the bisection algorithm.
+        """
+        sign = -1 if dn_a > 0 else 1
+
+        eps_0_b, dn_b, it, found = self._quick_exponential_find(
+            geom, n, curv, max_iter, max_restart_attemps, eps_0_a, dn_a, sign
+        )
+
+        # As a further attempt, try going with small increments/decrements up
+        # to when it changes sign
+
+        if it >= max_iter and not found:
+            eps_0_b, dn_b, it, found = self._slow_linear_find(
+                geom, n, curv, max_iter, eps_0_a, dn_a, sign
+            )
+
+        # As a final attempt, try to compute for a range a certain size
+        # respect eps_0_a and compute dn for each one of them, find the
+        # minimum hoping it is on the opposite side of dn_a.
+
+        if it >= max_iter and not found:
+            eps_0_b, dn_b, found, suffix = self._brute_force_find_minimum(
+                geom, n, curv, eps_0_a, dn_a, sign
+            )
+
         if it >= max_iter and not found:
             msg = (
                 'GenericSectionCalculator::'
                 '_prefind_range_curvature_equilibrium\n\t'
                 'Maximum number of iterations reached.'
-                f' Last iteration reached a unbalance of {dn_b:.3e}'
+                f' Last iteration reached a unbalance of {dn_b:.3e} {suffix}'
             )
             warnings.warn(
                 message=msg,
@@ -975,8 +1103,6 @@ def calculate_moment_curvature(
             # find the new position of neutral axis for mantaining equilibrium
             # store the information in the results object for the current
             # value of curvature
-            # DEBUG to remove
-            fact = 100 if i < 10 else 1
             with warnings.catch_warnings(record=True) as w:
                 warnings.simplefilter('always', NoConvergenceWarning)
                 strain = self.find_equilibrium_fixed_curvature(
@@ -984,7 +1110,7 @@ def calculate_moment_curvature(
                     n,
                     curv,
                     strain[0],
-                    max_iter=max_iter * fact,
+                    max_iter=max_iter,
                     tol=tol,
                 )
                 (
diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py
index 742ff248..8e4f065f 100644
--- a/tests/test_sections/test_generic_section.py
+++ b/tests/test_sections/test_generic_section.py
@@ -100,14 +100,6 @@ def test_rectangular_section():
         res_mc_fiber.m_y[-1], res_mc_fiber_same_curvature.m_y[-1]
     )
 
-    # Calculate moment-curvature for a given array of curvatures, but with
-    # significant axial compression. This should raise a ValueError since we
-    # cannot find equilibrium.
-    with pytest.raises(ValueError):
-        sec.section_calculator.calculate_moment_curvature(
-            theta=0, n=0.95 * n_min_fiber, chi=res_mc_fiber.chi_y
-        )
-
     # check if axial limit forces are the same for marin and fiber
     assert math.isclose(n_min_marin, n_min_fiber, rel_tol=2e-2)
     assert math.isclose(n_max_marin, n_max_fiber, rel_tol=2e-2)

From 5309ee1a72b82f3c1eb897d95a466fe6298f1d5b Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Mon, 20 Oct 2025 21:52:49 +0200
Subject: [PATCH 8/8] warning mm domain calculation and test

---
 structuralcodes/sections/_generic.py        | 41 ++++++++++++-------
 tests/test_sections/test_generic_section.py | 45 +++++++++++++++++++++
 2 files changed, 72 insertions(+), 14 deletions(-)

diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index 52b3ef03..3747bc88 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -1060,7 +1060,7 @@ def calculate_moment_curvature(
                 are disregarded, and the provided chi is used directly in the
                 calculations.
             max_iter (int): the maximum number of iterations in the iterative
-                process (default = 10).
+                process (default = 100).
             tol (float): the tolerance for convergence test in terms of
                 $\deltaN_a - \deltaN_b$ (default = 1e-2).
 
@@ -1600,13 +1600,21 @@ def calculate_nmm_interaction_domain(
         return res
 
     def calculate_mm_interaction_domain(
-        self, n: float = 0, num_theta: int = 32
+        self,
+        n: float = 0,
+        num_theta: int = 32,
+        max_iter: int = 100,
+        tol: float = 1e-2,
     ) -> s_res.MMInteractionDomain:
-        """Calculate the My-Mz interaction domain.
+        r"""Calculate the My-Mz interaction domain.
 
         Arguments:
             n (float): Axial force, default = 0.
             n_theta (int): Number of discretization for theta, default = 32.
+            max_iter (int): the maximum number of iterations in the iterative
+                process (default = 10).
+            tol (float): the tolerance for convergence test in terms of
+                $\deltaN_a - \deltaN_b$ (default = 1e-2).
 
         Return:
             MMInteractionDomain: The calculation results.
@@ -1622,18 +1630,23 @@ def calculate_mm_interaction_domain(
         # Compute strength for given angle of NA
         for i, th in enumerate(res.theta):
             # If for some reason there is the warning of non convergence,
-            # catch it and raise an Exception because the domain is not
-            # correctly computed
-            try:
-                with warnings.catch_warnings():
-                    warnings.simplefilter('error', NoConvergenceWarning)
-                    res_bend_strength = self.calculate_bending_strength(
-                        theta=th, n=n
+            # catch it and raise it with a custom message
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter('always', NoConvergenceWarning)
+                res_bend_strength = self.calculate_bending_strength(
+                    theta=th, n=n, max_iter=max_iter, tol=tol
+                )
+            if w:
+                # I should expect only one warning of this category
+                warn = w[0]
+                if issubclass(warn.category, NoConvergenceWarning):
+                    new_msg = (
+                        '\nGenericSectionCalculator::'
+                        'calculate_mm_interaction_domain\n'
+                        f'\tNo convergence during computation with theta {th}.'
+                        f'\n{warn.message}'
                     )
-            except NoConvergenceWarning as e:
-                raise RuntimeError(
-                    'Convergence cannot be found, the algorithm is stopped.'
-                ) from e
+                    warnings.warn(new_msg, NoConvergenceWarning)
             # Save forces
             res.forces[i, 0] = n
             res.forces[i, 1] = res_bend_strength.m_y
diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py
index 8e4f065f..54c26e91 100644
--- a/tests/test_sections/test_generic_section.py
+++ b/tests/test_sections/test_generic_section.py
@@ -1566,3 +1566,48 @@ def test_section_moment_curvature_warning(b, h, n_bars, diameter, fck, fyk):
         rtol=1e-5,
         atol=1e-8,
     )
+
+
+@pytest.mark.parametrize(
+    'b, h, n_bars, diameter, fck, fyk', [(200, 400, 4, 16, 30, 500)]
+)
+def test_section_mm_domain_warning(b, h, n_bars, diameter, fck, fyk):
+    """Test that a warning for no convergence is raised."""
+    # Create materials to use
+    concrete = ConcreteMC2010(fck)
+    steel = ReinforcementMC2010(fyk=fyk, Es=210000, ftk=fyk, epsuk=0.0675)
+
+    # Create the section
+    geo = RectangularGeometry(width=b, height=h, material=concrete)
+    geo = add_reinforcement_line(
+        geo,
+        coords_i=(-b / 2 + 40, -h / 2 + 40),
+        coords_j=(b / 2 - 40, -h / 2 + 40),
+        diameter=diameter,
+        material=steel,
+        n=n_bars,
+    )
+    geo = add_reinforcement_line(
+        geo,
+        coords_i=(-b / 2 + 40, h / 2 - 40),
+        coords_j=(b / 2 - 40, h / 2 - 40),
+        diameter=diameter,
+        material=steel,
+        n=n_bars,
+    )
+    section = GenericSection(geo)
+    n_u = section.section_calculator.calculate_limit_axial_load()[0]
+    # Compute moment curvature with no warning
+    res_good = section.section_calculator.calculate_mm_interaction_domain(
+        n=n_u * 0.5, max_iter=100
+    )
+
+    # Compute moment curvature without reaching equilibrium
+    # This is fictitiously tested forcing a low max_iter
+    with pytest.warns(NoConvergenceWarning):
+        res = section.section_calculator.calculate_mm_interaction_domain(
+            n=n_u * 0.5, max_iter=5
+        )
+
+    # Check that the results arrays have the same size
+    assert len(res.m_y) == len(res_good.m_y)