From ff8e6dbbba5502c57ef4ab768c7a97e97b48900c Mon Sep 17 00:00:00 2001 From: Morten Engen Date: Tue, 17 Mar 2026 22:53:24 +0100 Subject: [PATCH 1/3] Split base classes in separate modules Note that some methods on the SectionCalculator base class that were not general are removed. --- structuralcodes/core/base.py | 392 ------------------ structuralcodes/core/base/__init__.py | 14 + .../core/base/_constitutive_law.py | 168 ++++++++ structuralcodes/core/base/_geometry.py | 127 ++++++ structuralcodes/core/base/_material.py | 150 +++++++ structuralcodes/core/base/_section.py | 46 ++ 6 files changed, 505 insertions(+), 392 deletions(-) delete mode 100644 structuralcodes/core/base.py create mode 100644 structuralcodes/core/base/__init__.py create mode 100644 structuralcodes/core/base/_constitutive_law.py create mode 100644 structuralcodes/core/base/_geometry.py create mode 100644 structuralcodes/core/base/_material.py create mode 100644 structuralcodes/core/base/_section.py diff --git a/structuralcodes/core/base.py b/structuralcodes/core/base.py deleted file mode 100644 index 0940bc5f..00000000 --- a/structuralcodes/core/base.py +++ /dev/null @@ -1,392 +0,0 @@ -"""Abstract base classes.""" - -from __future__ import annotations # To have clean hints of ArrayLike in docs - -import abc -import typing as t - -import numpy as np -from numpy.typing import ArrayLike - -import structuralcodes.core._section_results as s_res - - -class Material(abc.ABC): - """Abstract base class for materials.""" - - _constitutive_law = None - _initial_strain: t.Optional[float] = None - _initial_stress: t.Optional[float] = None - _strain_compatibility: t.Optional[bool] = None - - def __init__( - self, - density: float, - initial_strain: t.Optional[float] = None, - initial_stress: t.Optional[float] = None, - strain_compatibility: t.Optional[bool] = None, - name: t.Optional[str] = None, - ) -> None: - """Initializes an instance of a new material. - - Args: - density (float): Density of the material in kg/m3. - - Keyword Args: - initial_strain (Optional[float]): Initial strain of the material. - initial_stress (Optional[float]): Initial stress of the material. - strain_compatibility (Optional[bool]): Only relevant if - initial_strain or initial_stress are different from zero. If - True, the material deforms with the geometry. If False, the - stress in the material upon loading is kept constant - corresponding to the initial strain. - name (Optional[str]): descriptive name of the material - - Raise: - ValueError: if both initial_strain and initial_stress are provided - """ - self._density = abs(density) - if initial_strain is not None and initial_stress is not None: - raise ValueError( - 'Both initial_strain and initial_stress cannot be provided.' - ) - self._initial_strain = initial_strain - self._initial_stress = initial_stress - self._strain_compatibility = strain_compatibility - self._name = name if name is not None else 'Material' - - @property - def constitutive_law(self) -> ConstitutiveLaw: - """Returns the ConstitutiveLaw of the object.""" - return self._constitutive_law - - @property - def name(self): - """Returns the name of the material.""" - return self._name - - @property - def density(self): - """Returns the density of the material in kg/m3.""" - return self._density - - @property - def initial_strain(self): - """Returns the initial strain of the material.""" - return self._initial_strain - - @property - def initial_stress(self): - """Returns the initial stress of the material.""" - return self._initial_stress - - @property - def strain_compatibility(self): - """Returns the strain compatibility of the material. - - If true (default), the strain compatibility is enforced - haveing the same strain as in all other materials of the - section at the same point. If false, the strain compatibility - is not enforced and the initial strain is applied to the section - independently. - """ - return self._strain_compatibility - - def _apply_initial_strain(self): - """Wraps the current constitutive law to apply initial strain.""" - strain_compatibility = ( - self._strain_compatibility - if self._strain_compatibility is not None - else True - ) - if self._initial_stress is not None: - # Specified a stress, compute the strain from it - self._initial_strain_from_stress() - if self._initial_strain is not None: - # Lazy import to avoid circular dependency - from structuralcodes.materials.constitutive_laws import ( # noqa: PLC0415 - InitialStrain, - ) - - if self._initial_stress is None: - # Compute the stress from the strain - self._initial_stress = self._constitutive_law.get_stress( - self._initial_strain - ) - - self._constitutive_law = InitialStrain( - self._constitutive_law, - self._initial_strain, - strain_compatibility, - ) - - def _initial_strain_from_stress(self): - """Computes the initial strain from the initial stress. - - This function is called internally so it assumes that the - initial stress is not None - """ - # Iteratively compute the initial strain that gives the desired - # initial stress. Note that the wrapped law can be nonlinear - tol = 1e-12 - max_iter = 100 - target_stress = self._initial_stress - strain = 0.0 - stress = self._constitutive_law.get_stress(strain) - d_stress = target_stress - stress - num_iter = 0 - while abs(d_stress) > tol and num_iter < max_iter: - tangent = self._constitutive_law.get_tangent(strain) - if tangent == 0: - raise ValueError( - 'Tangent modulus = 0 during initial strain computation.' - ) - d_strain = d_stress / tangent - strain += d_strain - stress = self._constitutive_law.get_stress(strain) - d_stress = target_stress - stress - num_iter += 1 - - if abs(d_stress) > tol: - raise RuntimeError('Failed to converge for given initial stress.') - - self._initial_strain = strain - - -class ConstitutiveLaw(abc.ABC): - """Abstract base class for constitutive laws.""" - - __materials__: t.Tuple[str] = () - constitutive_law_counter: t.ClassVar[int] = 0 - - def __init__(self, name: t.Optional[str] = None) -> None: - """Initialize a ConstitutiveLaw object.""" - self.id = self.constitutive_law_counter - self._name = name if name is not None else f'ConstitutiveLaw_{self.id}' - self._increase_global_counter() - - @property - def name(self): - """Returns the name of the constitutive law.""" - return self._name - - @classmethod - def _increase_global_counter(cls): - cls.constitutive_law_counter += 1 - - @abc.abstractmethod - def get_stress( - self, eps: t.Union[float, ArrayLike] - ) -> t.Union[float, ArrayLike]: - """Each constitutive law should provide a method to return the - stress given the strain level. - """ - - @abc.abstractmethod - def get_tangent( - self, eps: t.Union[float, ArrayLike] - ) -> t.Union[float, ArrayLike]: - """Each constitutive law should provide a method to return the - tangent at a given strain level. - """ - - @abc.abstractmethod - def get_ultimate_strain(self) -> t.Tuple[float, float]: - """Each constitutive law should provide a method to return the - ultimate strain (negative and positive). - """ - - def preprocess_strains_with_limits( - self, eps: t.Union[float, ArrayLike] - ) -> t.Union[float, ArrayLike]: - """Preprocess strain arrays setting those strains sufficiently - near to ultimate strain limits to exactly ultimate strain limit. - """ - eps = eps if np.isscalar(eps) else np.atleast_1d(eps) - eps_min, eps_max = self.get_ultimate_strain() - - if np.isscalar(eps): - if np.isclose(eps, eps_max, atol=1e-6): - return eps_max - if np.isclose(eps, eps_min, atol=1e-6): - return eps_min - return eps - idxs = np.isclose(eps, np.zeros_like(eps) + eps_max, atol=1e-6) - eps[idxs] = eps_max - idxs = np.isclose(eps, np.zeros_like(eps) + eps_min, atol=1e-6) - eps[idxs] = eps_min - - return eps - - def _discretize_law(self) -> ConstitutiveLaw: - """Discretize the law as a piecewise linear function.""" - - # Discretize the constitutive law in a "smart way" - def find_x_lim(x, y): - # Check if there are non-zero values for x > 0 - if np.any(y[0:] != 0): - # Find the last non-zero index for x > 0 - non_zero_indices = np.nonzero(y[0:])[0] - x_lim_index = 0 + non_zero_indices[-1] - return x[x_lim_index] - # All values are zero for x > 0 - return None - - eps_min, eps_max = self.get_ultimate_strain() - eps_max = min(eps_max, 1) - # Analise positive branch - eps = np.linspace(0, eps_max, 10000) - sig = self.get_stress(eps) - sig[(sig < np.max(sig) * 1e-6)] = 0 - eps_lim = find_x_lim(eps, sig) - # Now discretize the function in 10 steps for positive part - eps_pos = ( - np.linspace(0, -eps_min, 1) - if eps_lim is None - else np.linspace(0, eps_lim, 10) - ) - # Analise negative branch - eps = np.linspace(0, eps_min, 10000) - sig = -self.get_stress(eps) - sig[(sig < np.max(sig) * 1e-6)] = 0 - eps_lim = find_x_lim(-eps, sig) - # Now discretize the function in 10 steps for negative part - eps_neg = ( - np.linspace(eps_min, 0, 1, endpoint=False) - if eps_lim is None - else np.linspace(-eps_lim, 0, 10, endpoint=False) - ) - - eps = np.concatenate((eps_neg, eps_pos)) - sig = self.get_stress(eps) - from structuralcodes.materials.constitutive_laws import ( # noqa: PLC0415 - UserDefined, - ) - - return UserDefined(eps, sig) - - def __marin__(self, **kwargs): - """Function for getting the strain limits and coefficients - for marin integration. - - By default the law is discretized as a piecewise linear - function. Then marin coefficients are computed based on this - discretization. - """ - piecewise_law = self._discretize_law() - return piecewise_law.__marin__(**kwargs) - - def __marin_tangent__(self, **kwargs): - """Function for getting the strain limits and coefficients - for marin integration of tangent modulus. - - By default the law is discretized as a piecewise linear - function. Then marin coefficients are computed based on this - discretization. - """ - piecewise_law = self._discretize_law() - return piecewise_law.__marin_tangent__(**kwargs) - - def get_secant( - self, eps: t.Union[float, ArrayLike] - ) -> t.Union[float, ArrayLike]: - """Method to return the secant at a given strain level.""" - # Adjust eps if it is not scalar - eps = eps if np.isscalar(eps) else np.atleast_1d(eps) - - # Calculate secant for scalar eps - if np.isscalar(eps): - if eps != 0: - sig = self.get_stress(eps) - return sig / eps - return self.get_tangent(eps) - - # Calculate secant for array eps - secant = np.zeros_like(eps) - strain_is_zero = eps == 0 - strain_is_nonzero = eps != 0 - secant[strain_is_zero] = self.get_tangent(eps[strain_is_zero]) - secant[strain_is_nonzero] = ( - self.get_stress(eps[strain_is_nonzero]) / eps[strain_is_nonzero] - ) - return secant - - -class Section(abc.ABC): - """Abstract base class for a cross secion. - The section is defined by local axes y and z (mapped to x and y cartesian - plane respectively). - """ - - section_counter: t.ClassVar[int] = 0 - - def __init__(self, name: t.Optional[str] = None) -> None: - """Initialize a Section object.""" - self.id = self.section_counter - self._name = name if name is not None else f'Section_{self.id}' - self._increase_global_counter() - - @property - def name(self): - """Returns the name of the section.""" - return self._name - - @classmethod - def _increase_global_counter(cls): - cls.section_counter += 1 - - -class SectionCalculator(abc.ABC): - """Abstract class for SectionCalculators - defining the interface. - """ - - def __init__(self, section: Section) -> None: - """Initialization of SectionCalculator object, providing - a Section object. - """ - self.section = section - - @abc.abstractmethod - def _calculate_gross_section_properties(self) -> s_res.SectionProperties: - """Calculates the gross section properties of the section - This function is private and called when the section is created - It stores the result into the result object. - - Returns: - gross_section_properties (GrossSection) - """ - - @abc.abstractmethod - def calculate_bending_strength( - self, theta=0, n=0 - ) -> s_res.UltimateBendingMomentResults: - """Calculates the bending strength for given inclination of n.a. - and axial load. - - Input: - theta (float, default = 0): inclination of n.a. respect to y axis - n (float, default = 0): axial load applied to the section (+: tension, - -: compression) - - Return: - ultimate_bending_moment_result (UltimateBendingMomentResult) - """ - - @abc.abstractmethod - def calculate_moment_curvature( - self, theta=0, n=0 - ) -> s_res.MomentCurvatureResults: - """Calculates the moment-curvature relation for given inclination of - n.a. and axial load. - - Input: - theta (float, default = 0): inclination of n.a. respect to y axis - n (float, default = 0): axial load applied to the section - (+: tension, -: compression) - chi_incr (float, default = 1e-8): the curvature increment for the - analysis - - Return: - moment_curvature_result (MomentCurvatureResults) - """ diff --git a/structuralcodes/core/base/__init__.py b/structuralcodes/core/base/__init__.py new file mode 100644 index 00000000..54f8547c --- /dev/null +++ b/structuralcodes/core/base/__init__.py @@ -0,0 +1,14 @@ +"""Abstract base classes for the class hierarchy.""" + +from ._constitutive_law import ConstitutiveLaw +from ._geometry import Geometry +from ._material import Material +from ._section import Section, SectionCalculator + +__all__ = [ + 'ConstitutiveLaw', + 'Geometry', + 'Material', + 'Section', + 'SectionCalculator', +] diff --git a/structuralcodes/core/base/_constitutive_law.py b/structuralcodes/core/base/_constitutive_law.py new file mode 100644 index 00000000..badca37e --- /dev/null +++ b/structuralcodes/core/base/_constitutive_law.py @@ -0,0 +1,168 @@ +"""Abstract base class for constitutive laws.""" + +from __future__ import annotations # To have clean hints of ArrayLike in docs + +import abc +import typing as t + +import numpy as np +from numpy.typing import ArrayLike + + +class ConstitutiveLaw(abc.ABC): + """Abstract base class for constitutive laws.""" + + __materials__: t.Tuple[str] = () + constitutive_law_counter: t.ClassVar[int] = 0 + + def __init__(self, name: t.Optional[str] = None) -> None: + """Initialize a ConstitutiveLaw object.""" + self.id = self.constitutive_law_counter + self._name = name if name is not None else f'ConstitutiveLaw_{self.id}' + self._increase_global_counter() + + @property + def name(self): + """Returns the name of the constitutive law.""" + return self._name + + @classmethod + def _increase_global_counter(cls): + cls.constitutive_law_counter += 1 + + @abc.abstractmethod + def get_stress( + self, eps: t.Union[float, ArrayLike] + ) -> t.Union[float, ArrayLike]: + """Each constitutive law should provide a method to return the + stress given the strain level. + """ + + @abc.abstractmethod + def get_tangent( + self, eps: t.Union[float, ArrayLike] + ) -> t.Union[float, ArrayLike]: + """Each constitutive law should provide a method to return the + tangent at a given strain level. + """ + + @abc.abstractmethod + def get_ultimate_strain(self) -> t.Tuple[float, float]: + """Each constitutive law should provide a method to return the + ultimate strain (negative and positive). + """ + + def preprocess_strains_with_limits( + self, eps: t.Union[float, ArrayLike] + ) -> t.Union[float, ArrayLike]: + """Preprocess strain arrays setting those strains sufficiently + near to ultimate strain limits to exactly ultimate strain limit. + """ + eps = eps if np.isscalar(eps) else np.atleast_1d(eps) + eps_min, eps_max = self.get_ultimate_strain() + + if np.isscalar(eps): + if np.isclose(eps, eps_max, atol=1e-6): + return eps_max + if np.isclose(eps, eps_min, atol=1e-6): + return eps_min + return eps + idxs = np.isclose(eps, np.zeros_like(eps) + eps_max, atol=1e-6) + eps[idxs] = eps_max + idxs = np.isclose(eps, np.zeros_like(eps) + eps_min, atol=1e-6) + eps[idxs] = eps_min + + return eps + + def _discretize_law(self) -> ConstitutiveLaw: + """Discretize the law as a piecewise linear function.""" + + # Discretize the constitutive law in a "smart way" + def find_x_lim(x, y): + # Check if there are non-zero values for x > 0 + if np.any(y[0:] != 0): + # Find the last non-zero index for x > 0 + non_zero_indices = np.nonzero(y[0:])[0] + x_lim_index = 0 + non_zero_indices[-1] + return x[x_lim_index] + # All values are zero for x > 0 + return None + + eps_min, eps_max = self.get_ultimate_strain() + eps_max = min(eps_max, 1) + # Analise positive branch + eps = np.linspace(0, eps_max, 10000) + sig = self.get_stress(eps) + sig[(sig < np.max(sig) * 1e-6)] = 0 + eps_lim = find_x_lim(eps, sig) + # Now discretize the function in 10 steps for positive part + eps_pos = ( + np.linspace(0, -eps_min, 1) + if eps_lim is None + else np.linspace(0, eps_lim, 10) + ) + # Analise negative branch + eps = np.linspace(0, eps_min, 10000) + sig = -self.get_stress(eps) + sig[(sig < np.max(sig) * 1e-6)] = 0 + eps_lim = find_x_lim(-eps, sig) + # Now discretize the function in 10 steps for negative part + eps_neg = ( + np.linspace(eps_min, 0, 1, endpoint=False) + if eps_lim is None + else np.linspace(-eps_lim, 0, 10, endpoint=False) + ) + + eps = np.concatenate((eps_neg, eps_pos)) + sig = self.get_stress(eps) + from structuralcodes.materials.constitutive_laws import ( # noqa: PLC0415 + UserDefined, + ) + + return UserDefined(eps, sig) + + def __marin__(self, **kwargs): + """Function for getting the strain limits and coefficients + for marin integration. + + By default the law is discretized as a piecewise linear + function. Then marin coefficients are computed based on this + discretization. + """ + piecewise_law = self._discretize_law() + return piecewise_law.__marin__(**kwargs) + + def __marin_tangent__(self, **kwargs): + """Function for getting the strain limits and coefficients + for marin integration of tangent modulus. + + By default the law is discretized as a piecewise linear + function. Then marin coefficients are computed based on this + discretization. + """ + piecewise_law = self._discretize_law() + return piecewise_law.__marin_tangent__(**kwargs) + + def get_secant( + self, eps: t.Union[float, ArrayLike] + ) -> t.Union[float, ArrayLike]: + """Method to return the secant at a given strain level.""" + # Adjust eps if it is not scalar + eps = eps if np.isscalar(eps) else np.atleast_1d(eps) + + # Calculate secant for scalar eps + if np.isscalar(eps): + if eps != 0: + sig = self.get_stress(eps) + return sig / eps + return self.get_tangent(eps) + + # Calculate secant for array eps + secant = np.zeros_like(eps) + strain_is_zero = eps == 0 + strain_is_nonzero = eps != 0 + secant[strain_is_zero] = self.get_tangent(eps[strain_is_zero]) + secant[strain_is_nonzero] = ( + self.get_stress(eps[strain_is_nonzero]) / eps[strain_is_nonzero] + ) + return secant diff --git a/structuralcodes/core/base/_geometry.py b/structuralcodes/core/base/_geometry.py new file mode 100644 index 00000000..0ee62348 --- /dev/null +++ b/structuralcodes/core/base/_geometry.py @@ -0,0 +1,127 @@ +"""Abstract base class for geometries.""" + +from __future__ import annotations # To have clean hints of ArrayLike in docs + +import typing as t +from fnmatch import fnmatchcase + +from ._material import Material + + +class Geometry: + """Base class for a geometry object.""" + + section_counter: t.ClassVar[int] = 0 + + def __init__( + self, name: t.Optional[str] = None, group_label: t.Optional[str] = None + ) -> None: + """Initializes a geometry object. + + The name and grouplabel serve for filtering in a compound object. By + default it creates a new name each time. + + Arguments: + name (Optional(str)): The name to be given to the object. + group_label (Optional(str)): A label for grouping several objects. + """ + if name is not None: + self._name = name + else: + counter = Geometry.return_global_counter_and_increase() + self._name = f'Geometry_{counter}' + self._group_label = group_label + + @property + def name(self): + """Returns the name of the Geometry.""" + return self._name + + @property + def group_label(self): + """Returns the group_label fo the Geometry.""" + return self._group_label + + @classmethod + def _increase_global_counter(cls): + """Increases the global counter by one.""" + cls.section_counter += 1 + + @classmethod + def return_global_counter_and_increase(cls): + """Returns the current counter and increases it by one.""" + counter = cls.section_counter + cls._increase_global_counter() + return counter + + @staticmethod + def from_geometry( + geo: Geometry, + new_material: t.Optional[Material] = None, + ) -> Geometry: + """Create a new geometry with a different material.""" + raise NotImplementedError( + 'This method should be implemented by subclasses' + ) + + def _name_matches( + self, pattern: str, *, case_sensitive: bool = True + ) -> bool: + """Checks if the name matches a pattern. + + Arguments: + pattern (str): the string pattern to be checked + + Keyword Arguments: + case_sensitive (bool, optional): if True (default) the check is + case sensitive. + + Returns: + (bool): Returns True if the name matches the pattern. + + Note: + The matching permits to use: + - "*" any chars + - "?" single char + - "[abc]" character set + + Examples: + >>> geo.name_matches("nametos*") + >>> geo.name_matches("*pier*") + >>> geo.name_matches("Abutment??", case_senstive=False) + """ + if not case_sensitive: + return fnmatchcase(self.name.casefold(), pattern.casefold()) + return fnmatchcase(self.name, pattern) + + def _group_matches( + self, pattern: str, *, case_sensitive: bool = True + ) -> bool: + """Checks if the group_label matches a pattern. + + Arguments: + pattern (str): the string pattern to be checked + + Keyword Arguments: + case_sensitive (bool, optional): if True (default) the check is + case sensitive. + + Returns: + (bool): Returns True if the group_label matches the pattern. + + Note: + The matching permits to use: + - "*" any chars + - "?" single char + - "[abc]" character set + + Examples: + >>> geo.group_matches("nametos*") + >>> geo.group_matches("*pier*") + >>> geo.group_matches("Abutment??", case_senstive=False) + """ + if self.group_label is None: + return False + if not case_sensitive: + return fnmatchcase(self.group_label.casefold(), pattern.casefold()) + return fnmatchcase(self.group_label, pattern) diff --git a/structuralcodes/core/base/_material.py b/structuralcodes/core/base/_material.py new file mode 100644 index 00000000..27dd7a23 --- /dev/null +++ b/structuralcodes/core/base/_material.py @@ -0,0 +1,150 @@ +"""Abstract base class for materials.""" + +from __future__ import annotations # To have clean hints of ArrayLike in docs + +import abc +import typing as t + +from ._constitutive_law import ConstitutiveLaw + + +class Material(abc.ABC): + """Abstract base class for materials.""" + + _constitutive_law = None + _initial_strain: t.Optional[float] = None + _initial_stress: t.Optional[float] = None + _strain_compatibility: t.Optional[bool] = None + + def __init__( + self, + density: float, + initial_strain: t.Optional[float] = None, + initial_stress: t.Optional[float] = None, + strain_compatibility: t.Optional[bool] = None, + name: t.Optional[str] = None, + ) -> None: + """Initializes an instance of a new material. + + Args: + density (float): Density of the material in kg/m3. + + Keyword Args: + initial_strain (Optional[float]): Initial strain of the material. + initial_stress (Optional[float]): Initial stress of the material. + strain_compatibility (Optional[bool]): Only relevant if + initial_strain or initial_stress are different from zero. If + True, the material deforms with the geometry. If False, the + stress in the material upon loading is kept constant + corresponding to the initial strain. + name (Optional[str]): descriptive name of the material + + Raise: + ValueError: if both initial_strain and initial_stress are provided + """ + self._density = abs(density) + if initial_strain is not None and initial_stress is not None: + raise ValueError( + 'Both initial_strain and initial_stress cannot be provided.' + ) + self._initial_strain = initial_strain + self._initial_stress = initial_stress + self._strain_compatibility = strain_compatibility + self._name = name if name is not None else 'Material' + + @property + def constitutive_law(self) -> ConstitutiveLaw: + """Returns the ConstitutiveLaw of the object.""" + return self._constitutive_law + + @property + def name(self): + """Returns the name of the material.""" + return self._name + + @property + def density(self): + """Returns the density of the material in kg/m3.""" + return self._density + + @property + def initial_strain(self): + """Returns the initial strain of the material.""" + return self._initial_strain + + @property + def initial_stress(self): + """Returns the initial stress of the material.""" + return self._initial_stress + + @property + def strain_compatibility(self): + """Returns the strain compatibility of the material. + + If true (default), the strain compatibility is enforced + haveing the same strain as in all other materials of the + section at the same point. If false, the strain compatibility + is not enforced and the initial strain is applied to the section + independently. + """ + return self._strain_compatibility + + def _apply_initial_strain(self): + """Wraps the current constitutive law to apply initial strain.""" + strain_compatibility = ( + self._strain_compatibility + if self._strain_compatibility is not None + else True + ) + if self._initial_stress is not None: + # Specified a stress, compute the strain from it + self._initial_strain_from_stress() + if self._initial_strain is not None: + # Lazy import to avoid circular dependency + from structuralcodes.materials.constitutive_laws import ( # noqa: PLC0415 + InitialStrain, + ) + + if self._initial_stress is None: + # Compute the stress from the strain + self._initial_stress = self._constitutive_law.get_stress( + self._initial_strain + ) + + self._constitutive_law = InitialStrain( + self._constitutive_law, + self._initial_strain, + strain_compatibility, + ) + + def _initial_strain_from_stress(self): + """Computes the initial strain from the initial stress. + + This function is called internally so it assumes that the + initial stress is not None + """ + # Iteratively compute the initial strain that gives the desired + # initial stress. Note that the wrapped law can be nonlinear + tol = 1e-12 + max_iter = 100 + target_stress = self._initial_stress + strain = 0.0 + stress = self._constitutive_law.get_stress(strain) + d_stress = target_stress - stress + num_iter = 0 + while abs(d_stress) > tol and num_iter < max_iter: + tangent = self._constitutive_law.get_tangent(strain) + if tangent == 0: + raise ValueError( + 'Tangent modulus = 0 during initial strain computation.' + ) + d_strain = d_stress / tangent + strain += d_strain + stress = self._constitutive_law.get_stress(strain) + d_stress = target_stress - stress + num_iter += 1 + + if abs(d_stress) > tol: + raise RuntimeError('Failed to converge for given initial stress.') + + self._initial_strain = strain diff --git a/structuralcodes/core/base/_section.py b/structuralcodes/core/base/_section.py new file mode 100644 index 00000000..f57aae69 --- /dev/null +++ b/structuralcodes/core/base/_section.py @@ -0,0 +1,46 @@ +"""Abstract base classes for sections.""" + +from __future__ import annotations # To have clean hints of ArrayLike in docs + +import abc +import typing as t + +from ._geometry import Geometry + + +class Section(abc.ABC): + """Abstract base class for a cross secion. + The section is defined by local axes y and z (mapped to x and y cartesian + plane respectively). + """ + + geometry: Geometry + section_counter: t.ClassVar[int] = 0 + section_calculator: SectionCalculator + + def __init__(self, name: t.Optional[str] = None) -> None: + """Initialize a Section object.""" + self.id = self.section_counter + self._name = name if name is not None else f'Section_{self.id}' + self._increase_global_counter() + + @property + def name(self): + """Returns the name of the section.""" + return self._name + + @classmethod + def _increase_global_counter(cls): + cls.section_counter += 1 + + +class SectionCalculator(abc.ABC): + """Abstract class for SectionCalculators + defining the interface. + """ + + def __init__(self, section: Section) -> None: + """Initialization of SectionCalculator object, providing + a Section object. + """ + self.section = section From cd0fca80599458dcca3ad005fa2013f2dd9be5ee Mon Sep 17 00:00:00 2001 From: Morten Engen Date: Tue, 17 Mar 2026 22:55:11 +0100 Subject: [PATCH 2/3] Update imports of base classes --- structuralcodes/geometry/_circular.py | 3 +- structuralcodes/geometry/_geometry.py | 167 ++++-------------- structuralcodes/geometry/_rectangular.py | 3 +- structuralcodes/geometry/_reinforcement.py | 3 +- .../geometry/profiles/_base_profile.py | 2 +- .../materials/concrete/_concrete.py | 2 +- .../materials/reinforcement/_reinforcement.py | 2 +- structuralcodes/sections/_generic.py | 4 +- .../section_integrators/_fiber_integrator.py | 2 +- .../section_integrators/_marin_integrator.py | 2 +- 10 files changed, 45 insertions(+), 145 deletions(-) diff --git a/structuralcodes/geometry/_circular.py b/structuralcodes/geometry/_circular.py index b5c69e44..7114d47b 100644 --- a/structuralcodes/geometry/_circular.py +++ b/structuralcodes/geometry/_circular.py @@ -11,8 +11,7 @@ from numpy.typing import ArrayLike from shapely import Point -from structuralcodes.core.base import Material - +from ..core.base import Material from ._geometry import SurfaceGeometry diff --git a/structuralcodes/geometry/_geometry.py b/structuralcodes/geometry/_geometry.py index 054566d9..900700fa 100644 --- a/structuralcodes/geometry/_geometry.py +++ b/structuralcodes/geometry/_geometry.py @@ -4,7 +4,6 @@ import typing as t import warnings -from fnmatch import fnmatchcase from math import atan2 import numpy as np @@ -21,10 +20,11 @@ ) from shapely.ops import split -from structuralcodes.core.base import Material from structuralcodes.materials.basic import ElasticMaterial from structuralcodes.materials.concrete import Concrete +from ..core.base import Geometry, Material + def _mirror_about_axis_matrix(axis: LineString) -> np.ndarray: if not isinstance(axis, LineString): @@ -59,136 +59,6 @@ def _mirror_about_axis_matrix(axis: LineString) -> np.ndarray: return T_inv @ R_inv @ M @ R @ T -class Geometry: - """Base class for a geometry object.""" - - section_counter: t.ClassVar[int] = 0 - - def __init__( - self, name: t.Optional[str] = None, group_label: t.Optional[str] = None - ) -> None: - """Initializes a geometry object. - - The name and grouplabel serve for filtering in a compound object. By - default it creates a new name each time. - - Arguments: - name (Optional(str)): The name to be given to the object. - group_label (Optional(str)): A label for grouping several objects. - """ - if name is not None: - self._name = name - else: - counter = Geometry.return_global_counter_and_increase() - self._name = f'Geometry_{counter}' - self._group_label = group_label - - @property - def name(self): - """Returns the name of the Geometry.""" - return self._name - - @property - def group_label(self): - """Returns the group_label fo the Geometry.""" - return self._group_label - - @classmethod - def _increase_global_counter(cls): - """Increases the global counter by one.""" - cls.section_counter += 1 - - @classmethod - def return_global_counter_and_increase(cls): - """Returns the current counter and increases it by one.""" - counter = cls.section_counter - cls._increase_global_counter() - return counter - - @staticmethod - def from_geometry( - geo: Geometry, - new_material: t.Optional[Material] = None, - ) -> Geometry: - """Create a new geometry with a different material.""" - raise NotImplementedError( - 'This method should be implemented by subclasses' - ) - - def __add__(self, other: Geometry) -> CompoundGeometry: - """Add operator "+" for geometries. - - Arguments: - other (Geometry): The other geometry to add. - - Returns: - CompoundGeometry: A new CompoundGeometry. - """ - return CompoundGeometry([self, other]) - - def _name_matches( - self, pattern: str, *, case_sensitive: bool = True - ) -> bool: - """Checks if the name matches a pattern. - - Arguments: - pattern (str): the string pattern to be checked - - Keyword Arguments: - case_sensitive (bool, optional): if True (default) the check is - case sensitive. - - Returns: - (bool): Returns True if the name matches the pattern. - - Note: - The matching permits to use: - - "*" any chars - - "?" single char - - "[abc]" character set - - Examples: - >>> geo.name_matches("nametos*") - >>> geo.name_matches("*pier*") - >>> geo.name_matches("Abutment??", case_senstive=False) - """ - if not case_sensitive: - return fnmatchcase(self.name.casefold(), pattern.casefold()) - return fnmatchcase(self.name, pattern) - - def _group_matches( - self, pattern: str, *, case_sensitive: bool = True - ) -> bool: - """Checks if the group_label matches a pattern. - - Arguments: - pattern (str): the string pattern to be checked - - Keyword Arguments: - case_sensitive (bool, optional): if True (default) the check is - case sensitive. - - Returns: - (bool): Returns True if the group_label matches the pattern. - - Note: - The matching permits to use: - - "*" any chars - - "?" single char - - "[abc]" character set - - Examples: - >>> geo.group_matches("nametos*") - >>> geo.group_matches("*pier*") - >>> geo.group_matches("Abutment??", case_senstive=False) - """ - if self.group_label is None: - return False - if not case_sensitive: - return fnmatchcase(self.group_label.casefold(), pattern.casefold()) - return fnmatchcase(self.group_label, pattern) - - class PointGeometry(Geometry): """Class for a point geometry with material. @@ -395,6 +265,17 @@ def from_geometry( group_label=geo._group_label, ) + def __add__(self, other: Geometry) -> CompoundGeometry: + """Add operator "+" for geometries. + + Arguments: + other (Geometry): The other geometry to add. + + Returns: + CompoundGeometry: A new CompoundGeometry. + """ + return CompoundGeometry([self, other]) + def create_line_point_angle( point: t.Union[Point, t.Tuple[float, float]], @@ -699,6 +580,17 @@ def split_two_lines( # get the intersection return self.polygon.intersection(lines_polygon) + def __add__(self, other: Geometry) -> CompoundGeometry: + """Add operator "+" for geometries. + + Arguments: + other (Geometry): The other geometry to add. + + Returns: + CompoundGeometry: A new CompoundGeometry. + """ + return CompoundGeometry([self, other]) + def __sub__(self, other: Geometry) -> SurfaceGeometry: """Add operator "-" for geometries. @@ -1088,6 +980,17 @@ def mirror(self, axis: LineString) -> CompoundGeometry: processed_geoms.append(pg.mirror(axis)) return CompoundGeometry(geometries=processed_geoms) + def __add__(self, other: Geometry) -> CompoundGeometry: + """Add operator "+" for geometries. + + Arguments: + other (Geometry): The other geometry to add. + + Returns: + CompoundGeometry: A new CompoundGeometry. + """ + return CompoundGeometry([self, other]) + def __sub__(self, other: Geometry) -> CompoundGeometry: """Add operator "-" for geometries. diff --git a/structuralcodes/geometry/_rectangular.py b/structuralcodes/geometry/_rectangular.py index 21974945..058d0516 100644 --- a/structuralcodes/geometry/_rectangular.py +++ b/structuralcodes/geometry/_rectangular.py @@ -11,8 +11,7 @@ from numpy.typing import ArrayLike from shapely import Polygon -from structuralcodes.core.base import Material - +from ..core.base import Material from ._geometry import SurfaceGeometry diff --git a/structuralcodes/geometry/_reinforcement.py b/structuralcodes/geometry/_reinforcement.py index 05faf7b5..f67e467d 100644 --- a/structuralcodes/geometry/_reinforcement.py +++ b/structuralcodes/geometry/_reinforcement.py @@ -6,8 +6,7 @@ import numpy as np from shapely import Point -from structuralcodes.core.base import Material - +from ..core.base import Material from ._geometry import CompoundGeometry, PointGeometry, SurfaceGeometry diff --git a/structuralcodes/geometry/profiles/_base_profile.py b/structuralcodes/geometry/profiles/_base_profile.py index e1816b21..366b8be7 100644 --- a/structuralcodes/geometry/profiles/_base_profile.py +++ b/structuralcodes/geometry/profiles/_base_profile.py @@ -9,7 +9,7 @@ from shapely.affinity import rotate, translate from shapely.ops import split -from structuralcodes.core._marin_integration import marin_integration +from ...core._marin_integration import marin_integration class BaseProfile: diff --git a/structuralcodes/materials/concrete/_concrete.py b/structuralcodes/materials/concrete/_concrete.py index 4e9890a3..8d2d9b75 100644 --- a/structuralcodes/materials/concrete/_concrete.py +++ b/structuralcodes/materials/concrete/_concrete.py @@ -3,7 +3,7 @@ import abc import typing as t -from structuralcodes.core.base import ConstitutiveLaw, Material +from ...core.base import ConstitutiveLaw, Material class Concrete(Material): diff --git a/structuralcodes/materials/reinforcement/_reinforcement.py b/structuralcodes/materials/reinforcement/_reinforcement.py index 1e829e8e..d388c00c 100644 --- a/structuralcodes/materials/reinforcement/_reinforcement.py +++ b/structuralcodes/materials/reinforcement/_reinforcement.py @@ -3,7 +3,7 @@ import abc import typing as t -from structuralcodes.core.base import ConstitutiveLaw, Material +from ...core.base import ConstitutiveLaw, Material class Reinforcement(Material): diff --git a/structuralcodes/sections/_generic.py b/structuralcodes/sections/_generic.py index ccfe7ff9..cfb89002 100644 --- a/structuralcodes/sections/_generic.py +++ b/structuralcodes/sections/_generic.py @@ -13,8 +13,6 @@ from shapely.ops import unary_union import structuralcodes.core._section_results as s_res -from structuralcodes.core.base import Section, SectionCalculator -from structuralcodes.core.errors import NoConvergenceWarning from structuralcodes.geometry import ( CompoundGeometry, PointGeometry, @@ -22,6 +20,8 @@ ) from structuralcodes.materials.basic import ElasticMaterial +from ..core.base import Section, SectionCalculator +from ..core.errors import NoConvergenceWarning from .section_integrators import SectionIntegrator, integrator_factory diff --git a/structuralcodes/sections/section_integrators/_fiber_integrator.py b/structuralcodes/sections/section_integrators/_fiber_integrator.py index 6303e4ee..46ecbfab 100644 --- a/structuralcodes/sections/section_integrators/_fiber_integrator.py +++ b/structuralcodes/sections/section_integrators/_fiber_integrator.py @@ -9,9 +9,9 @@ from numpy.typing import ArrayLike, NDArray from shapely import Polygon -from structuralcodes.core.base import ConstitutiveLaw from structuralcodes.geometry import CompoundGeometry, SurfaceGeometry +from ...core.base import ConstitutiveLaw from ._section_integrator import SectionIntegrator diff --git a/structuralcodes/sections/section_integrators/_marin_integrator.py b/structuralcodes/sections/section_integrators/_marin_integrator.py index 58e4d1dd..c195e168 100644 --- a/structuralcodes/sections/section_integrators/_marin_integrator.py +++ b/structuralcodes/sections/section_integrators/_marin_integrator.py @@ -11,13 +11,13 @@ from shapely.geometry.polygon import orient from structuralcodes.core._marin_integration import marin_integration -from structuralcodes.core.base import ConstitutiveLaw from structuralcodes.geometry import ( CompoundGeometry, SurfaceGeometry, create_line_point_angle, ) +from ...core.base import ConstitutiveLaw from ._section_integrator import SectionIntegrator From 3cdc99de79aa224815cac9d8187358786f9ceb20 Mon Sep 17 00:00:00 2001 From: Morten Engen Date: Tue, 17 Mar 2026 22:55:40 +0100 Subject: [PATCH 3/3] Annotate result classes with base classes for geometry and section --- 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 3dbde253..29cc4df5 100644 --- a/structuralcodes/core/_section_results.py +++ b/structuralcodes/core/_section_results.py @@ -10,6 +10,8 @@ from numpy.typing import ArrayLike, NDArray from shapely import Point +from .base import Geometry, Section + @dataclass class SectionProperties: @@ -154,7 +156,7 @@ def isclose(self, other, rtol=1e-5, atol=1e-8): def _matching_geometries( - section, # For now we can't type hint the Section or even GenericSection due to circular import. #noqa: E501 + section: Section, name: t.Optional[str] = None, group_label: t.Optional[str] = None, case_sensitive: bool = True, @@ -186,7 +188,7 @@ def _strain_from_kinematics(eps_a, chi_y, chi_z, y: float, z: float): return eps_a - chi_z * y + chi_y * z -def _point_inside_geometry_surface(geo, y: float, z: float) -> bool: +def _point_inside_geometry_surface(geo: Geometry, y: float, z: float) -> bool: """Check if a point is inside or on the boundary of a surface geometry.""" p = Point(y, z) return geo.polygon.contains(p) or geo.polygon.touches(p) @@ -205,7 +207,7 @@ def _point_matches_geometry_point(gp, y: float, z: float) -> bool: def _get_point_response( - section, # For now we can't type hint the Section or even GenericSection due to circular import. #noqa: E501 + section: Section, eps_a: float, chi_y: float, chi_z: float, @@ -282,7 +284,7 @@ class MomentCurvatureResults: m_y: ArrayLike = None # the moment m_z: ArrayLike = None # the moment - section = None + section: Section = None detailed_result: SectionDetailedResultState = None seed: int = None