diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py
index 81da6d85..b7c9e84b 100644
--- a/structuralcodes/core/_section_results.py
+++ b/structuralcodes/core/_section_results.py
@@ -160,8 +160,186 @@ 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
+    seed: int = None
+    current_step: int = None
+    num_points: int = None
+
+    def create_detailed_result(self, num_points=1000):
+        """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],
+            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,
+            seed=self.seed,
+        )
+        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,
+                seed=self.seed,
+            )
+
+    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,
+                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,
+            )
+
+
+class SectionDetailedResultState:
+    """Something."""
+
+    def __init__(
+        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
+        self.n = n
+        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):
+        # 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=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(
+                {
+                    '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
@@ -178,6 +356,84 @@ 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
+
+    detailed_result: SectionDetailedResultState = 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_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_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/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py
index 6398b8a2..6d96c261 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,101 @@ def polygon(self) -> Polygon:
         """Returns the Shapely Polygon."""
         return self._polygon
 
+    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
+            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([])
+        # 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]
+            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 * Ax
+            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
+            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
+                + (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.
 
@@ -800,6 +896,33 @@ def area(self) -> float:
             area += geo.area
         return area
 
+    def get_point_geometries(
+        self, group_label: t.Optional[str] = None
+    ) -> 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:
+            Tuple[ndarray, List[Material]]: Coordinates array of shape (n, 2)
+            and list of materials.
+
+        """
+        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)
+
+        if len(coords) == 0:
+            return None, None
+        return np.array(coords), materials
+
     def calculate_extents(self) -> t.Tuple[float, float, float, float]:
         """Calculate extents of CompundGeometry.
 
diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py
index fd3ecc3c..13f9f221 100644
--- a/structuralcodes/sections/_generic.py
+++ b/structuralcodes/sections/_generic.py
@@ -784,6 +784,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
 
@@ -826,6 +827,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)
diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index 8f59ba92..f38c711b 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
 
@@ -28,7 +28,9 @@
     ElasticPlastic,
     ParabolaRectangle,
 )
-from structuralcodes.materials.reinforcement import ReinforcementMC2010
+from structuralcodes.materials.reinforcement import (
+    ReinforcementMC2010,
+)
 
 
 # Test create line
@@ -735,3 +737,190 @@ 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
+    xs1, ys1 = geo.random_points_within(50, seed=42)
+
+    # Reset seed to same value
+    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_point_geometries():
+    """Test get_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_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_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_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_point_geometries(group_label='prestress')
+    assert wrong_coords is None
+
+    # Verify that the sum of filtered results equals total
+    assert (
+        reinforcement_coords.shape[0] + prestress_coords.shape[0]
+        == all_coords.shape[0]
+    )