From 4fee8f82391292cd00bfda3a13a22c136adcfe76 Mon Sep 17 00:00:00 2001 From: Diego Talledo <38036285+talledodiego@users.noreply.github.com> Date: Fri, 13 Mar 2026 14:34:18 +0100 Subject: [PATCH 1/2] Results for calculate strain profile --- structuralcodes/core/_section_results.py | 180 +++++++++++++++++- structuralcodes/sections/_generic.py | 53 +++++- tests/test_sections/test_generic_section.py | 20 +- .../test_rc_sections/test_rc_results.py | 76 ++++++++ 4 files changed, 311 insertions(+), 18 deletions(-) diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py index d5e3bf76..60cf12c8 100644 --- a/structuralcodes/core/_section_results.py +++ b/structuralcodes/core/_section_results.py @@ -7,7 +7,7 @@ from dataclasses import dataclass, field, fields import numpy as np -from numpy.typing import ArrayLike +from numpy.typing import ArrayLike, NDArray from shapely import Point @@ -800,6 +800,184 @@ def get_point_stress( ) +@dataclass(slots=True) +class StrainProfileResult: + """Class for storing the results from calculate_strain_profile method.""" + + # Solved generalized strains + eps_a: float = 0.0 # the axial strain at 0, 0 + chi_y: float = 0.0 # the curvature respect y axes + chi_z: float = 0.0 # the curvature respect z axes + + # target external loads assigned by user + n_ext: float = 0.0 # Target axial load acting at 0, 0 + m_y_ext: float = 0.0 # Target moment My + m_z_ext: float = 0.0 # Target moment Mz + + # Actual integrated loads at the converged strain state + n: float = 0.0 # Axial load acting at 0, 0 + m_y: float = 0.0 # Bending moment My + m_z: float = 0.0 # Bending moment Mz + + # Solver settings and diagnostics + tolerance: float = 0.0 + max_iter: int = 0 + used_initial_tangent: bool = False + iterations: int = 0 + converged: bool = False + residual: NDArray[np.float64] = field( + default_factory=lambda: np.zeros(3, dtype=float) + ) + + # Iteration history + residual_history: list[NDArray[np.float64]] = field(default_factory=list) + strain_history: list[NDArray[np.float64]] = field(default_factory=list) + + # For context store the section + section: t.Any = None # Note for future: if I want to type this I also have problem of circular import? #noqa E501 + # The detailed result data structure + detailed_result: SectionDetailedResultState = None + + @property + def residual_norm_history(self) -> t.List: + """Returns the history of residual norm.""" + return [float(np.linalg.norm(x)) for x in self.residual_history] + + @property + def delta_strain_history(self) -> t.List: + """Returns as a list the history of delta_strain.""" + return [ + self.strain_history[i] - self.strain_history[i - 1] + for i in range(1, len(self.strain_history)) + ] + + @property + def delta_strain_norm_history(self) -> t.List: + """Returns as a list the history of norm of delta_strain.""" + return [float(np.linalg.norm(x)) for x in self.delta_strain_history] + + @property + def response_history(self) -> t.List: + """Returns as a list the response (i.e. internal forces) history.""" + loads = np.array([self.n_ext, self.m_y_ext, self.m_z_ext]) + return [(loads - x) for x in self.residual_history] + + @property + def strain_plane(self) -> NDArray[np.float64]: + """Returns the strain profile as a numpy array.""" + return np.array([self.eps_a, self.chi_y, self.chi_z]) + + @property + def residual_norm(self) -> float: + """Returns the norm of the residual at last iteration.""" + return float(np.linalg.norm(self.residual)) + + def to_list(self) -> t.List: + """Returns the strain profile coefficients in a list.""" + return [self.eps_a, self.chi_y, self.chi_z] + + def create_detailed_result(self, num_points=1000): + """Create the detailed result object. + + Arguments: + num_points (int): Number of random points to sample for each + surface geometry (default = 1000). + """ + self.detailed_result = SectionDetailedResultState( + section=self.section, + eps_a=self.eps_a, + chi_y=self.chi_y, + chi_z=self.chi_z, + n=self.n, + m_y=self.m_y, + m_z=self.m_z, + num_points=num_points, + ) + + def get_point_strain( + self, + y: float, + z: float, + name: t.Optional[str] = None, + group_label: t.Optional[str] = None, + case_sensitive: bool = True, + all_results: bool = False, + ) -> float: + """Return the strain at a given point (y,z). + + Arguments: + y (float): The y-coordinate of the point. + z (float): The z-coordinate of the point. + name (str, optional): The name of the surface geometry to check. + group_label (str, optional): The group label of the surface + geometry to check. + case_sensitive (bool, optional): If True (default) the matching is + case sensitive. + all_results (bool): If True, return the strain for all geometries + that matches the filters, otherwise return the strain for the + first geometry that matches the filters (default False). + + Returns: + float: The strain at the given point, or None if the point is not + within any of the geometries that match the filters. + """ + return _get_point_response( + section=self.section, + eps_a=self.eps_a, + chi_y=self.chi_y, + chi_z=self.chi_z, + y=y, + z=z, + response_type='strain', + name=name, + group_label=group_label, + case_sensitive=case_sensitive, + all_results=all_results, + ) + + def get_point_stress( + self, + y: float, + z: float, + name: t.Optional[str] = None, + group_label: t.Optional[str] = None, + case_sensitive: bool = True, + all_results: bool = False, + ) -> float: + """Return the stress at a given point (y,z). + + Arguments: + y (float): The y-coordinate of the point. + z (float): The z-coordinate of the point. + name (str, optional): The pattern for filtering the geometries by + their name. + group_label (str, optional): The pattern for filtering the + geometries by their group_label. + case_sensitive (bool, optional): If True (default) the matching is + case sensitive. + all_results (bool): If True, return the stress for all geometries + that matches the filters, otherwise return the stress for the + first geometry that matches the filters (default False). + + Returns: + float: The strain at the given point, or None if the point is not + within any of the geometries that match the filters. + """ + return _get_point_response( + section=self.section, + eps_a=self.eps_a, + chi_y=self.chi_y, + chi_z=self.chi_z, + y=y, + z=z, + response_type='stress', + name=name, + group_label=group_label, + case_sensitive=case_sensitive, + all_results=all_results, + ) + + @dataclass class InteractionDomain: """Class for storing common data on all interaction domain results. diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py index eb9939c6..a26254ed 100644 --- a/structuralcodes/sections/_generic.py +++ b/structuralcodes/sections/_generic.py @@ -1438,7 +1438,7 @@ def calculate_strain_profile( initial: bool = False, max_iter: int = 10, tol: float = 1e-6, - ) -> t.List[float]: + ) -> s_res.StrainProfileResult: """Get the strain plane for a given axial force and biaxial bending. Args: @@ -1453,8 +1453,10 @@ def calculate_strain_profile( increment. Returns: - List(float): 3 floats: Axial strain at (0,0), and curvatures of the - section around y and z axes. + StrainProfileResult: A custom object of class StrainProfileResult + that contains the results. Note that the to_list() method + returns as a list only the strain profile coefficients. These + can be for instanced passed to `integrate_strain_profile`. """ # Get the gometry geom = self.section.geometry @@ -1477,7 +1479,12 @@ def calculate_strain_profile( # Calculate strain plane with Newton Rhapson Iterative method num_iter = 0 - strain = np.zeros(3) + strain = np.zeros(3, dtype=float) + residual = np.zeros(3, dtype=float) + converged = False + + residual_history = [] + strain_history = [] # Factorize once the stiffness matrix if using initial if initial: @@ -1495,6 +1502,10 @@ def calculate_strain_profile( response = np.array(self.integrate_strain_profile(strain=strain)) residual = loads - response + # Append to history variables + residual_history.append(residual.copy()) + strain_history.append(strain.copy()) + if initial: # Solve using the decomposed matrix delta_strain = lu_solve((lu, piv), residual) @@ -1512,16 +1523,42 @@ def calculate_strain_profile( # Solve using the current tangent stiffness delta_strain = np.linalg.solve(stiffness_tangent, residual) + # Check for convergence: + if np.linalg.norm(delta_strain) < tol: + converged = True + break + # Update the strain strain += delta_strain num_iter += 1 - # Check for convergence: - if np.linalg.norm(delta_strain) < tol: - break + # Create the results object + res = s_res.StrainProfileResult() + res.eps_a = float(strain[0]) + res.chi_y = float(strain[1]) + res.chi_z = float(strain[2]) + + res.n_ext = float(n) + res.m_y_ext = float(my) + res.m_z_ext = float(mz) + + res.n = float(response[0]) + res.m_y = float(response[1]) + res.m_z = float(response[2]) + + res.max_iter = max_iter + res.tolerance = tol + res.used_initial_tangent = initial + res.iterations = num_iter + 1 + res.converged = converged + res.residual = residual + res.residual_history = residual_history + res.strain_history = strain_history + + res.section = self.section if num_iter >= max_iter: raise StopIteration('Maximum number of iterations reached.') - return strain.tolist() + return res diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py index 044f2fe8..3ddbe6be 100644 --- a/tests/test_sections/test_generic_section.py +++ b/tests/test_sections/test_generic_section.py @@ -927,10 +927,10 @@ def test_strain_plane_calculation_elastic_Nmm(n, my, mz, Ec, b, h): # Compare assert np.allclose( - np.array(strain_fiber), expected_strain, rtol=1e-2, atol=1e-8 + strain_fiber.strain_plane, expected_strain, rtol=1e-2, atol=1e-8 ) assert np.allclose( - np.array(strain_marin), expected_strain, rtol=1e-2, atol=1e-8 + strain_marin.strain_plane, expected_strain, rtol=1e-2, atol=1e-8 ) @@ -1001,10 +1001,10 @@ def test_strain_plane_calculation_elastic_kNm(n, my, mz, Ec, b, h): # Compare assert np.allclose( - np.array(strain_fiber), expected_strain, rtol=5e-2, atol=1e-5 + strain_fiber.strain_plane, expected_strain, rtol=5e-2, atol=1e-5 ) assert np.allclose( - np.array(strain_marin), expected_strain, rtol=5e-2, atol=1e-5 + strain_marin.strain_plane, expected_strain, rtol=5e-2, atol=1e-5 ) @@ -1056,24 +1056,26 @@ def test_strain_plane_calculation_rectangular_rc(n, my, mz, fck, b, h): strain_marin = section.section_calculator.calculate_strain_profile( n, my, mz, tol=1e-7 ) - strain_marin = np.array(strain_marin) + strain_marin_array = np.array(strain_marin.to_list()) section = GenericSection(geo, integrator='fiber', mesh_size=0.0001) strain_fiber = section.section_calculator.calculate_strain_profile( n, my, mz, tol=1e-7 ) - strain_fiber = np.array(strain_fiber) + strain_fiber_array = np.array(strain_fiber.to_list()) # check that initial tangent gives the same solution at the end strain_fiber_initial = section.section_calculator.calculate_strain_profile( n, my, mz, tol=1e-7, initial=True, max_iter=80 ) - strain_fiber_initial = np.array(strain_fiber_initial) + strain_fiber_initial_array = np.array(strain_fiber_initial.to_list()) - assert np.allclose(strain_marin, strain_fiber, rtol=2e-2, atol=1e-6) assert np.allclose( - strain_fiber, strain_fiber_initial, rtol=2e-2, atol=1e-6 + strain_marin_array, strain_fiber_array, rtol=2e-2, atol=1e-6 + ) + assert np.allclose( + strain_fiber_array, strain_fiber_initial_array, rtol=2e-2, atol=1e-6 ) diff --git a/tests/test_sections/test_rc_sections/test_rc_results.py b/tests/test_sections/test_rc_sections/test_rc_results.py index 3535bcef..1d1fb8e9 100644 --- a/tests/test_sections/test_rc_sections/test_rc_results.py +++ b/tests/test_sections/test_rc_sections/test_rc_results.py @@ -1,5 +1,7 @@ """Testing section results for RC cross sections.""" +import math + import numpy as np import pytest @@ -275,7 +277,9 @@ def test_moment_curvature_result_point(simple_rc_section): # We know that the point (160, 160) is in the top reinforcement line y2, z2 = 160.0, 160.0 stress = res.get_point_stress(y=y2, z=z2, group_label='bottom') + strain = res.get_point_strain(y=y2, z=z2, group_label='bottom') assert stress is None + assert strain is None stress = res.get_point_stress(y=y2, z=z2, group_label='top') assert isinstance(stress, np.ndarray) @@ -340,3 +344,75 @@ def test_bending_strength_result(simple_rc_section): ) assert np.isclose(stress_top, res.detailed_result.point_data['stress'][-1]) + + +@pytest.mark.parametrize( + 'max_iter, tol, initial', + [ + (10, 1e-6, False), + (20, 1e-8, False), + (100, 1e-6, True), + ], +) +def test_calc_strain_profile_results( + simple_rc_section, max_iter, tol, initial +): + """Test that the calculate_strain_profile result.""" + section = simple_rc_section + + # External forces + n = 0 + my = -175 * 1e6 + mz = 0 + f_ext = np.array([n, my, mz]) + + # Get the bending strength result + strain_res = section.section_calculator.calculate_strain_profile( + *f_ext, max_iter=max_iter, tol=tol, initial=initial + ) + + # Check from result that the residual = f_ext - f_int + f_int = np.array([strain_res.n, strain_res.m_y, strain_res.m_z]) + assert np.allclose( + strain_res.residual, f_ext - f_int, rtol=1e-7, atol=1e-7 + ) + + # Check the the correct data are stored from convergence + assert math.isclose(strain_res.tolerance, tol) + assert strain_res.max_iter == max_iter + assert strain_res.used_initial_tangent == initial + + # check that length of residual_history is the same as iterations + assert len(strain_res.residual_history) == strain_res.iterations + + assert len(strain_res.response_history) == strain_res.iterations + assert len(strain_res.strain_history) == strain_res.iterations + assert len(strain_res.residual_norm_history) == strain_res.iterations + + # Check that length of delta_strain is iterations - 1 + assert len(strain_res.delta_strain_history) == strain_res.iterations - 1 + assert ( + len(strain_res.delta_strain_norm_history) == strain_res.iterations - 1 + ) + + # Check point stress and strain at the bottom reinforcement line + stress_bottom = strain_res.get_point_stress( + y=-160.0, z=-160.0, group_label='bottom' + ) + strain_bottom = strain_res.get_point_strain( + y=-160.0, z=-160.0, group_label='bottom' + ) + + # Get detailed result + strain_res.create_detailed_result() + + assert np.isclose( + stress_bottom, strain_res.detailed_result.point_data['stress'][0] + ) + assert np.isclose( + stress_bottom, + strain_res.detailed_result.point_data['material'][0].fyd(), + ) + assert np.isclose( + strain_bottom, strain_res.detailed_result.point_data['strain'][0] + ) From 2614700f862daff5a25f4597ec63e68bd3c96dd3 Mon Sep 17 00:00:00 2001 From: Diego Talledo <38036285+talledodiego@users.noreply.github.com> Date: Fri, 13 Mar 2026 14:50:26 +0100 Subject: [PATCH 2/2] fix linting --- structuralcodes/sections/_generic.py | 43 +++++++++++++--------------- 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py index a26254ed..e64611ea 100644 --- a/structuralcodes/sections/_generic.py +++ b/structuralcodes/sections/_generic.py @@ -1534,29 +1534,26 @@ def calculate_strain_profile( num_iter += 1 # Create the results object - res = s_res.StrainProfileResult() - res.eps_a = float(strain[0]) - res.chi_y = float(strain[1]) - res.chi_z = float(strain[2]) - - res.n_ext = float(n) - res.m_y_ext = float(my) - res.m_z_ext = float(mz) - - res.n = float(response[0]) - res.m_y = float(response[1]) - res.m_z = float(response[2]) - - res.max_iter = max_iter - res.tolerance = tol - res.used_initial_tangent = initial - res.iterations = num_iter + 1 - res.converged = converged - res.residual = residual - res.residual_history = residual_history - res.strain_history = strain_history - - res.section = self.section + res = s_res.StrainProfileResult( + eps_a=float(strain[0]), + chi_y=float(strain[1]), + chi_z=float(strain[2]), + n_ext=float(n), + m_y_ext=float(my), + m_z_ext=float(mz), + n=float(response[0]), + m_y=float(response[1]), + m_z=float(response[2]), + max_iter=max_iter, + tolerance=tol, + used_initial_tangent=initial, + iterations=num_iter + 1, + converged=converged, + residual=residual.copy(), + residual_history=residual_history.copy(), + strain_history=strain_history.copy(), + section=self.section, + ) if num_iter >= max_iter: raise StopIteration('Maximum number of iterations reached.')