From 7ab139c5df657651879e88d13912a3aca2ec82b7 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 18:27:33 +0200
Subject: [PATCH 01/10] experimental idea: add methods to calculate fiber
 strains and stresses in UltimateBendingMomentResults

---
 structuralcodes/core/_section_results.py | 64 ++++++++++++++++++++++++
 structuralcodes/sections/_generic.py     |  1 +
 2 files changed, 65 insertions(+)

diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index 4a355b7e..a02d59a9 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -4,6 +4,7 @@
 
 from dataclasses import dataclass, field, fields
 
+import numpy as np
 from numpy.typing import ArrayLike
 
 
@@ -159,6 +160,69 @@ class UltimateBendingMomentResults:
     chi_z: float = 0  # the curvature corresponding to the ultimate moment
     eps_a: float = 0  # the axial strain at 0,0 corresponding to Mult
 
+    section = None
+
+    def get_fibers_strains(self, idx: int = None) -> ArrayLike:
+        """Return a dataset with the strains in each fiber."""
+        if self.section is None:
+            return None
+
+        strain = [self.eps_a, self.chi_y, self.chi_z]
+
+        y = []
+        z = []
+        A = []
+        eps = []
+        if idx is not None:
+            tr = self.section.section_calculator.triangulated_data[idx]
+            return (
+                tr[0],
+                tr[1],
+                tr[2],
+                strain[0] - strain[2] * tr[0] + strain[1] * tr[1],
+            )
+        for tr in self.section.section_calculator.triangulated_data:
+            # All have the same material
+            y.append(tr[0])
+            z.append(tr[1])
+            A.append(tr[2])
+            strains = strain[0] - strain[2] * tr[0] + strain[1] * tr[1]
+            eps.append(strains)
+
+        return np.hstack(y), np.hstack(z), np.hstack(A), np.hstack(eps)
+
+    def get_fibers_stresses(self, idx: int = None) -> ArrayLike:
+        """Return a dataset with the stresses in each fiber."""
+        if self.section is None:
+            return None
+
+        strain = [self.eps_a, self.chi_y, self.chi_z]
+
+        y = []
+        z = []
+        A = []
+        sig = []
+        if idx is not None:
+            tr = self.section.section_calculator.triangulated_data[idx]
+            return (
+                tr[0],
+                tr[1],
+                tr[2],
+                tr[3].get_stress(
+                    strain[0] - strain[2] * tr[0] + strain[1] * tr[1]
+                ),
+            )
+        for tr in self.section.section_calculator.triangulated_data:
+            # All have the same material
+            y.append(tr[0])
+            z.append(tr[1])
+            A.append(tr[2])
+            strains = strain[0] - strain[2] * tr[0] + strain[1] * tr[1]
+            stresses = tr[3].get_stress(strains)
+            sig.append(stresses)
+
+        return np.hstack(y), np.hstack(z), np.hstack(A), np.hstack(sig)
+
 
 @dataclass
 class InteractionDomain:
diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index 75a0440f..87406826 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -770,6 +770,7 @@ def calculate_bending_strength(
         res.eps_a = strain[0]
         res.m_y = M[0, 0]
         res.m_z = M[1, 0]
+        res.section = self.section
 
         return res
 

From d80fea9cc3a62e8df035739166119e695d471c15 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 20:58:36 +0200
Subject: [PATCH 02/10] get random points from Surface Geometry

---
 structuralcodes/geometry/_geometry.py | 87 +++++++++++++++++++++++++++
 1 file changed, 87 insertions(+)

diff --git a/structuralcodes/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py
index 6398b8a2..30acde04 100644
--- a/structuralcodes/geometry/_geometry.py
+++ b/structuralcodes/geometry/_geometry.py
@@ -7,6 +7,7 @@
 from math import atan2
 
 import numpy as np
+import triangle
 from numpy.typing import ArrayLike
 from shapely import affinity
 from shapely.geometry import (
@@ -454,6 +455,92 @@ def polygon(self) -> Polygon:
         """Returns the Shapely Polygon."""
         return self._polygon
 
+    def random_points_within(self, num_points: int = 100) -> np.ndarray:
+        """Returns coordinates of random points within the polygon.
+
+        Arguments:
+            num_points (int): Number of random points to generate.
+
+        Returns:
+            x, y (ndarray, ndarray): Arrays with the x and y coordinates of the
+            random points within the polygon.
+        """
+
+        def _prepare_triangulation_data(poly: Polygon) -> dict[str:ArrayLike]:
+            # Create the tri dictionary
+            tri: dict[str:ArrayLike] = {}
+            # 1. External boundary process
+            # 1a. Get vertices, skipping the last one
+            vertices = np.column_stack(poly.exterior.xy)[:-1, :]
+            n_vertices = vertices.shape[0]
+            # 1b. Create segments
+            node_i = np.arange(n_vertices)
+            node_j = np.roll(node_i, -1)
+            segments = np.column_stack((node_i, node_j))
+
+            # 2. Process holes
+            holes = []
+            for interior in poly.interiors:
+                # 2a. Get vertices, skipping the last one
+                vertices_int = np.column_stack(interior.xy)[:-1, :]
+                n_vertices_int = vertices_int.shape[0]
+                # 2b. Create segments
+                node_i = np.arange(n_vertices_int) + n_vertices
+                node_j = np.roll(node_i, -1)
+                segments_int = np.column_stack((node_i, node_j))
+                c = Polygon(interior)
+                holes.append([c.centroid.x, c.centroid.y])
+                # Append to the global arrays
+                vertices = np.vstack((vertices, vertices_int))
+                segments = np.vstack((segments, segments_int))
+                n_vertices += n_vertices_int
+            # Return the dictionary with data for triangulate
+            tri['vertices'] = vertices
+            tri['segments'] = segments
+            if len(holes) > 0:
+                tri['holes'] = holes
+            return tri
+
+        tri = _prepare_triangulation_data(self.polygon)
+        triangles = triangle.triangulate(tri, 'p')
+
+        xs = np.array([])
+        ys = np.array([])
+        for tr in triangles['triangles']:
+            # Get vertices for the triangle
+            Ax = triangles['vertices'][tr[0]][0]
+            Ay = triangles['vertices'][tr[0]][1]
+            Bx = triangles['vertices'][tr[1]][0]
+            By = triangles['vertices'][tr[1]][1]
+            Cx = triangles['vertices'][tr[2]][0]
+            Cy = triangles['vertices'][tr[2]][1]
+            # area of the triangle
+            a = Ax * By - Ay * Bx
+            a += Bx * Cy - By * Cx
+            a += Cx * Ay - Cy * Cx
+            a = abs(a) * 0.5
+            # number of points in this triangle (at least 1)
+            n = max(1, int(num_points * a / self.area))
+            # generate random points in the triangle
+            r1 = np.random.uniform(0, 1, n)
+            r2 = np.random.uniform(0, 1, n)
+            x = (
+                (1 - np.sqrt(r1)) * Ax
+                + (np.sqrt(r1) * (1 - r2)) * Bx
+                + (r2 * np.sqrt(r1)) * Cx
+            )
+            y = (
+                (1 - np.sqrt(r1)) * Ay
+                + (np.sqrt(r1) * (1 - r2)) * By
+                + (r2 * np.sqrt(r1)) * Cy
+            )
+
+            # Concatenate the new points
+            xs = np.concatenate((xs, x))
+            ys = np.concatenate((ys, y))
+
+        return xs, ys
+
     def calculate_extents(self) -> t.Tuple[float, float, float, float]:
         """Calculate extents of SurfaceGeometry.
 

From 91c2f7f72964fba117392c1cb641fd9a52e5a239 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 21:06:48 +0200
Subject: [PATCH 03/10] bugfix in random_points_within

---
 structuralcodes/geometry/_geometry.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/structuralcodes/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py
index 30acde04..b5d3649f 100644
--- a/structuralcodes/geometry/_geometry.py
+++ b/structuralcodes/geometry/_geometry.py
@@ -517,7 +517,7 @@ def _prepare_triangulation_data(poly: Polygon) -> dict[str:ArrayLike]:
             # area of the triangle
             a = Ax * By - Ay * Bx
             a += Bx * Cy - By * Cx
-            a += Cx * Ay - Cy * Cx
+            a += Cx * Ay - Cy * Ax
             a = abs(a) * 0.5
             # number of points in this triangle (at least 1)
             n = max(1, int(num_points * a / self.area))

From a407fcd9d632b75219246d16e16b4961b7a06885 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 21:18:41 +0200
Subject: [PATCH 04/10] included tests random points within

---
 tests/test_geometry/test_geometry.py | 101 ++++++++++++++++++++++++++-
 1 file changed, 100 insertions(+), 1 deletion(-)

diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index 8f59ba92..3d8ab151 100644
--- a/tests/test_geometry/test_geometry.py
+++ b/tests/test_geometry/test_geometry.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 import pytest
-from shapely import LineString, MultiLineString, MultiPolygon, Polygon
+from shapely import LineString, MultiLineString, MultiPolygon, Point, Polygon
 from shapely.affinity import translate
 from shapely.testing import assert_geometries_equal
 
@@ -735,3 +735,102 @@ def test_mirror_geometry():
         geo_m.point_geometries[0].x,
         abs_tol=1e-5,
     )
+
+
+def test_random_points_within():
+    """Test random_points_within method generates points inside geometry."""
+    # Create a simple rectangular geometry
+    width = 100
+    height = 50
+    concrete = ElasticMaterial(E=30_000, density=2400)
+
+    geo = RectangularGeometry(width, height, concrete, True)
+
+    # Test with different numbers of points
+    for num_points in [10, 50, 100]:
+        xs, ys = geo.random_points_within(num_points)
+
+        # Check that we get arrays
+        assert isinstance(xs, np.ndarray)
+        assert isinstance(ys, np.ndarray)
+
+        # Check that arrays have the same length
+        assert len(xs) == len(ys)
+
+        # Check that we got roughly the expected number of points
+        # (allow some variation)
+        assert len(xs) > 0
+        # Allow some overshoot due to triangulation
+        assert len(xs) <= num_points * 1.5
+
+        # Check that all points are within the geometry bounds
+        min_x, max_x, min_y, max_y = geo.calculate_extents()
+
+        assert np.all(xs >= min_x)
+        assert np.all(xs <= max_x)
+        assert np.all(ys >= min_y)
+        assert np.all(ys <= max_y)
+
+        # Check that points are actually inside the polygon using shapely
+        for x, y in zip(xs, ys):
+            point = Point(x, y)
+            assert geo.polygon.contains(point) or geo.polygon.touches(point)
+
+
+def test_random_points_within_complex_geometry():
+    """Test random_points_within with a geometry containing holes."""
+    concrete = ElasticMaterial(E=30_000, density=2400)
+
+    # Create outer rectangle
+    outer_vertices = [(-50, -25), (50, -25), (50, 25), (-50, 25)]
+    outer_polygon = Polygon(outer_vertices)
+
+    # Create inner hole (smaller rectangle)
+    hole_vertices = [(-20, -10), (20, -10), (20, 10), (-20, 10)]
+    hole_polygon = Polygon(hole_vertices)
+
+    # Create polygon with hole
+    polygon_with_hole = Polygon(
+        outer_polygon.exterior, [hole_polygon.exterior]
+    )
+
+    geo = SurfaceGeometry(polygon_with_hole, concrete, True)
+
+    # Generate random points
+    num_points = 100
+    xs, ys = geo.random_points_within(num_points)
+
+    # Check basic properties
+    assert isinstance(xs, np.ndarray)
+    assert isinstance(ys, np.ndarray)
+    assert len(xs) == len(ys)
+    assert len(xs) > 0
+
+    # Check that all points are within the outer boundary but not in the hole
+    for x, y in zip(xs, ys):
+        point = Point(x, y)
+        # Point should be in the geometry (which excludes the hole)
+        assert geo.polygon.contains(point) or geo.polygon.touches(point)
+        # Point should not be inside the hole
+        assert not hole_polygon.contains(point)
+
+
+def test_random_points_within_reproducibility():
+    """Test that random_points_within can be made reproducible with seed."""
+    width = 100
+    height = 50
+    concrete = ElasticMaterial(E=30_000, density=2400)
+
+    geo = RectangularGeometry(width, height, concrete, True)
+
+    # Set seed for reproducibility
+    np.random.seed(42)
+    xs1, ys1 = geo.random_points_within(50)
+
+    # Reset seed to same value
+    np.random.seed(42)
+    xs2, ys2 = geo.random_points_within(50)
+
+    # Results should be identical
+    np.testing.assert_array_equal(xs1, xs2)
+    np.testing.assert_array_equal(ys1, ys2)

From 7634c9bf76971be07104db4e1405c4f3db92967a Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 21:39:38 +0200
Subject: [PATCH 05/10] feat: add method to retrieve coordinates of point
 geometries as numpy array

---
 structuralcodes/geometry/_geometry.py | 24 +++++++
 tests/test_geometry/test_geometry.py  | 97 ++++++++++++++++++++++++++-
 2 files changed, 120 insertions(+), 1 deletion(-)

diff --git a/structuralcodes/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py
index b5d3649f..386dd3aa 100644
--- a/structuralcodes/geometry/_geometry.py
+++ b/structuralcodes/geometry/_geometry.py
@@ -887,6 +887,30 @@ def area(self) -> float:
             area += geo.area
         return area
 
+    def get_coordinates_point_geometries(
+        self, group_label: t.Optional[str] = None
+    ) -> np.ndarray:
+        """Returns the coordinates of the point geometries as an array.
+
+        Arguments:
+            group_label (Optional(str)): If provided, only point geometries
+                with the given group_label are returned.
+
+        Returns:
+            ndarray: An array of shape (n, 2) with the coordinates of the
+            point geometries.
+        """
+        coords = np.zeros((len(self.point_geometries), 3))
+        for i, pg in enumerate(self.point_geometries):
+            if group_label is not None and pg.group_label != group_label:
+                continue
+            coords[i, 0] = pg.x
+            coords[i, 1] = pg.y
+            coords[i, 2] = 1
+        # Only keep rows where the last column is 1
+        filtered_coords = coords[coords[:, 2] == 1]
+        return filtered_coords[:, 0:2]
+
     def calculate_extents(self) -> t.Tuple[float, float, float, float]:
         """Calculate extents of CompundGeometry.
 
diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index 3d8ab151..2d514015 100644
--- a/tests/test_geometry/test_geometry.py
+++ b/tests/test_geometry/test_geometry.py
@@ -28,7 +28,9 @@
     ElasticPlastic,
     ParabolaRectangle,
 )
-from structuralcodes.materials.reinforcement import ReinforcementMC2010
+from structuralcodes.materials.reinforcement import (
+    ReinforcementMC2010,
+)
 
 
 # Test create line
@@ -834,3 +836,96 @@ def test_random_points_within_reproducibility():
     # Results should be identical
     np.testing.assert_array_equal(xs1, xs2)
     np.testing.assert_array_equal(ys1, ys2)
+
+
+def test_get_coordinates_point_geometries():
+    """Test get_coordinates_point_geometries method with filtering."""
+    # Create materials
+    concrete = ConcreteMC2010(45)
+    reinforcement = ReinforcementMC2010(450, 200000, 450, 0.075)
+    prestress_reinf = ReinforcementMC2010(
+        fyk=1620, Es=200000, ftk=1620, epsuk=0.02, initial_strain=0.0054
+    )
+
+    # Create the polygon
+    poly = Polygon(
+        [
+            (250, 0),
+            (250, 140),
+            (70, 200),
+            (70, 920),
+            (250, 980),
+            (250, 1120),
+            (-250, 1120),
+            (-250, 980),
+            (-70, 920),
+            (-70, 200),
+            (-250, 140),
+            (-250, 0),
+        ]
+    )
+
+    # Create surface geometry
+    geo = SurfaceGeometry(poly=poly, material=concrete)
+
+    # Add reinforcement lines
+    position = [40, 110, 1010, 1080]
+    number = [4, 2, 2, 4]
+    for y, n in zip(position, number):
+        geo = add_reinforcement_line(
+            geo=geo,
+            coords_i=(-210, y),
+            coords_j=(210, y),
+            diameter=10,
+            n=n,
+            material=reinforcement,
+            group_label='reinforcement',
+        )
+
+    # Add prestress reinforcement
+    area = 840
+    d = (area * 4 / np.pi) ** 0.5
+    geo = add_reinforcement_line(
+        geo=geo,
+        coords_i=(0, 70),
+        coords_j=(0, 210),
+        diameter=d,
+        n=2,
+        material=prestress_reinf,
+        group_label='prestress_reinforcement',
+    )
+
+    # Test without group_label filter - should return all points
+    all_coords = geo.get_coordinates_point_geometries()
+    total_reinforcement_points = sum(number)  # 12 regular reinforcement
+    total_prestress_points = 2  # 2 prestress reinforcement
+    expected_total = total_reinforcement_points + total_prestress_points
+    assert all_coords.shape[0] == expected_total
+    assert all_coords.shape[1] == 2  # x, y coordinates
+
+    # Test with "reinforcement" filter
+    reinforcement_coords = geo.get_coordinates_point_geometries(
+        group_label='reinforcement'
+    )
+    assert reinforcement_coords.shape[0] == total_reinforcement_points
+    assert reinforcement_coords.shape[1] == 2
+
+    # Test with "prestress_reinforcement" filter
+    prestress_coords = geo.get_coordinates_point_geometries(
+        group_label='prestress_reinforcement'
+    )
+    assert prestress_coords.shape[0] == total_prestress_points
+    assert prestress_coords.shape[1] == 2
+
+    # Test with wrong filter - should return empty array
+    wrong_coords = geo.get_coordinates_point_geometries(
+        group_label='prestress'
+    )
+    assert wrong_coords.shape[0] == 0
+    assert wrong_coords.shape[1] == 2  # Still 2 columns even if empty
+
+    # Verify that the sum of filtered results equals total
+    assert (
+        reinforcement_coords.shape[0] + prestress_coords.shape[0]
+        == all_coords.shape[0]
+    )

From bca0d738d56a87233741f4afc21141463144ac07 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 22:48:01 +0200
Subject: [PATCH 06/10] attempt_detailed_results

---
 structuralcodes/core/_section_results.py | 82 ++++++++++++++++++++++++
 structuralcodes/geometry/_geometry.py    | 31 ++++-----
 2 files changed, 98 insertions(+), 15 deletions(-)

diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index a02d59a9..ad7ec35e 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -146,6 +146,74 @@ class MomentCurvatureResults:
     # The stresses can be recomputed if needed on the fly? Or storing them?
 
 
+class SectionDetailedResultState:
+    """Something."""
+
+    def __init__(self, section, eps_a, chi_y, chi_z, n, m_y, m_z):
+        self.eps_a = eps_a
+        self.chi_y = chi_y
+        self.chi_z = chi_z
+        self.n = n
+        self.m_y = m_y
+        self.m_z = m_z
+        self.section = section
+        self._create_data_structure()
+
+    def _create_data_structure(self):
+        # Data containers
+        surface_data = []
+        point_data = {}
+
+        # SurfaceGeometry random points
+        for surf in self.section.geometry.geometries:
+            # Use surf.random_points_within to get random points
+            y, z = surf.random_points_within(num_points=1000)
+            strain = self.eps_a - self.chi_z * y + self.chi_y * z
+            stress = surf.material.constitutive_law.get_stress(strain)
+            surface_data.append(
+                {
+                    'y': y,
+                    'z': z,
+                    'strain': strain,
+                    'stress': stress,
+                    'name': getattr(surf, 'name', None),
+                    'group_label': getattr(surf, 'group_label', None),
+                }
+            )
+
+        # Collect all group_label values from point_geometries and
+        # find unique ones
+        unique_labels = set(
+            [pg.group_label for pg in self.section.geometry.point_geometries]
+        )
+
+        for label in unique_labels:
+            # There should always be at least one label, which is None
+            coords, mats = self.section.geometry.get_point_geometries(
+                group_label=label
+            )
+            strain = (
+                self.eps_a
+                - self.chi_z * coords[:, 0]
+                + self.chi_y * coords[:, 1]
+            )
+            # Find indices where all mats are the same
+            unique_mats = set(tuple(mat for mat in mats))
+            for mat in unique_mats:
+                mask = [m == mat for m in mats]
+                selected_coords = coords[mask]
+                stress = mat.constitutive_law.get_stress(strain[mask])
+                point_data[label] = {
+                    'coords': selected_coords,
+                    'strain': strain[mask],
+                    'stress': stress,
+                }
+
+        # Store results for later use (e.g., plotting)
+        self.surface_data = surface_data
+        self.point_data = point_data
+
+
 @dataclass
 class UltimateBendingMomentResults:
     """Class for storing the ultimate bending moment computation for a given
@@ -162,6 +230,20 @@ class UltimateBendingMomentResults:
 
     section = None
 
+    detailed_result: SectionDetailedResultState = None
+
+    def create_detailed_result(self):
+        """Create the detailed result object."""
+        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,
+        )
+
     def get_fibers_strains(self, idx: int = None) -> ArrayLike:
         """Return a dataset with the strains in each fiber."""
         if self.section is None:
diff --git a/structuralcodes/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py
index 386dd3aa..29f9726b 100644
--- a/structuralcodes/geometry/_geometry.py
+++ b/structuralcodes/geometry/_geometry.py
@@ -887,29 +887,30 @@ def area(self) -> float:
             area += geo.area
         return area
 
-    def get_coordinates_point_geometries(
+    def get_point_geometries(
         self, group_label: t.Optional[str] = None
-    ) -> np.ndarray:
-        """Returns the coordinates of the point geometries as an array.
+    ) -> t.Tuple[np.ndarray, t.List[Material]]:
+        """Returns the coordinates and the material of the point geometries.
+         The coordinates are returned as a numpy array.
 
         Arguments:
             group_label (Optional(str)): If provided, only point geometries
                 with the given group_label are returned.
 
         Returns:
-            ndarray: An array of shape (n, 2) with the coordinates of the
-            point geometries.
+            Tuple[ndarray, List[Material]]: Coordinates array of shape (n, 2)
+            and list of materials.
+
         """
-        coords = np.zeros((len(self.point_geometries), 3))
-        for i, pg in enumerate(self.point_geometries):
-            if group_label is not None and pg.group_label != group_label:
-                continue
-            coords[i, 0] = pg.x
-            coords[i, 1] = pg.y
-            coords[i, 2] = 1
-        # Only keep rows where the last column is 1
-        filtered_coords = coords[coords[:, 2] == 1]
-        return filtered_coords[:, 0:2]
+        coords = []
+        materials = []
+
+        for pg in self.point_geometries:
+            if group_label is None or pg.group_label == group_label:
+                coords.append([pg.x, pg.y])
+                materials.append(pg.material)
+
+        return np.array(coords), materials
 
     def calculate_extents(self) -> t.Tuple[float, float, float, float]:
         """Calculate extents of CompundGeometry.

From ab16f599d557e84f87353851ab138f8850f262fb Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 23:15:22 +0200
Subject: [PATCH 07/10] possible specifying number of points

---
 structuralcodes/core/_section_results.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index ad7ec35e..26770380 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -149,7 +149,9 @@ class MomentCurvatureResults:
 class SectionDetailedResultState:
     """Something."""
 
-    def __init__(self, section, eps_a, chi_y, chi_z, n, m_y, m_z):
+    def __init__(
+        self, section, eps_a, chi_y, chi_z, n, m_y, m_z, num_points=1000
+    ):
         self.eps_a = eps_a
         self.chi_y = chi_y
         self.chi_z = chi_z
@@ -157,9 +159,9 @@ def __init__(self, section, eps_a, chi_y, chi_z, n, m_y, m_z):
         self.m_y = m_y
         self.m_z = m_z
         self.section = section
-        self._create_data_structure()
+        self._create_data_structure(num_points=num_points)
 
-    def _create_data_structure(self):
+    def _create_data_structure(self, num_points=1000):
         # Data containers
         surface_data = []
         point_data = {}
@@ -167,7 +169,7 @@ def _create_data_structure(self):
         # SurfaceGeometry random points
         for surf in self.section.geometry.geometries:
             # Use surf.random_points_within to get random points
-            y, z = surf.random_points_within(num_points=1000)
+            y, z = surf.random_points_within(num_points=num_points)
             strain = self.eps_a - self.chi_z * y + self.chi_y * z
             stress = surf.material.constitutive_law.get_stress(strain)
             surface_data.append(

From e40b097ef75c9833cf29707546ae0f51d7293b4e Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sat, 4 Oct 2025 23:18:31 +0200
Subject: [PATCH 08/10] bugfix num_points argument

---
 structuralcodes/core/_section_results.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index 26770380..a8404345 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -234,7 +234,7 @@ class UltimateBendingMomentResults:
 
     detailed_result: SectionDetailedResultState = None
 
-    def create_detailed_result(self):
+    def create_detailed_result(self, num_points=1000):
         """Create the detailed result object."""
         self.detailed_result = SectionDetailedResultState(
             section=self.section,
@@ -244,6 +244,7 @@ def create_detailed_result(self):
             n=self.n,
             m_y=self.m_y,
             m_z=self.m_z,
+            num_points=num_points,
         )
 
     def get_fibers_strains(self, idx: int = None) -> ArrayLike:

From 29ea9673e5eef4992d5bb53e35509ba5c4b800a3 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sun, 5 Oct 2025 10:58:27 +0200
Subject: [PATCH 09/10] implement detailed results to moment_curvature

---
 structuralcodes/core/_section_results.py | 56 +++++++++++++++++++++++-
 structuralcodes/sections/_generic.py     |  1 +
 2 files changed, 55 insertions(+), 2 deletions(-)

diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index a8404345..b31d60d6 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -142,8 +142,60 @@ class MomentCurvatureResults:
     m_y: ArrayLike = None  # the moment
     m_z: ArrayLike = None  # the moment
 
-    # The strains can be reconstructed at each step from chi and eps_axial
-    # The stresses can be recomputed if needed on the fly? Or storing them?
+    section = None
+
+    detailed_result: SectionDetailedResultState = None
+    current_step: int = None
+    num_points: int = None
+
+    def create_detailed_result(self, num_points=1000):
+        """Create the detailed result object."""
+        self.detailed_result = SectionDetailedResultState(
+            section=self.section,
+            eps_a=self.eps_axial[0],
+            chi_y=self.chi_y[0],
+            chi_z=self.chi_z[0],
+            n=self.n,
+            m_y=self.m_y[0],
+            m_z=self.m_z[0],
+            num_points=num_points,
+        )
+        self.current_step = 0
+        self.num_points = num_points
+
+    def next_step(self):
+        """Advance to the next step in the detailed result."""
+        if self.detailed_result is None:
+            return
+        if self.current_step < len(self.m_y) - 1:
+            self.current_step += 1
+            self.detailed_result = SectionDetailedResultState(
+                section=self.section,
+                eps_a=self.eps_axial[self.current_step],
+                chi_y=self.chi_y[self.current_step],
+                chi_z=self.chi_z[self.current_step],
+                n=self.n,
+                m_y=self.m_y[self.current_step],
+                m_z=self.m_z[self.current_step],
+                num_points=self.num_points,
+            )
+
+    def previous_step(self):
+        """Go back to the previous step in the detailed result."""
+        if self.detailed_result is None:
+            return
+        if self.current_step > 0:
+            self.current_step -= 1
+            self.detailed_result = SectionDetailedResultState(
+                section=self.section,
+                eps_a=self.eps_axial[self.current_step],
+                chi_y=self.chi_y[self.current_step],
+                chi_z=self.chi_z[self.current_step],
+                n=self.n,
+                m_y=self.m_y[self.current_step],
+                m_z=self.m_z[self.current_step],
+                num_points=self.num_points,
+            )
 
 
 class SectionDetailedResultState:
diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index 87406826..786ec979 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -813,6 +813,7 @@ def calculate_moment_curvature(
 
         # Create an empty response object
         res = s_res.MomentCurvatureResults()
+        res.section = self.section
         res.n = n
         # Rotate the section of angle theta
         rotated_geom = self.section.geometry.rotate(-theta)

From 20631a40112f4505544e1c19a2e76d7043132ae1 Mon Sep 17 00:00:00 2001
From: Diego Talledo <38036285+talledodiego@users.noreply.github.com>
Date: Sun, 5 Oct 2025 18:00:32 +0200
Subject: [PATCH 10/10] feat: add seed parameter for reproducibility in random
 point generation

---
 structuralcodes/core/_section_results.py | 62 ++++++++++++++++++++++--
 structuralcodes/geometry/_geometry.py    | 19 ++++++--
 tests/test_geometry/test_geometry.py     | 23 ++++-----
 3 files changed, 83 insertions(+), 21 deletions(-)

diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index b31d60d6..142fbbb5 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -145,11 +145,19 @@ class MomentCurvatureResults:
     section = None
 
     detailed_result: SectionDetailedResultState = None
+    seed: int = None
     current_step: int = None
     num_points: int = None
 
     def create_detailed_result(self, num_points=1000):
-        """Create the detailed result object."""
+        """Create the detailed result object.
+
+        Arguments:
+            num_points (int): Number of random points to sample in each
+                surface geometry (default = 1000).
+        """
+        if self.seed is None:
+            self.seed = np.random.randint(1, 100, 1)[0].item()
         self.detailed_result = SectionDetailedResultState(
             section=self.section,
             eps_a=self.eps_axial[0],
@@ -159,6 +167,7 @@ def create_detailed_result(self, num_points=1000):
             m_y=self.m_y[0],
             m_z=self.m_z[0],
             num_points=num_points,
+            seed=self.seed,
         )
         self.current_step = 0
         self.num_points = num_points
@@ -178,6 +187,7 @@ def next_step(self):
                 m_y=self.m_y[self.current_step],
                 m_z=self.m_z[self.current_step],
                 num_points=self.num_points,
+                seed=self.seed,
             )
 
     def previous_step(self):
@@ -195,6 +205,25 @@ def previous_step(self):
                 m_y=self.m_y[self.current_step],
                 m_z=self.m_z[self.current_step],
                 num_points=self.num_points,
+                seed=self.seed,
+            )
+
+    def set_step(self, step: int):
+        """Set the detailed result to a specific step."""
+        if self.detailed_result is None:
+            return
+        if 0 <= step < len(self.m_y):
+            self.current_step = step
+            self.detailed_result = SectionDetailedResultState(
+                section=self.section,
+                eps_a=self.eps_axial[self.current_step],
+                chi_y=self.chi_y[self.current_step],
+                chi_z=self.chi_z[self.current_step],
+                n=self.n,
+                m_y=self.m_y[self.current_step],
+                m_z=self.m_z[self.current_step],
+                num_points=self.num_points,
+                seed=self.seed,
             )
 
 
@@ -202,8 +231,32 @@ class SectionDetailedResultState:
     """Something."""
 
     def __init__(
-        self, section, eps_a, chi_y, chi_z, n, m_y, m_z, num_points=1000
+        self,
+        section,
+        eps_a,
+        chi_y,
+        chi_z,
+        n,
+        m_y,
+        m_z,
+        num_points=1000,
+        seed=None,
     ):
+        """Create the SectionDetailedResult.
+
+        Arguments:
+            section: The section object.
+            eps_a: The axial strain at the section centroid.
+            chi_y: The curvature about the y-axis.
+            chi_z: The curvature about the z-axis.
+            n: The axial force.
+            m_y: The bending moment about the y-axis.
+            m_z: The bending moment about the z-axis.
+            num_points (int): Number of random points to sample in each surface
+                geometry (default = 1000).
+            seed (int): Seed for random number generator to ensure
+                reproducibility of random points.
+        """
         self.eps_a = eps_a
         self.chi_y = chi_y
         self.chi_z = chi_z
@@ -211,6 +264,7 @@ def __init__(
         self.m_y = m_y
         self.m_z = m_z
         self.section = section
+        self.seed = seed
         self._create_data_structure(num_points=num_points)
 
     def _create_data_structure(self, num_points=1000):
@@ -221,7 +275,9 @@ def _create_data_structure(self, num_points=1000):
         # SurfaceGeometry random points
         for surf in self.section.geometry.geometries:
             # Use surf.random_points_within to get random points
-            y, z = surf.random_points_within(num_points=num_points)
+            y, z = surf.random_points_within(
+                num_points=num_points, seed=self.seed
+            )
             strain = self.eps_a - self.chi_z * y + self.chi_y * z
             stress = surf.material.constitutive_law.get_stress(strain)
             surface_data.append(
diff --git a/structuralcodes/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py
index 29f9726b..6d96c261 100644
--- a/structuralcodes/geometry/_geometry.py
+++ b/structuralcodes/geometry/_geometry.py
@@ -455,11 +455,14 @@ def polygon(self) -> Polygon:
         """Returns the Shapely Polygon."""
         return self._polygon
 
-    def random_points_within(self, num_points: int = 100) -> np.ndarray:
+    def random_points_within(
+        self, num_points: int = 100, seed: int = None
+    ) -> np.ndarray:
         """Returns coordinates of random points within the polygon.
 
         Arguments:
             num_points (int): Number of random points to generate.
+            seed (int): Seed for the random number generator.
 
         Returns:
             x, y (ndarray, ndarray): Arrays with the x and y coordinates of the
@@ -506,7 +509,12 @@ def _prepare_triangulation_data(poly: Polygon) -> dict[str:ArrayLike]:
 
         xs = np.array([])
         ys = np.array([])
-        for tr in triangles['triangles']:
+        # Manage reproducibility of random points:
+        # seeds must be random in the triangles loop
+        n_tr = len(triangles['triangles'])
+        rng = np.random.default_rng(seed)
+        seeds = rng.integers(1, 300, n_tr)
+        for tr, s in zip(triangles['triangles'], seeds):
             # Get vertices for the triangle
             Ax = triangles['vertices'][tr[0]][0]
             Ay = triangles['vertices'][tr[0]][1]
@@ -522,8 +530,9 @@ def _prepare_triangulation_data(poly: Polygon) -> dict[str:ArrayLike]:
             # number of points in this triangle (at least 1)
             n = max(1, int(num_points * a / self.area))
             # generate random points in the triangle
-            r1 = np.random.uniform(0, 1, n)
-            r2 = np.random.uniform(0, 1, n)
+            rng = np.random.default_rng(s)
+            r1 = rng.uniform(0, 1, n)
+            r2 = rng.uniform(0, 1, n)
             x = (
                 (1 - np.sqrt(r1)) * Ax
                 + (np.sqrt(r1) * (1 - r2)) * Bx
@@ -910,6 +919,8 @@ def get_point_geometries(
                 coords.append([pg.x, pg.y])
                 materials.append(pg.material)
 
+        if len(coords) == 0:
+            return None, None
         return np.array(coords), materials
 
     def calculate_extents(self) -> t.Tuple[float, float, float, float]:
diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index 2d514015..f38c711b 100644
--- a/tests/test_geometry/test_geometry.py
+++ b/tests/test_geometry/test_geometry.py
@@ -826,20 +826,18 @@ def test_random_points_within_reproducibility():
     geo = RectangularGeometry(width, height, concrete, True)
 
     # Set seed for reproducibility
-    np.random.seed(42)
-    xs1, ys1 = geo.random_points_within(50)
+    xs1, ys1 = geo.random_points_within(50, seed=42)
 
     # Reset seed to same value
-    np.random.seed(42)
-    xs2, ys2 = geo.random_points_within(50)
+    xs2, ys2 = geo.random_points_within(50, seed=42)
 
     # Results should be identical
     np.testing.assert_array_equal(xs1, xs2)
     np.testing.assert_array_equal(ys1, ys2)
 
 
-def test_get_coordinates_point_geometries():
-    """Test get_coordinates_point_geometries method with filtering."""
+def test_get_point_geometries():
+    """Test get_point_geometries method with filtering."""
     # Create materials
     concrete = ConcreteMC2010(45)
     reinforcement = ReinforcementMC2010(450, 200000, 450, 0.075)
@@ -896,7 +894,7 @@ def test_get_coordinates_point_geometries():
     )
 
     # Test without group_label filter - should return all points
-    all_coords = geo.get_coordinates_point_geometries()
+    all_coords, _ = geo.get_point_geometries()
     total_reinforcement_points = sum(number)  # 12 regular reinforcement
     total_prestress_points = 2  # 2 prestress reinforcement
     expected_total = total_reinforcement_points + total_prestress_points
@@ -904,25 +902,22 @@ def test_get_coordinates_point_geometries():
     assert all_coords.shape[1] == 2  # x, y coordinates
 
     # Test with "reinforcement" filter
-    reinforcement_coords = geo.get_coordinates_point_geometries(
+    reinforcement_coords, _ = geo.get_point_geometries(
         group_label='reinforcement'
     )
     assert reinforcement_coords.shape[0] == total_reinforcement_points
     assert reinforcement_coords.shape[1] == 2
 
     # Test with "prestress_reinforcement" filter
-    prestress_coords = geo.get_coordinates_point_geometries(
+    prestress_coords, _ = geo.get_point_geometries(
         group_label='prestress_reinforcement'
     )
     assert prestress_coords.shape[0] == total_prestress_points
     assert prestress_coords.shape[1] == 2
 
     # Test with wrong filter - should return empty array
-    wrong_coords = geo.get_coordinates_point_geometries(
-        group_label='prestress'
-    )
-    assert wrong_coords.shape[0] == 0
-    assert wrong_coords.shape[1] == 2  # Still 2 columns even if empty
+    wrong_coords, _ = geo.get_point_geometries(group_label='prestress')
+    assert wrong_coords is None
 
     # Verify that the sum of filtered results equals total
     assert (