Get stress at point yz

structuralcodes/core/_section_results.py

@@ -2,6 +2,7 @@
 
 from __future__ import annotations  # To have clean hints of ArrayLike in docs
 
+import typing as t
 from dataclasses import dataclass, field, fields
 
 import numpy as np
@@ -264,3 +265,25 @@ class MMInteractionDomain(InteractionDomain):
 
     num_theta: float = 0  # number of discretizations along the angle
     theta: ArrayLike = None  # Array with shape (n,) containing the angle of NA
+
+
+@dataclass
+class PointStressResult:
+    """Class for storing stress results at a specific point in a section.
+
+    Attributes:
+        y (float): Y-coordinate of the point.
+        z (float): Z-coordinate of the point.
+        strain (float): Strain at the point.
+        stress (float): Stress at the point.
+        geometry_name (str, optional): Name of the geometry containing the
+            point.
+        geometry_group_label (str, optional): Group label of the geometry
+    """
+
+    y: float = 0.0
+    z: float = 0.0
+    strain: float = 0.0
+    stress: float = 0.0
+    geometry_name: t.Optional[str] = None
+    geometry_group_label: t.Optional[str] = None

structuralcodes/sections/_generic.py

@@ -10,6 +10,7 @@
 from numpy.typing import ArrayLike, NDArray
 from scipy.linalg import lu_factor, lu_solve
 from shapely import MultiPolygon
+from shapely.geometry import Point
 from shapely.ops import unary_union
 
 import structuralcodes.core._section_results as s_res
@@ -1493,3 +1494,82 @@ def calculate_strain_profile(
             raise StopIteration('Maximum number of iterations reached.')
 
         return strain.tolist()
+
+    def get_stress_point(
+        self,
+        y: float,
+        z: float,
+        eps_a: float,
+        chi_y: float,
+        chi_z: float,
+        where: t.Literal[0, 1, 2] = 0,
+    ) -> t.List[s_res.PointStressResult]:
+        """Get the stress at a given point (y, z) inside the cross-section.
+
+        Evaluates all geometries that contain the point and returns a list
+        of results. The `where` parameter filters which geometry types to
+        evaluate.
+
+        Args:
+            y (float): Y-coordinate of the point.
+            z (float): Z-coordinate of the point.
+            eps_a (float): Strain at (0,0).
+            chi_y (float): Curvature around the Y-axis (1/mm).
+            chi_z (float): Curvature around the Z-axis (1/mm).
+            where (int): 0: all geometries; 1: surface geometries only;
+                2: point geometries only.
+
+        Returns:
+            List[PointStressResult]: List of stress results at the point
+                (y, z). Returns an empty list if the point is not found
+                within any geometry.
+        """
+        if y is None or z is None:
+            return []
+
+        geom = self.section.geometry
+        pt = Point(y, z)
+        results = []
+
+        # Calculate strain at the point
+        # According to the codebase convention:
+        # eps = eps_a + chi_y * z - chi_z * y
+        strain = eps_a + z * chi_y - y * chi_z
+
+        # Check point geometries (reinforcement)
+        if where in [0, 2]:
+            for g in geom.point_geometries:
+                center = g._point
+                r = g._diameter / 2
+                poly = center.buffer(r)
+                if poly.contains(pt) or poly.touches(pt):
+                    stress = g.material.constitutive_law.get_stress(strain)
+                    results.append(
+                        s_res.PointStressResult(
+                            y=y,
+                            z=z,
+                            strain=strain,
+                            stress=stress,
+                            geometry_name=g.name,
+                            geometry_group_label=g.group_label,
+                        )
+                    )
+
+        # Check surface geometries
+        if where in [0, 1]:
+            for g in geom.geometries:
+                poly = g.polygon
+                if poly.contains(pt) or poly.touches(pt):
+                    stress = g.material.constitutive_law.get_stress(strain)
+                    results.append(
+                        s_res.PointStressResult(
+                            y=y,
+                            z=z,
+                            strain=strain,
+                            stress=stress,
+                            geometry_name=g.name,
+                            geometry_group_label=g.group_label,
+                        )
+                    )
+
+        return results

tests/test_sections/test_generic_section.py

@@ -6,6 +6,7 @@
 import pytest
 from shapely import Polygon
 
+from structuralcodes import set_design_code
 from structuralcodes.codes.ec2_2004 import reinforcement_duct_props
 from structuralcodes.geometry import (
     CircularGeometry,
@@ -20,7 +21,11 @@
     ElasticPlasticMaterial,
     GenericMaterial,
 )
-from structuralcodes.materials.concrete import ConcreteEC2_2004, ConcreteMC2010
+from structuralcodes.materials.concrete import (
+    ConcreteEC2_2004,
+    ConcreteMC2010,
+    create_concrete,
+)
 from structuralcodes.materials.constitutive_laws import (
     Elastic,
     ElasticPlastic,
@@ -33,6 +38,7 @@
 from structuralcodes.materials.reinforcement import (
     ReinforcementEC2_2004,
     ReinforcementMC2010,
+    create_reinforcement,
 )
 from structuralcodes.sections import (
     GenericSection,
@@ -1658,3 +1664,112 @@ def test_issue_cracked_properties():
     assert cracked_properties_before.isclose(
         cracked_properties_after, rtol=1e-3, atol=1e-6
     )
+
+
+def test_get_stress_point():
+    """Test get_stress_point method using validated results from notebook."""
+    # Set the active design code
+    set_design_code('ec2_2004')
+    # Create materials
+    concrete = create_concrete(fck=45)
+    reinforcement = create_reinforcement(
+        fyk=500, Es=200000, ftk=550, epsuk=0.07
+    )
+
+    # Create rectangular geometry
+    width = 250
+    height = 500
+    polygon = Polygon(
+        [
+            (-width / 2, -height / 2),
+            (width / 2, -height / 2),
+            (width / 2, height / 2),
+            (-width / 2, height / 2),
+        ]
+    )
+    geometry = SurfaceGeometry(
+        poly=polygon, material=concrete, group_label='concrete'
+    )
+
+    # Add reinforcement
+    diameter_reinf = 25
+    cover = 50
+    geometry = add_reinforcement(
+        geometry,
+        (
+            -width / 2 + cover + diameter_reinf / 2,
+            -height / 2 + cover + diameter_reinf / 2,
+        ),
+        diameter_reinf,
+        reinforcement,
+        group_label='bar 1',
+    )
+    geometry = add_reinforcement(
+        geometry,
+        (
+            width / 2 - cover - diameter_reinf / 2,
+            -height / 2 + cover + diameter_reinf / 2,
+        ),
+        diameter_reinf,
+        reinforcement,
+        group_label='bar 2',
+    )
+
+    # Create section
+    section = GenericSection(geometry)
+
+    # Calculate strain profile
+    strain = section.section_calculator.calculate_strain_profile(
+        n=-200.0 * 1e3,
+        my=-200.0 * 1e6,
+        mz=0.0 * 1e6,
+    )
+
+    # Test point at reinforcement location (validated results)
+    y_bar = width / 2 - cover - diameter_reinf / 2  # 62.5
+    z_bar = -height / 2 + cover + diameter_reinf / 2  # -187.5
+
+    results = section.section_calculator.get_stress_point(
+        y=y_bar,
+        z=z_bar,
+        eps_a=strain[0],
+        chi_y=strain[1],
+        chi_z=strain[2],
+        where=0,
+    )
+
+    # Find bar and concrete results
+    bar_result = next(
+        (r for r in results if r.geometry_group_label == 'bar 2'), None
+    )
+    concrete_result = next(
+        (r for r in results if r.geometry_group_label == 'concrete'), None
+    )
+
+    # Validate against notebook results
+    assert bar_result is not None
+    assert math.isclose(bar_result.stress, 434.4545058908614, rel_tol=1e-3)
+    assert math.isclose(bar_result.strain, 0.002172272529454307, rel_tol=1e-6)
+
+    assert concrete_result is not None
+    assert math.isclose(concrete_result.stress, 0.0, abs_tol=1e-6)
+    assert math.isclose(
+        concrete_result.strain, 0.002172272529454307, rel_tol=1e-6
+    )
+
+    # Test point at top of section (validated results)
+    results_top = section.section_calculator.get_stress_point(
+        y=0,
+        z=height / 2,  # 250.0
+        eps_a=strain[0],
+        chi_y=strain[1],
+        chi_z=strain[2],
+        where=0,
+    )
+
+    concrete_top = results_top[0]
+    assert concrete_top.geometry_group_label == 'concrete'
+    assert math.isclose(concrete_top.stress, -26.341523317880416, rel_tol=1e-3)
+    assert math.isclose(
+        concrete_top.strain, -0.0013015754221469022, rel_tol=1e-6
+    )