feat: detailed results for calculate strain profile (#325)

* Results for calculate strain profile

* fix linting

Author: talledodiego

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.

structuralcodes/sections/_generic.py

@@ -1712,7 +1712,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:
@@ -1727,8 +1727,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
@@ -1751,7 +1753,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:
@@ -1769,6 +1776,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)
@@ -1786,14 +1797,37 @@ 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(
+            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:
             warnings.warn(
@@ -1804,4 +1838,4 @@ def calculate_strain_profile(
                 category=NoConvergenceWarning,
             )
 
-        return strain.tolist()
+        return res

tests/test_sections/test_generic_section.py

@@ -920,10 +920,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
     )
 
 
@@ -994,10 +994,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
     )
 
 
@@ -1049,24 +1049,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
     )