diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py index 81da6d85..784a08b8 100644 --- a/structuralcodes/core/_section_results.py +++ b/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 diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py index fd3ecc3c..b30a3923 100644 --- a/structuralcodes/sections/_generic.py +++ b/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 diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py index c2206059..6e12539e 100644 --- a/tests/test_sections/test_generic_section.py +++ b/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 + )