diff --git a/structuralcodes/core/_section_results.py b/structuralcodes/core/_section_results.py index 60cf12c8..3dbde253 100644 --- a/structuralcodes/core/_section_results.py +++ b/structuralcodes/core/_section_results.py @@ -978,6 +978,154 @@ def get_point_stress( ) +@dataclass(slots=True) +class IntegrateStrainStiffnessResult: + """Class for storing the results from integrating modulus.""" + + # Input strain profile + 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 + + # Tangent + tangent: NDArray[np.float64] = field( + default_factory=lambda: np.zeros((3, 3), dtype=float) + ) + + def asarray(self, dtype=np.float64) -> NDArray: + """Return an array representation of the tangent.""" + return np.asarray(self.tangent, dtype=dtype) + + +@dataclass(slots=True) +class IntegrateStrainForceResult: + """Class for storing the results from integrating stresses.""" + + # Input strain profile + 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 + + # Integrated loads + 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 + + # 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 + + def asarray(self, dtype=np.float64) -> NDArray: + """Return an array representation of the forces.""" + return np.array([self.n, self.m_y, self.m_z], dtype=dtype) + + def astuple(self) -> tuple[float, float, float]: + """Return a tuple representation of the forces.""" + return (self.n, self.m_y, self.m_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. diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py index 588d5d21..ccfe7ff9 100644 --- a/structuralcodes/sections/_generic.py +++ b/structuralcodes/sections/_generic.py @@ -7,7 +7,7 @@ import warnings import numpy as np -from numpy.typing import ArrayLike, NDArray +from numpy.typing import ArrayLike from scipy.linalg import lu_factor, lu_solve from shapely import MultiPolygon from shapely.ops import unary_union @@ -835,11 +835,23 @@ def _rotate_integration_data(self, theta: float): ) self.integration_data = rotated_integration_data + @t.overload + def integrate_strain_profile( + self, strain: ArrayLike, integrate: t.Literal['stress'] = 'stress' + ) -> s_res.IntegrateStrainForceResult: ... + + @t.overload + def integrate_strain_profile( + self, strain: ArrayLike, integrate: t.Literal['modulus'] + ) -> s_res.IntegrateStrainStiffnessResult: ... + def integrate_strain_profile( self, strain: ArrayLike, integrate: t.Literal['stress', 'modulus'] = 'stress', - ) -> t.Union[t.Tuple[float, float, float], NDArray]: + ) -> t.Union[ + s_res.IntegrateStrainForceResult, s_res.IntegrateStrainStiffnessResult + ]: """Integrate a strain profile returning stress resultants or tangent section stiffness matrix. @@ -854,16 +866,19 @@ def integrate_strain_profile( tangent section stiffness matrix (default is 'stress'). Returns: - Union(Tuple(float, float, float),NDArray): N, My and Mz when - `integrate='stress'`, or a numpy array representing the stiffness - matrix then `integrate='modulus'`. + Returns: + IntegrateStrainForcesResult + When ``integrate="stress"``. + + IntegrateStrainStiffnessResult + When ``integrate="modulus"``. Examples: result = self.integrate_strain_profile(strain,integrate='modulus') - # `result` will be the tangent stiffness matrix (a 3x3 numpy array) + # `result` will contain the tangent stiffness matrix (a 3x3 array) result = self.integrate_strain_profile(strain) - # `result` will be a tuple containing section forces (N, My, Mz) + # `result` will contain section forces (N, My, Mz) Raises: ValueError: If a unkown value is passed to the `integrate` @@ -888,10 +903,21 @@ def integrate_strain_profile( # b. (NDArray, integration_data) # (i.e. section stiff matrix, integration_data) # We need to return only forces or stifness - # (without triangultion data) - if len(result) == 2: - return result[0] - return result[:-1] + # (without triangulation data) + if integrate == 'stress': + return s_res.IntegrateStrainForceResult( + eps_a=strain[0], + chi_y=strain[1], + chi_z=strain[2], + n=result[0], + m_y=result[1], + m_z=result[2], + ) + if integrate == 'modulus': + return s_res.IntegrateStrainStiffnessResult(tangent=result[0]) + raise Exception( + 'integrate argument not valid. Valid options: stress and modulus' + ) def calculate_bending_strength( self, theta=0, n=0, max_iter: int = 100, tol: float = 1e-2 @@ -1773,7 +1799,7 @@ def calculate_strain_profile( break # Calculate response and residuals - response = np.array(self.integrate_strain_profile(strain=strain)) + response = self.integrate_strain_profile(strain=strain).asarray() residual = loads - response # Append to history variables @@ -1785,14 +1811,9 @@ def calculate_strain_profile( delta_strain = lu_solve((lu, piv), residual) else: # Calculate the current tangent stiffness - stiffness_tangent, _ = ( - self.integrator.integrate_strain_response_on_geometry( - geom, - strain, - integrate='modulus', - integration_data=self.integration_data, - ) - ) + stiffness_tangent = self.integrate_strain_profile( + strain=strain, integrate='modulus' + ).asarray() # Solve using the current tangent stiffness delta_strain = np.linalg.solve(stiffness_tangent, residual) diff --git a/tests/test_sections/test_generic_section.py b/tests/test_sections/test_generic_section.py index 405c65de..e862a08e 100644 --- a/tests/test_sections/test_generic_section.py +++ b/tests/test_sections/test_generic_section.py @@ -73,7 +73,7 @@ def test_rectangular_section(): # Use integrate_strain_response N, My, Mz = sec.section_calculator.integrate_strain_profile( (res_marin.eps_a, res_marin.chi_y, res_marin.chi_z) - ) + ).astuple() assert math.isclose(N, res_marin.n) assert math.isclose(My, res_marin.m_y) @@ -174,7 +174,7 @@ def test_rectangular_section_tangent_stiffness(b, h, E, integrator): # compute stiffness matrix stiffness = sec.section_calculator.integrate_strain_profile( [0, 0, 0], 'modulus' - ) + ).asarray() assert stiffness.shape == (3, 3) stiffness /= E @@ -256,7 +256,7 @@ def test_rectangular_rc_section_initial_tangent_stiffness( # compute initial stiffness matrix (gross) stiffness = sec.section_calculator.integrate_strain_profile( [0, 0, 0], 'modulus' - ) + ).asarray() assert stiffness.shape == (3, 3) assert math.isclose(EA_gross, stiffness[0, 0], rel_tol=1e-3) @@ -351,7 +351,7 @@ def test_rectangular_rc_section_tangent_stiffness( stiffness = sec.section_calculator.integrate_strain_profile( [eps_a, chi_y, 0], 'modulus' - ) + ).asarray() assert stiffness.shape == (3, 3) @@ -396,7 +396,7 @@ def test_rectangular_rc_section_tangent_stiffness( stiffness2 = sec.section_calculator.integrate_strain_profile( [0, 0, 0], 'modulus' - ) + ).asarray() assert np.allclose( stiffness, @@ -1505,7 +1505,9 @@ def test_section_parallel_material_elastic(): ) section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 1e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 1e-5, 0) + ).asarray() M = res[1] # The same with a GenericMaterial with two parallel materials @@ -1518,7 +1520,9 @@ def test_section_parallel_material_elastic(): ) section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 1e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 1e-5, 0) + ).asarray() M_p = res[1] # Check they are the same @@ -1534,7 +1538,9 @@ def test_section_parallel_material_elastic(): ) section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 1e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 1e-5, 0) + ).asarray() M_p = res[1] # Check they are the same @@ -1549,7 +1555,9 @@ def test_section_parallel_material_elasticplastic(): geo = RectangularGeometry(width=100, height=100, material=mat) section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 3e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 3e-5, 0) + ).asarray() M = res[1] # The same with a GenericMaterial with two parallel materials @@ -1560,7 +1568,9 @@ def test_section_parallel_material_elasticplastic(): geo = RectangularGeometry(width=100, height=100, material=mat) section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 3e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 3e-5, 0) + ).asarray() M_p = res[1] # Check they are the same @@ -1575,7 +1585,9 @@ def test_section_parallel_material_elasticplastic(): geo = RectangularGeometry(width=100, height=100, material=mat) section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 3e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 3e-5, 0) + ).asarray() M_p = res[1] # Check they are the same @@ -1597,13 +1609,17 @@ def test_section_parallel_marin_concrete_tension(): geometry=geo, integrator='fiber', mesh_size=0.0001 ) - res = section.section_calculator.integrate_strain_profile((0, 4e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 4e-5, 0) + ).asarray() M_f = res[1] # Use marin integrator section = GenericSection(geometry=geo) - res = section.section_calculator.integrate_strain_profile((0, 4e-5, 0)) + res = section.section_calculator.integrate_strain_profile( + (0, 4e-5, 0) + ).asarray() M_m = res[1] # Check they are the same