diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml
index 482242ed..7dfbba64 100644
--- a/.github/workflows/build.yaml
+++ b/.github/workflows/build.yaml
@@ -6,6 +6,7 @@ on:
pull_request:
branches:
- 'dev'
+ - 'implement-shell-section'
jobs:
build:
runs-on: ubuntu-latest
diff --git a/structuralcodes/geometry/__init__.py b/structuralcodes/geometry/__init__.py
index f8b13953..1bffb0b4 100644
--- a/structuralcodes/geometry/__init__.py
+++ b/structuralcodes/geometry/__init__.py
@@ -15,6 +15,7 @@
add_reinforcement_circle,
add_reinforcement_line,
)
+from ._shell_geometry import ShellGeometry, ShellReinforcement
__all__ = [
'Geometry',
@@ -28,4 +29,6 @@
'CircularGeometry',
'add_reinforcement_circle',
'RectangularGeometry',
+ 'ShellGeometry',
+ 'ShellReinforcement',
]
diff --git a/structuralcodes/geometry/_shell_geometry.py b/structuralcodes/geometry/_shell_geometry.py
new file mode 100644
index 00000000..b4ae44ce
--- /dev/null
+++ b/structuralcodes/geometry/_shell_geometry.py
@@ -0,0 +1,437 @@
+"""A concrete implementation of a shell geometry."""
+
+import typing as t
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+from ..core.base import Material
+from ._geometry import Geometry
+
+
+class ShellReinforcement(Geometry):
+ """A class for representing reinforcement in a shell geometry."""
+
+ _z: float
+ _n_bars: float
+ _cc_bars: float
+ _diameter_bar: float
+ _material: Material
+ _phi: float
+ _area: t.Optional[float]
+ _T: t.Optional[ArrayLike]
+ _local_stress_to_global_force: t.Optional[ArrayLike]
+ _local_modulus_to_global_stiffness: t.Optional[ArrayLike]
+
+ def __init__(
+ self,
+ z: float,
+ n_bars: float,
+ cc_bars: float,
+ diameter_bar: float,
+ material: Material,
+ phi: float,
+ name: t.Optional[str] = None,
+ group_label: t.Optional[str] = None,
+ ) -> None:
+ """Initialize a shell reinforcement."""
+ super().__init__(name, group_label)
+
+ if not isinstance(material, Material):
+ raise TypeError(
+ 'Material should be a valid structuralcodes.base.Material'
+ f' object. {repr(material)}'
+ )
+ self._z = z
+ self._n_bars = n_bars
+ self._cc_bars = cc_bars
+ self._diameter_bar = diameter_bar
+ self._material = material
+ self._phi = phi
+ self._area = None
+ self._T = None
+ self._local_stress_to_global_force = None
+ self._local_modulus_to_global_stiffness = None
+
+ @property
+ def z(self) -> float:
+ """Return the reinforcement position over the thickness."""
+ return self._z
+
+ @property
+ def n_bars(self) -> float:
+ """Return the number of bars per bundle."""
+ return self._n_bars
+
+ @property
+ def cc_bars(self) -> float:
+ """Return the spacing between bars."""
+ return self._cc_bars
+
+ @property
+ def diameter_bar(self) -> float:
+ """Return the bar diameter."""
+ return self._diameter_bar
+
+ @property
+ def material(self) -> Material:
+ """Return the material of the reinforcement."""
+ return self._material
+
+ @property
+ def phi(self) -> float:
+ """Return the orientation angle of the reinforcement."""
+ return self._phi
+
+ @property
+ def T(self) -> np.ndarray:
+ """Return the transformation matrix for the reinforcement."""
+ if self._T is None:
+ c, s = np.cos(self.phi), np.sin(self.phi)
+ self._T = np.array(
+ [
+ [c * c, s * s, c * s],
+ [s * s, c * c, -c * s],
+ [-2 * c * s, 2 * c * s, c * c - s * s],
+ ]
+ )
+ return self._T
+
+ @property
+ def local_stress_to_global_force(self) -> np.ndarray:
+ """Return the array to multiply with a local stress to obtain
+ contribution the global stress resultant.
+ """
+ if self._local_stress_to_global_force is None:
+ self._local_stress_to_global_force = self.area * self.T.T[:, 0]
+ return self._local_stress_to_global_force
+
+ @property
+ def local_modulus_to_global_stiffness(self) -> np.ndarray:
+ """Return the matrix to multiply with the local modulus to obtain the
+ contribution to the global stiffness.
+ """
+ if self._local_modulus_to_global_stiffness is None:
+ self._local_modulus_to_global_stiffness = (
+ self.T.T @ np.diag([1, 0, 0]) @ self.T * self.area
+ )
+ return self._local_modulus_to_global_stiffness
+
+ @property
+ def area(self) -> float:
+ """Return the reinforcement area per unit width."""
+ if self._area is None:
+ self._area = (
+ self.n_bars
+ * np.pi
+ * (self.diameter_bar / 2) ** 2
+ / self.cc_bars
+ )
+ return self._area
+
+ def _repr_svg_(self) -> str:
+ """Returns the svg representation."""
+ raise NotImplementedError
+
+
+class ShellGeometry(Geometry):
+ """A class for a shell with a thickness and material."""
+
+ _reinforcement: t.List[ShellReinforcement]
+
+ def __init__(
+ self,
+ thickness: float,
+ material: Material,
+ name: t.Optional[str] = None,
+ group_label: t.Optional[str] = None,
+ ) -> None:
+ """Initialize a shell geometry."""
+ super().__init__(name=name, group_label=group_label)
+
+ if thickness <= 0:
+ raise ValueError('Shell thickness must be positive.')
+
+ if not isinstance(material, Material):
+ raise TypeError(
+ 'Material should be a valid structuralcodes.base.Material'
+ f' object. {repr(material)}'
+ )
+
+ self._thickness = thickness
+ self._material = material
+
+ self._reinforcement = []
+
+ @property
+ def thickness(self) -> float:
+ """Return the shell thickness."""
+ return self._thickness
+
+ @property
+ def material(self) -> Material:
+ """Return the material of the shell."""
+ return self._material
+
+ @property
+ def reinforcement(self) -> t.List[ShellReinforcement]:
+ """Return all reinforcement layers."""
+ return self._reinforcement
+
+ def add_reinforcement(
+ self,
+ reinforcement: t.Union[ShellReinforcement, t.List[ShellReinforcement]],
+ ) -> None:
+ """Add reinforcement to the shell geometry."""
+ if isinstance(reinforcement, ShellReinforcement):
+ self._reinforcement.append(reinforcement)
+ to_validate = [reinforcement]
+ elif isinstance(reinforcement, list):
+ self._reinforcement.extend(reinforcement)
+ to_validate = reinforcement
+ else:
+ raise TypeError(
+ 'Reinforcement must be a ShellReinforcement or a list of '
+ 'ShellReinforcement.'
+ )
+
+ # Validate each reinforcement layer
+ for r in to_validate:
+ half_thickness = self._thickness / 2
+ if not (
+ -half_thickness + r.diameter_bar / 2
+ <= r.z
+ <= half_thickness - r.diameter_bar / 2
+ ):
+ raise ValueError(
+ f'Reinforcement at z = {r.z:.2f} mm is outside the '
+ f'range [-{half_thickness:.2f}, {half_thickness:.2f}] mm.'
+ )
+
+ def _repr_svg_(self) -> str: # noqa: PLR0915
+ """Returns the svg representation."""
+ # overall drawing area dimensions
+ w, h = 1000, 1000
+ t_val = self.thickness
+ # half-width/height for centering
+ hw, hh = w / 2, h / 2
+ # single view dimensions
+ sw, sh = w + 200, h + 200
+ # side view height based on thickness
+ sh_side = int(t_val + 200)
+ # main and side view ViewBox definitions
+ vb_main = f'{-hw - 100} {-hh - 100} {w + 200} {h + 200}'
+ vb_side = f'{-hw - 100} {-t_val / 2 - 100} {w + 200} {t_val + 200}'
+ # scaling and spacing between views
+ scale, gap = 0.6, 50
+
+ def draw_rebars(view: str) -> str:
+ # draw parallel rebar lines for top/bottom plan views
+ lines: t.List[str] = []
+ for layer in self.reinforcement:
+ # select only top or bottom layers
+ top = view == 'top' and layer.z >= 0
+ bot = view == 'bottom' and layer.z < 0
+ if not (top or bot):
+ continue
+ # orientation and spacing
+ phi, sp = layer.phi, layer.cc_bars
+ c, s = np.cos(phi), np.sin(phi)
+ # perpendicular offset vector
+ px, py = -s, c
+ # color by orientation
+ col = (
+ 'red'
+ if np.isclose(phi % np.pi, 0)
+ else 'blue'
+ if np.isclose(phi % np.pi, np.pi / 2)
+ else 'green'
+ )
+ n = int(w / sp) + 3 # linecount with safety margin to angles
+ L = 2000 # line length extension
+ for i in range(-n // 2, n // 2 + 1):
+ ox, oy = i * sp * px, i * sp * py
+ # line endpoints
+ x1, y1 = ox - L * c, oy - L * s
+ x2, y2 = ox + L * c, oy + L * s
+ for _ in range(int(layer.n_bars)):
+ # actual bar offset
+ lines.append(
+ f''
+ )
+ return ''.join(lines)
+
+ def build_plan(view: str) -> str:
+ # define clipping area for concrete
+ clip = (
+ ''
+ f''
+ ''
+ )
+ # concrete background
+ bg = (
+ f''
+ )
+ # concrete outline
+ out = (
+ f''
+ )
+ # view label
+ lbl = (
+ f'' + f'{view} view'
+ )
+ # rebar lines clipped to concrete area
+ body = (
+ ''
+ + draw_rebars(view)
+ + ''
+ )
+ return clip + bg + body + out + lbl
+
+ def build_side(ax: str) -> str:
+ # side view perpendicular direction
+ span = w if ax == 'x' else h
+ half = span / 2
+ t_ht = self.thickness
+ cid = 'clipSide' + ax
+ # clipping region for shell thickness
+ clip = (
+ ''
+ f''
+ ''
+ )
+ # shell background
+ bg = (
+ f''
+ )
+ # shell outline
+ out = (
+ f''
+ )
+ # side view label
+ lbl = (
+ f'' + f'Side {ax.upper()}-view'
+ )
+ elems: t.List[str] = []
+ for layer in self.reinforcement:
+ # select bars not parallel to this side-plane
+ perp = not (
+ (
+ ax == 'x'
+ and np.isclose(layer.phi % np.pi, np.pi / 2, atol=1e-2)
+ )
+ or (
+ ax == 'y'
+ and np.isclose(layer.phi % np.pi, 0, atol=1e-2)
+ )
+ )
+ if perp:
+ # draw circles for perp layers
+ step = layer.cc_bars or span
+ kmax = int(np.floor(half / step))
+ xs = [k * step for k in range(-kmax, kmax + 1)]
+ col = (
+ 'red'
+ if np.isclose(layer.phi % np.pi, 0)
+ else 'blue'
+ if np.isclose(layer.phi % np.pi, np.pi / 2)
+ else 'green'
+ )
+ for x in xs:
+ for j in range(int(layer.n_bars)):
+ ex = (
+ j - (layer.n_bars - 1) / 2
+ ) * layer.diameter_bar
+ cy = -layer.z
+ elems.append(
+ f''
+ )
+ else:
+ # draw lines for parallel layers
+ col = (
+ 'red'
+ if np.isclose(layer.phi % np.pi, 0)
+ else 'blue'
+ if np.isclose(layer.phi % np.pi, np.pi / 2)
+ else 'green'
+ )
+ y = -layer.z
+ for j in range(int(layer.n_bars)):
+ ex = (j - (layer.n_bars - 1) / 2) * layer.diameter_bar
+ x1 = -half + ex
+ x2 = half + ex
+ elems.append(
+ f''
+ )
+ body = (
+ '' + ''.join(elems) + ''
+ )
+ return clip + bg + body + out + lbl
+
+ # assemble all SVG fragments
+ svg = [
+ f'')
+ return ''.join(svg)
diff --git a/structuralcodes/materials/constitutive_laws/__init__.py b/structuralcodes/materials/constitutive_laws/__init__.py
index bacf6b75..4b018984 100644
--- a/structuralcodes/materials/constitutive_laws/__init__.py
+++ b/structuralcodes/materials/constitutive_laws/__init__.py
@@ -4,7 +4,16 @@
from ...core.base import ConstitutiveLaw, Material
from ._bilinearcompression import BilinearCompression
+from ._concrete_smeared_cracking import (
+ ConcreteSmearedCracking,
+ ConstantPoissonReduction,
+ GeneralVecchioCollins,
+ NoTension,
+ calculate_principal_strains,
+ establish_strain_transformation_matrix,
+)
from ._elastic import Elastic
+from ._elastic_2d import Elastic2D
from ._elasticplastic import ElasticPlastic
from ._initial_strain import InitialStrain
from ._parabolarectangle import ParabolaRectangle
@@ -14,6 +23,7 @@
__all__ = [
'Elastic',
+ 'Elastic2D',
'ElasticPlastic',
'ParabolaRectangle',
'BilinearCompression',
@@ -23,10 +33,17 @@
'InitialStrain',
'get_constitutive_laws_list',
'create_constitutive_law',
+ 'ConcreteSmearedCracking',
+ 'NoTension',
+ 'ConstantPoissonReduction',
+ 'GeneralVecchioCollins',
+ 'calculate_principal_strains',
+ 'establish_strain_transformation_matrix',
]
CONSTITUTIVE_LAWS: t.Dict[str, ConstitutiveLaw] = {
'elastic': Elastic,
+ 'elastic2d': Elastic2D,
'elasticplastic': ElasticPlastic,
'elasticperfectlyplastic': ElasticPlastic,
'bilinearcompression': BilinearCompression,
diff --git a/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/__init__.py b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/__init__.py
new file mode 100644
index 00000000..92a3adf0
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/__init__.py
@@ -0,0 +1,25 @@
+"""Classes for concrete smeared cracking."""
+
+from ._core import ConcreteSmearedCracking
+from ._cracking_criterion import CrackingCriterion, NoTension
+from ._poisson_reduction import ConstantPoissonReduction, PoissonReduction
+from ._strength_reduction_lateral_cracking import (
+ GeneralVecchioCollins,
+ StrengthReductionLateralCracking,
+)
+from ._utils import (
+ calculate_principal_strains,
+ establish_strain_transformation_matrix,
+)
+
+__all__ = [
+ 'ConcreteSmearedCracking',
+ 'CrackingCriterion',
+ 'NoTension',
+ 'PoissonReduction',
+ 'ConstantPoissonReduction',
+ 'StrengthReductionLateralCracking',
+ 'GeneralVecchioCollins',
+ 'calculate_principal_strains',
+ 'establish_strain_transformation_matrix',
+]
diff --git a/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_core.py b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_core.py
new file mode 100644
index 00000000..104bc820
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_core.py
@@ -0,0 +1,185 @@
+"""The core class for representing concrete smeared cracking."""
+
+from __future__ import annotations # To have clean hints of ArrayLike in docs
+
+import typing as t
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+from structuralcodes.core.base import ConstitutiveLaw
+
+from ._cracking_criterion import CrackingCriterion, NoTension
+from ._poisson_reduction import PoissonReduction
+from ._strength_reduction_lateral_cracking import (
+ StrengthReductionLateralCracking,
+)
+from ._utils import (
+ calculate_principal_strains,
+ establish_strain_transformation_matrix,
+)
+
+
+class ConcreteSmearedCracking:
+ """A plane stress model representing smeared cracking in reinforced
+ concrete.
+ """
+
+ _initial_modulus_compression: t.Optional[float] = None
+ _initial_modulus_tension: t.Optional[float] = None
+ _uniaxial_compression: ConstitutiveLaw
+ _strength_reduction_lateral_cracking: StrengthReductionLateralCracking
+ _poisson_reduction: PoissonReduction
+ _cracking_criterion: CrackingCriterion
+
+ def __init__(
+ self,
+ uniaxial_compression: ConstitutiveLaw,
+ strength_reduction_lateral_cracking: StrengthReductionLateralCracking,
+ poisson_reduction: PoissonReduction,
+ cracking_criterion: t.Optional[CrackingCriterion] = None,
+ ):
+ """Initialize the model."""
+ self._uniaxial_compression = uniaxial_compression
+ self._strength_reduction_lateral_cracking = (
+ strength_reduction_lateral_cracking
+ )
+ self._poisson_reduction = poisson_reduction
+ self._cracking_criterion = cracking_criterion or NoTension()
+
+ def get_stress(
+ self, eps: ArrayLike, return_global: bool = True
+ ) -> ArrayLike:
+ """Calculate the stresses based on strains.
+
+ Arguments:
+ eps (ArrayLike): The strains epsilon_x, epsilon_y and gamma_xy.
+
+ """
+ # Calculate principal strains and principal strain direction
+ eps_p, phi = calculate_principal_strains(strain=eps)
+
+ # Establish strain transformation matrix
+ T = establish_strain_transformation_matrix(phi=phi)
+
+ # Include the influence of Poisson's ratio on the principal strains.
+ eps_pf = self.calculate_effective_principal_strains(eps_p)
+
+ # Compressive-strength reduction factor due to lateral tension.
+ beta = self.strength_reduction_lateral_cracking.reduction(eps_pf)
+
+ # Compute the principal stresses from the uniaxial compression law.
+ sig_p = self.uniaxial_compression.get_stress(eps_pf)
+
+ # Apply the compressive-strength reduction factor
+ sig_p[sig_p < 0] *= beta
+
+ if not return_global:
+ return np.array([*sig_p, 0])
+
+ # Transform back to global coords
+ return T.T @ np.array([*sig_p, 0])
+
+ def get_secant(
+ self, eps: ArrayLike, return_global: bool = True
+ ) -> np.ndarray:
+ # Calculate principal strains and principal strain direction
+ eps_p, phi = calculate_principal_strains(strain=eps)
+
+ # Establish strain transformation matrix
+ T = establish_strain_transformation_matrix(phi=phi)
+
+ # Poisson correction
+ eps_pf = self.calculate_effective_principal_strains(eps_p)
+ cracked = self.cracking_criterion.cracked(eps_p=eps_pf)
+
+ # Establish the diagonal of material stiffness matrix
+ D = self.uniaxial_compression.get_secant(eps_pf)
+
+ # Compressive-strength reduction factor due to lateral tension.
+ beta = self.strength_reduction_lateral_cracking.reduction(eps_pf)
+
+ # Scale down compressive stiffness with strength reduction factor
+ D[eps_pf < 0] *= beta
+
+ # Correct for the poisson effect
+ C = self.poisson_reduction.poisson_matrix(cracked=cracked) @ np.diag(D)
+
+ # Ensure symmetry
+ C = 0.5 * (C + C.T)
+
+ # Shear modulus
+ G_value = D.mean() / (
+ 2 * (1 + self.poisson_reduction.nu(cracked=cracked))
+ )
+ G_value = max(G_value, 1e-4 * self.initial_modulus_compression)
+ G_12 = np.array([[G_value]]) # Shape (1, 1)
+
+ Z_21 = np.zeros((2, 1)) # Shape (2, 1)
+
+ # 3x3 secant stiffness matrix
+ Cp = np.block(
+ [
+ [C, Z_21],
+ [Z_21.T, G_12],
+ ]
+ )
+ if not return_global:
+ return Cp
+
+ return T.T @ Cp @ T
+
+ def get_tangent():
+ pass
+
+ def calculate_effective_principal_strains(
+ self, eps_p: ArrayLike
+ ) -> ArrayLike:
+ """Compute the effective principal strains to include the influence of
+ Poisson's ratio. Taken from 'Nonlinear Analysis of Reinforced-Concrete
+ Shells' by M. A. Polak and F. J. Vecchio (1993).
+ """
+ # Calculate the effective principal strains assuming uncracked
+ eps_pf_uncracked = (
+ self.poisson_reduction.poisson_matrix(cracked=False) @ eps_p
+ )
+
+ # If the assumption of uncracked holds after calculating effective
+ # strains, return the effective uncracked principal strains
+ if not self.cracking_criterion.cracked(eps_p=eps_pf_uncracked):
+ return eps_pf_uncracked
+
+ # Else, return the effective principal strains assuming cracked
+ return self.poisson_reduction.poisson_matrix(cracked=True) @ eps_p
+
+ @property
+ def uniaxial_compression(self) -> ConstitutiveLaw:
+ """Return the constitutive law for uniaxial compression."""
+ return self._uniaxial_compression
+
+ @property
+ def strength_reduction_lateral_cracking(
+ self,
+ ) -> StrengthReductionLateralCracking:
+ """Return the model for strength reduction due to lateral cracking."""
+ return self._strength_reduction_lateral_cracking
+
+ @property
+ def poisson_reduction(self) -> PoissonReduction:
+ """Return the model for reduction of Poisson ratio due to cracking."""
+ return self._poisson_reduction
+
+ @property
+ def cracking_criterion(self) -> CrackingCriterion:
+ """Return the cracking criterion."""
+ return self._cracking_criterion
+
+ @property
+ def initial_modulus_compression(self) -> float:
+ """Return the initial modulus in uniaxial compression."""
+ if self._initial_modulus_compression is None:
+ self._initial_modulus_compression = (
+ self.uniaxial_compression.get_tangent(0.0)
+ )
+
+ return self._initial_modulus_compression
diff --git a/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_cracking_criterion.py b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_cracking_criterion.py
new file mode 100644
index 00000000..567cdd09
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_cracking_criterion.py
@@ -0,0 +1,24 @@
+"""Cracking criteria for concrete smeared cracking."""
+
+from __future__ import annotations # To have clean hints of ArrayLike in docs
+
+from abc import ABC, abstractmethod
+
+from numpy.typing import ArrayLike
+
+
+class CrackingCriterion(ABC):
+ """An abstract based class for modelling a cracking criterion."""
+
+ @abstractmethod
+ def cracked(self, eps_p: ArrayLike) -> bool:
+ """Check if the material is cracked based on principal strains."""
+ raise NotImplementedError
+
+
+class NoTension(CrackingCriterion):
+ """Cracking criterion to represent a material with no tensile strength."""
+
+ def cracked(self, eps_p: ArrayLike) -> bool:
+ """Check if the material is cracked based on the principal strains."""
+ return max(eps_p) > 0
diff --git a/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_poisson_reduction.py b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_poisson_reduction.py
new file mode 100644
index 00000000..1eec7e5b
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_poisson_reduction.py
@@ -0,0 +1,79 @@
+"""Reduction of Poisson ratio in concrete smeared cracking."""
+
+from __future__ import annotations # To have clean hints of ArrayLike in docs
+
+import typing as t
+from abc import ABC, abstractmethod
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+
+class PoissonReduction(ABC):
+ """An abstract base class for modelling reduction of the Poisson ratio due
+ to cracking.
+ """
+
+ @abstractmethod
+ def nu(
+ self,
+ cracked: bool,
+ eps: ArrayLike,
+ eps_p: ArrayLike,
+ eps_pf: ArrayLike,
+ ):
+ raise NotImplementedError
+
+ @abstractmethod
+ def poisson_matrix(
+ self,
+ cracked: bool,
+ eps: ArrayLike,
+ eps_p: ArrayLike,
+ eps_pf: ArrayLike,
+ ):
+ raise NotImplementedError
+
+
+class ConstantPoissonReduction(PoissonReduction):
+ """A material with a constant reduction in Poisson ratio due to
+ cracking.
+ """
+
+ _initial_nu: float
+ _cracked_nu: float
+ _poisson_matrix_cracked: t.Optional[np.ndarray]
+ _poisson_matrix_uncracked: t.Optional[np.ndarray]
+
+ def __init__(self, initial_nu: float, cracked_nu: float = 0.0):
+ """Initialize the model."""
+ self._initial_nu = initial_nu
+ self._cracked_nu = cracked_nu
+ self._poisson_matrix_cracked = None
+ self._poisson_matrix_uncracked = None
+
+ def nu(self, cracked: bool, *args, **kwargs):
+ """Return the Poisson ratio."""
+ del args, kwargs
+ return self._initial_nu if not cracked else self._cracked_nu
+
+ def poisson_matrix(self, cracked: bool, *args, **kwargs) -> np.ndarray:
+ """Return the 2x2 Poisson matrix."""
+ del args, kwargs
+
+ # Return Poisson matrix
+ if cracked:
+ if self._poisson_matrix_cracked is None:
+ self._poisson_matrix_cracked = np.eye(2)
+ return self._poisson_matrix_cracked
+
+ if self._poisson_matrix_uncracked is None:
+ self._poisson_matrix_uncracked = (
+ 1 / (1 - self._initial_nu**2)
+ ) * np.array(
+ [
+ [1, self._initial_nu],
+ [self._initial_nu, 1],
+ ]
+ )
+ return self._poisson_matrix_uncracked
diff --git a/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_strength_reduction_lateral_cracking.py b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_strength_reduction_lateral_cracking.py
new file mode 100644
index 00000000..4a47f69a
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_strength_reduction_lateral_cracking.py
@@ -0,0 +1,59 @@
+"""Strength reduction due to lateral cracking for concrete smeared cracking."""
+
+from __future__ import annotations # To have clean hints of ArrayLike in docs
+
+from abc import ABC, abstractmethod
+
+from numpy.typing import ArrayLike
+
+
+class StrengthReductionLateralCracking(ABC):
+ """Base class for models for accounting for strength reduction due to
+ lateral cracking.
+ """
+
+ @abstractmethod
+ def reduction(self, strains: ArrayLike) -> float:
+ """Calculate the strength reduction based on a set of strains."""
+ raise NotImplementedError
+
+
+class GeneralVecchioCollins(StrengthReductionLateralCracking):
+ """General model for accounting for strength reduction due to lateral
+ cracking based on the work of Vecchio and Collins.
+ """
+
+ _c_1: float
+ _c_2: float
+ _minimum_value: float
+
+ def __init__(self, c_1: float, c_2: float, minimum_value: float = 0.0):
+ """Initialize the model."""
+ self._c_1 = c_1
+ self._c_2 = c_2
+ self._minimum_value = minimum_value
+
+ @property
+ def c_1(self) -> float:
+ """Return the coefficient c_1."""
+ return self._c_1
+
+ @property
+ def c_2(self) -> float:
+ """Return the coefficient c_2."""
+ return self._c_2
+
+ @property
+ def minimum_value(self) -> float:
+ """Return the minimum value."""
+ return self._minimum_value
+
+ def reduction(self, strains: ArrayLike) -> float:
+ """Calculate the strength reduction based on a set of principal
+ strains.
+
+ beta = 1 / (c_1 + c_2 * eps_1)
+ """
+ beta = 1 / (self.c_1 + self.c_2 * max(strains))
+
+ return max(min(beta, 1), self.minimum_value)
diff --git a/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_utils.py b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_utils.py
new file mode 100644
index 00000000..e059f9fe
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_concrete_smeared_cracking/_utils.py
@@ -0,0 +1,40 @@
+"""Utility functions for concrete smeared cracking."""
+
+from __future__ import annotations # To have clean hints of ArrayLike in docs
+
+import typing as t
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+
+def calculate_principal_strains(
+ strain: ArrayLike,
+) -> t.Tuple[ArrayLike, float]:
+ """Calculate the principal strains in 2D."""
+ eps = np.atleast_1d(strain)
+ # Use arctan2 to obtain the angle in the correct quadrant
+ phi = 0.5 * np.arctan2(eps[2], eps[0] - eps[1])
+ eps_p = np.array(
+ [
+ (eps[0] + eps[1]) / 2
+ + (((eps[0] - eps[1]) / 2) ** 2 + 0.25 * eps[2] ** 2) ** 0.5,
+ (eps[0] + eps[1]) / 2
+ - (((eps[0] - eps[1]) / 2) ** 2 + 0.25 * eps[2] ** 2) ** 0.5,
+ ]
+ )
+
+ return eps_p, phi
+
+
+def establish_strain_transformation_matrix(phi: float) -> ArrayLike:
+ """Establish the strain transformation matrix for a given angle."""
+ c = np.cos(phi)
+ s = np.sin(phi)
+ return np.array(
+ [
+ [c * c, s * s, s * c],
+ [s * s, c * c, -s * c],
+ [-2 * s * c, 2 * s * c, c * c - s * s],
+ ]
+ )
diff --git a/structuralcodes/materials/constitutive_laws/_elastic_2d.py b/structuralcodes/materials/constitutive_laws/_elastic_2d.py
new file mode 100644
index 00000000..9548dc6c
--- /dev/null
+++ b/structuralcodes/materials/constitutive_laws/_elastic_2d.py
@@ -0,0 +1,86 @@
+import typing as t
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+from ._elastic import Elastic
+
+
+class Elastic2D(Elastic):
+ """Class for elastic constitutive law for 2D operations."""
+
+ __materials__: t.Tuple[str] = (
+ 'concrete',
+ 'steel',
+ 'rebars',
+ )
+
+ def __init__(
+ self,
+ E: float,
+ nu: float,
+ name: t.Optional[str] = None,
+ ) -> None:
+ """Initialize an Elastic2D Material.
+
+ Arguments:
+ E (float): The elastic modulus.
+ nu (float): Poisson's ratio.
+ name (str, optional): A descriptive name for the constitutive law.
+ """
+ super().__init__(E=E, name=name)
+ self._nu = nu
+
+ self._C_matrix: t.Optional[ArrayLike] = None
+
+ @property
+ def E(self) -> float:
+ """Return the elastic modulus."""
+ return self._E
+
+ @property
+ def nu(self) -> float:
+ """Return Poisson's ratio.
+
+ Note:
+ Including this property allows future checks or transformations
+ (e.g., ensuring 0 ≤ nu < 0.5).
+ """
+ return self._nu
+
+ @property
+ def C_matrix(self) -> np.ndarray:
+ """Return the 2D constitutive matrix."""
+ if self._C_matrix is None:
+ E = self.E
+ nu = self.nu
+
+ self._C_matrix = (
+ E
+ / (1 - nu**2)
+ * np.array(
+ [
+ [1.0, nu, 0.0],
+ [nu, 1.0, 0.0],
+ [0.0, 0.0, (1 - nu) / 2.0],
+ ]
+ )
+ )
+ return self._C_matrix
+
+ def get_stress(self, eps: ArrayLike) -> np.ndarray:
+ """Return a 2D stress vector [sigma_x, sigma_y, tau_xy]
+ given a 2D strain vector [eps_x, epx_y, gamma_xy].
+ """
+ eps = np.atleast_1d(eps)
+ try:
+ return self.C_matrix @ eps
+ except ValueError as e:
+ raise ValueError(
+ 'The input strain vector must have a length of 3.'
+ ) from e
+
+ def get_secant(self, *args, **kwargs) -> np.ndarray:
+ """Return the 2D secant stiffness matrix."""
+ del args, kwargs
+ return self.C_matrix
diff --git a/structuralcodes/sections/__init__.py b/structuralcodes/sections/__init__.py
index 6a302251..3dc1184f 100644
--- a/structuralcodes/sections/__init__.py
+++ b/structuralcodes/sections/__init__.py
@@ -2,10 +2,12 @@
from ._generic import GenericSection, GenericSectionCalculator
from ._rc_utils import calculate_elastic_cracked_properties
+from ._shell_section import ShellSection, ShellSectionCalculator
from .section_integrators import (
FiberIntegrator,
MarinIntegrator,
SectionIntegrator,
+ ShellFiberIntegrator,
integrator_factory,
marin_integration,
)
@@ -19,4 +21,7 @@
'integrator_factory',
'marin_integration',
'calculate_elastic_cracked_properties',
+ 'ShellFiberIntegrator',
+ 'ShellSection',
+ 'ShellSectionCalculator',
]
diff --git a/structuralcodes/sections/_shell_section.py b/structuralcodes/sections/_shell_section.py
new file mode 100644
index 00000000..a0a112a0
--- /dev/null
+++ b/structuralcodes/sections/_shell_section.py
@@ -0,0 +1,243 @@
+import typing as t
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+from scipy.linalg import lu_factor, lu_solve
+
+from ..core.base import Section, SectionCalculator
+from ..geometry import ShellGeometry
+from .section_integrators import ShellFiberIntegrator
+
+
+class ShellSection(Section):
+ """This is the implementation of the shell class section.
+
+ Attributes:
+ geometry ShellGeometry: The geometry of
+ the section.
+ name (str): The name of the section.
+ """
+
+ def __init__(
+ self,
+ geometry: ShellGeometry,
+ name: t.Optional[str] = None,
+ **kwargs,
+ ) -> None:
+ """Initialize a ShellSection.
+
+ Note:
+ The ShellSection uses a ShellSectionCalculator for all
+ calculations. The ShellSectionCalculator uses a
+ ShellFiberIntegrator for integrating over the thickness.
+ """
+ if name is None:
+ name = 'ShellSection'
+ super().__init__(name)
+ self._geometry = geometry
+ self.section_calculator = ShellSectionCalculator(
+ section=self, **kwargs
+ )
+
+ @property
+ def geometry(self):
+ """Return the geometry of the section."""
+ return self._geometry
+
+
+class ShellSectionCalculator(SectionCalculator):
+ """A calculator for shell sections."""
+
+ section: ShellSection
+ mesh_size: float
+
+ def __init__(
+ self,
+ section: ShellSection,
+ **kwargs,
+ ) -> None:
+ """Initialize the ShellSectionCalculator.
+
+ Arguments:
+ section (ShellSection): The section object.
+ """
+ super().__init__(section=section)
+ self.integrator = ShellFiberIntegrator()
+ self.mesh_size = kwargs.get('mesh_size', 0.01)
+ self.layers: t.Optional[t.Tuple] = None
+
+ def _calculate_gross_section_properties(self):
+ pass
+
+ def calculate_bending_strength(self):
+ pass
+
+ def calculate_moment_curvature(self):
+ pass
+
+ def integrate_strain_profile(
+ self,
+ strain: ArrayLike,
+ integrate: t.Literal['stress', 'modulus'] = 'stress',
+ ) -> t.Union[t.Tuple[float, float, float, float, float, float], NDArray]:
+ """Integrate a strain profile returning stress resultants or tangent
+ section stiffness matrix.
+
+ Arguments:
+ strain (ArrayLike): Represents the deformation plane. The strain
+ should have six entries representing respectively: eps_x,
+ eps_y, gamma_xy, kappa_x, kappa_y, kappa_xy.
+ integrate (str): a string indicating the quantity to integrate over
+ the section. It can be 'stress' or 'modulus'. When 'stress'
+ is selected, the return value will be the stress resultants Nx,
+ Ny, Nxy, Mx, My and Mxy, while if 'modulus' is selected, the
+ return will be the tangent section stiffness matrix (default is
+ 'stress').
+
+ Returns:
+ Union(Tuple(float, float, float, float, float, float),NDArray):
+ Nx, Ny, Nxy, Mx, My and Mxy My when `integrate='stress'`, or a
+ numpy array representing the stiffness matrix then
+ `integrate='modulus'`.
+
+ Examples:
+ result = self.integrate_strain_profile(strain,integrate='modulus')
+ # `result` will be the tangent stiffness matrix (a 6x6 numpy array)
+
+ result = self.integrate_strain_profile(strain)
+ # `result` will be a tuple containing section forces (Nx, Ny, Nxy,
+ Mx, My, Mxy)
+
+ Raises:
+ ValueError: If a unkown value is passed to the `integrate`
+ parameter.
+ """
+ result = self.integrator.integrate_strain_response_on_geometry(
+ geo=self.section.geometry,
+ strain=strain,
+ integrate=integrate,
+ mesh_size=self.mesh_size,
+ layers=self.layers,
+ )
+
+ # Save layers for future use
+ if self.layers is None and result[-1] is not None:
+ self.layers = result[-1]
+
+ # Return the results without layers
+ if len(result) == 2:
+ # Return only the stiffness matrix
+ return result[0]
+
+ # Return only the stress resultants
+ return result[:-1]
+
+ def calculate_strain_profile(
+ self,
+ nx: float,
+ ny: float,
+ nxy: float,
+ mx: float,
+ my: float,
+ mxy: float,
+ initial: bool = False,
+ max_iter: int = 10,
+ tol: float = 1e-6,
+ ) -> t.List[float]:
+ """Computes the strain profile for a given set of stress resultants.
+
+ Args:
+ nx (float): Membrane force in x-direction.
+ ny (float): Membrane force in y-direction.
+ nxy (float): Membrane shear.
+ mx (float): Plate moment giving stresses in x-direction.
+ my (float): Plate moment giving stresses in y-direction.
+ mxy (float): Plate moment giving shear stresses in the xy-plane.
+
+ Returns:
+ list(float): The strain response eps_x, eps_y, gamma_xy, chi_x,
+ chi_y and chi_xy.
+ """
+ # Get the gometry
+ geom = self.section.geometry
+
+ # Collect loads in a numpy array
+ loads = np.array([nx, ny, nxy, mx, my, mxy])
+
+ # Compute initial tangent stiffness matrix
+ stiffness, layers = (
+ self.integrator.integrate_strain_response_on_geometry(
+ geom,
+ [0, 0, 0, 0, 0, 0],
+ integrate='modulus',
+ mesh_size=self.mesh_size,
+ layers=self.layers,
+ )
+ )
+
+ # Save layers if needed
+ if self.layers is None and layers is not None:
+ self.layers = layers
+
+ # Calculate strain plane with Newton Rhapson Iterative method
+ num_iter = 0
+ strain = np.zeros(6)
+
+ # Factorize once the stiffness matrix if using initial
+ if initial:
+ # LU factorization
+ lu, piv = lu_factor(stiffness)
+
+ first_norm = None
+
+ # Do Newton loops
+ while True:
+ # Check if number of iterations exceeds the maximum
+ if num_iter > max_iter:
+ raise StopIteration('Maximum number of iterations reached.')
+
+ # Calculate response and residuals
+ response = np.array(self.integrate_strain_profile(strain=strain))
+ residual = loads - response
+ if initial:
+ # If the initial stiffness is used, we follow a regular
+ # Newton-Raphson scheme where we calculate the strain increment
+ # from the residual and the initial tangent stiffness matrix
+
+ # Solve using the decomposed matrix
+ delta_strain = lu_solve((lu, piv), residual)
+
+ else:
+ # The the current stiffness is used, we use a secant stiffness
+ # and calculate a new total strain from which we derive the
+ # strain increment
+
+ # Calculate the current stiffness
+ stiffness, _ = (
+ self.integrator.integrate_strain_response_on_geometry(
+ geom,
+ strain,
+ integrate='modulus',
+ layers=self.layers,
+ )
+ )
+
+ # Solve using the current stiffness
+ delta_strain = np.linalg.solve(stiffness, loads) - strain
+
+ # Update the strain
+ strain += delta_strain
+
+ num_iter += 1
+
+ # Check for convergence:
+ if first_norm is None:
+ first_norm = np.dot(residual, delta_strain)
+ current_norm = first_norm
+ else:
+ current_norm = np.dot(residual, delta_strain)
+
+ if abs(current_norm / first_norm) < tol and num_iter > 1:
+ break
+
+ return strain.tolist()
diff --git a/structuralcodes/sections/section_integrators/__init__.py b/structuralcodes/sections/section_integrators/__init__.py
index 1f57ffc0..0fec7a53 100644
--- a/structuralcodes/sections/section_integrators/__init__.py
+++ b/structuralcodes/sections/section_integrators/__init__.py
@@ -6,6 +6,7 @@
from ._fiber_integrator import FiberIntegrator
from ._marin_integrator import MarinIntegrator
from ._section_integrator import SectionIntegrator
+from ._shell_integrator import ShellFiberIntegrator
__all__ = [
'integrator_factory',
@@ -13,4 +14,5 @@
'MarinIntegrator',
'SectionIntegrator',
'marin_integration',
+ 'ShellFiberIntegrator',
]
diff --git a/structuralcodes/sections/section_integrators/_shell_integrator.py b/structuralcodes/sections/section_integrators/_shell_integrator.py
new file mode 100644
index 00000000..b619bc30
--- /dev/null
+++ b/structuralcodes/sections/section_integrators/_shell_integrator.py
@@ -0,0 +1,206 @@
+"""A concrete implementation of a shell integrator using fiber
+discretization.
+"""
+
+import math
+import typing as t
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from ...geometry._shell_geometry import ShellGeometry
+from ._section_integrator import SectionIntegrator
+
+
+class ShellFiberIntegrator(SectionIntegrator):
+ """A concrete implementation of a fiber integrator for shell sections."""
+
+ def prepare_input(
+ self,
+ geo: ShellGeometry,
+ strain: ArrayLike,
+ integrate: t.Literal['stress', 'modulus'] = 'stress',
+ **kwargs,
+ ) -> t.Tuple[t.List[t.Tuple[np.ndarray, np.ndarray]], None]:
+ """Prepare general input to the integration of stress or material
+ modulus in the shell section.
+ Calculate the stress resultants or secant section stiffness based on
+ strains at discrete fibers along the shell thickness.
+
+ Arguments:
+ geo (ShellGeometry): The shell geometry object.
+ strain (ArrayLike): The strains and curvatures of the shell section
+ [eps_x, eps_y, gamma_xy, kappa_x, kappa_y, kappa_xy].
+ integrate (str): A string indicating the quantity to integrate over
+ the shell section. It can be 'stress' or 'modulus'. When
+ 'stress' is selected, the return value will be the generalised
+ stress resultants: membrane forces and bending moments Nx, Ny,
+ Nxy, Mx, My, Mxy. When 'modulus' is selected, the return value
+ will be the section stiffness matrix K (6x6), which relates
+ generalised strains to stress resultants (default is 'stress').
+
+ Keyword Arguments:
+ z_coords (ArrayLike): The z-coordinates of the layers in the shell.
+ mesh_size (float): fraction of the total shell thickness for each
+ layer ([0,1]). Default is 0.01.
+
+ Returns:
+ Tuple: (prepared_input, z_coords)
+ """
+ strain = np.atleast_1d(strain)
+
+ layers = kwargs.get('layers')
+ z_coords, dz = (None, None) if layers is None else layers
+ t_total = geo.thickness
+ if z_coords is None and dz is None:
+ mesh_size = kwargs.get('mesh_size', 0.01)
+ if not (0 < mesh_size <= 1):
+ raise ValueError('mesh_size must be [0,1].')
+ n_layers = max(1, math.ceil(1 / mesh_size))
+ dz = t_total / n_layers
+ z_coords = np.linspace(
+ -t_total / 2 + dz / 2, t_total / 2 - dz / 2, n_layers
+ )
+ material = geo.material
+ prepared_input = []
+ IA = []
+ z_list = []
+
+ # A positive in-plane strain gives tension
+ # A positive curvature gives a negative strain for a positive z-value
+ for z in z_coords:
+ fiber_strain = strain[:3] - z * strain[3:]
+ if integrate == 'stress':
+ integrand = material.constitutive_law.get_stress(fiber_strain)
+ elif integrate == 'modulus':
+ integrand = material.constitutive_law.get_secant(fiber_strain)
+ else:
+ raise ValueError(f'Unknown integrate type: {integrate}')
+
+ IA.append(integrand * dz)
+ z_list.append(z)
+
+ for r in geo.reinforcement:
+ z_r = r.z
+ material = r.material
+
+ fiber_strain = strain[:3] - z_r * strain[3:]
+ eps_sj = r.T @ fiber_strain
+
+ if integrate == 'stress':
+ sig_sj = material.constitutive_law.get_stress(eps_sj[0])
+ integrand = r.local_stress_to_global_force * sig_sj
+ elif integrate == 'modulus':
+ mod = material.constitutive_law.get_secant(eps_sj[0])
+ integrand = r.local_modulus_to_global_stiffness * mod
+ else:
+ raise ValueError(f'Unknown integrate type: {integrate}')
+
+ IA.append(integrand)
+ z_list.append(z_r)
+
+ MA = np.stack(IA, axis=0)
+ prepared_input = [(z_list, MA)]
+
+ return prepared_input, (z_coords, dz)
+
+ def integrate_stress(
+ self,
+ prepared_input: t.List[t.Tuple[np.ndarray, np.ndarray]],
+ ) -> t.Tuple[float, float, float, float, float, float]:
+ """Integrate stresses over the shell thickness.
+
+ Arguments:
+ prepared_input (List): The prepared input from .prepare_input().
+
+ Returns:
+ Tuple(float, float, float, float, float, float):
+ The stress resultants Nx, Ny, Nxy, Mx, My, Mxy
+ """
+ z, fiber_stress = prepared_input[0]
+ Nx = np.sum(fiber_stress[:, 0])
+ Ny = np.sum(fiber_stress[:, 1])
+ Nxy = np.sum(fiber_stress[:, 2])
+
+ # A positive stress at a positive z-value gives a negative moment
+ # contribution
+ Mx = -np.sum(fiber_stress[:, 0] * z)
+ My = -np.sum(fiber_stress[:, 1] * z)
+ Mxy = -np.sum(fiber_stress[:, 2] * z)
+ return Nx, Ny, Nxy, Mx, My, Mxy
+
+ def integrate_modulus(
+ self,
+ prepared_input: t.List[t.Tuple[np.ndarray, np.ndarray]],
+ ) -> NDArray:
+ """Integrate material modulus over shell thickness to obtain shell
+ section stiffness.
+
+ Arguments:
+ prepared_input (List): The prepared input from .prepare_input().
+
+ Returns:
+ NDArray: Section stiffness matrix (6x6).
+ """
+ z, MA = prepared_input[0]
+ A = np.zeros((3, 3))
+ B = np.zeros((3, 3))
+ D = np.zeros((3, 3))
+ for C_layer, z_i in zip(MA, z):
+ A += C_layer
+ B -= z_i * C_layer
+ D += z_i**2 * C_layer
+
+ return np.block([[A, B], [B, D]])
+
+ def integrate_strain_response_on_geometry(
+ self,
+ geo: ShellGeometry,
+ strain: ArrayLike,
+ integrate: t.Literal['stress', 'modulus'] = 'stress',
+ **kwargs,
+ ) -> t.Union[
+ t.Tuple[float, float, float, float, float, float], np.ndarray
+ ]:
+ """Integrate strain profile through the shell thickness.
+ This method evaluates the response of a shell section subjected to a
+ generalised strain state, either by integrating the stresses to obtain
+ resultant forces and moments, or by integrating the secant stiffness
+ to obtain the section stiffness matrix.
+
+ Arguments:
+ geo (ShellGeometry): The shell geometry including thickness and
+ material.
+ strain (ArrayLike): Generalised strain vector of the form:
+ [eps_x, eps_y, gamma_xy, kappa_x, kappa_y, kappa_xy]
+ integrate (str): Type of integration to perform. Must be one of:
+ - 'stress': returns the stress resultants (Nx, Ny, Nxy, Mx, My,
+ Mxy)
+ - 'modulus': returns the section stiffness matrix (6x6)
+
+ Returns:
+ Union[Tuple[float, float, float, float, float, float], np.ndarray]:
+ - If integrate = 'stress': returns 6 stress resultants.
+ - If integrate = 'modulus': returns the 6x6 stiffness matrix.
+
+ Example:
+ result = integrate_strain_response_on_geometry(geo, strain,
+ integrate='modulus')
+ 'result will be the 6x6 stiffness matrix of the shell.
+
+ Raises:
+ ValueError: If `integrate` is not 'stress' or 'modulus'.
+ """
+ # Prepare the general input based on the geometry and the input strains
+ prepared_input, layers = self.prepare_input(
+ geo=geo,
+ strain=strain,
+ integrate=integrate,
+ **kwargs,
+ )
+ # Return the calculated response
+ if integrate == 'stress':
+ return *self.integrate_stress(prepared_input), layers
+ if integrate == 'modulus':
+ return self.integrate_modulus(prepared_input), layers
+ raise ValueError(f'Unknown integrate type: {integrate}')
diff --git a/tests/test_geometry/test_shell.py b/tests/test_geometry/test_shell.py
new file mode 100644
index 00000000..60e0c1bf
--- /dev/null
+++ b/tests/test_geometry/test_shell.py
@@ -0,0 +1,183 @@
+"""Tests for the Geometry."""
+
+import numpy as np
+import pytest
+
+from structuralcodes.geometry import (
+ PointGeometry,
+ ShellGeometry,
+ ShellReinforcement,
+)
+from structuralcodes.materials.basic import GenericMaterial
+from structuralcodes.materials.concrete import ConcreteEC2_2004
+from structuralcodes.materials.constitutive_laws import (
+ Elastic2D,
+)
+from structuralcodes.materials.reinforcement import ReinforcementEC2_2004
+
+
+def test_shell_geometry():
+ """Test the ShellGeometry class."""
+ # Choose a constitutive law to use
+ const = Elastic2D(E=200000, nu=0.2)
+
+ # Create a material to use
+ concrete = ConcreteEC2_2004(fck=35, constitutive_law=const)
+
+ shell = ShellGeometry(
+ thickness=200, material=concrete, name='Shell', group_label='Group1'
+ )
+ assert shell.thickness == 200
+ assert isinstance(shell.material, ConcreteEC2_2004)
+ assert shell.name == 'Shell'
+ assert shell.group_label == 'Group1'
+
+
+def test_negative_thickness_raises():
+ """Test that a negative thickness raises a ValueError."""
+ material = GenericMaterial(
+ density=7850, constitutive_law=Elastic2D(E=200000, nu=0.2)
+ )
+ with pytest.raises(ValueError):
+ ShellGeometry(thickness=-200, material=material)
+
+
+def test_shell_reinforcement():
+ """Test the ShellReinforcement class."""
+ z = -60
+ n_bars = 4
+ cc_bars = 500
+ d = 16
+ material = GenericMaterial(
+ density=7850, constitutive_law=Elastic2D(E=200000, nu=0.3)
+ )
+ phi = np.pi / 4
+ shell_reinforcement = ShellReinforcement(
+ z=z,
+ n_bars=n_bars,
+ cc_bars=cc_bars,
+ diameter_bar=d,
+ material=material,
+ phi=phi,
+ name='rebars_1',
+ group_label='group_A',
+ )
+ assert shell_reinforcement.z == z
+ assert shell_reinforcement.n_bars == n_bars
+ assert shell_reinforcement.cc_bars == cc_bars
+ assert shell_reinforcement.diameter_bar == d
+ assert shell_reinforcement.material == material
+ assert shell_reinforcement.phi == phi
+ assert shell_reinforcement.name == 'rebars_1'
+ assert shell_reinforcement.group_label == 'group_A'
+
+
+def test_add_reinforcement():
+ """Test the add_reinforcement function."""
+ # Create a shell geometry
+ material = GenericMaterial(
+ density=7850, constitutive_law=Elastic2D(E=200000, nu=0.3)
+ )
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500, Es=200000, ftk=500, epsuk=3e-2
+ )
+ shell = ShellGeometry(thickness=200, material=material)
+
+ # Create a reinforcement
+ reinf_1 = ShellReinforcement(
+ z=-60,
+ n_bars=4,
+ cc_bars=500,
+ diameter_bar=16,
+ material=reinforcement,
+ phi=0,
+ )
+
+ reinf_2 = ShellReinforcement(
+ z=60,
+ n_bars=6,
+ cc_bars=600,
+ diameter_bar=20,
+ material=reinforcement,
+ phi=np.pi / 4,
+ )
+
+ # Add a single reinforcement
+ shell.add_reinforcement(reinf_1)
+ assert len(shell.reinforcement) == 1
+ assert shell.reinforcement[0] is reinf_1
+
+ # Add reinforcement as a list
+ shell.add_reinforcement([reinf_1, reinf_2])
+ assert len(shell.reinforcement) == 3
+
+
+def test_add_reinforcement_invalid_type():
+ """Test that adding a reinforcement of invalid type raises an error."""
+ # Create a shell geometry
+ material = GenericMaterial(
+ density=7850, constitutive_law=Elastic2D(E=200000, nu=0.2)
+ )
+ shell = ShellGeometry(thickness=200, material=material)
+
+ # Create a reinforcement
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500, ftk=500, Es=200000, epsuk=3e-2
+ )
+ reinf_1 = PointGeometry(np.array([2, 3]), 12, reinforcement, name='Rebar')
+ with pytest.raises(TypeError):
+ shell.add_reinforcement(reinf_1)
+
+
+@pytest.mark.parametrize(
+ 'invalid_z',
+ [
+ 95, # barely outside of shell (half_thickness=100, d/2=6-> max=94)
+ 250, # clearly outside the shell
+ ],
+)
+def test_add_reinforcement_invalid_z_value(invalid_z):
+ """Test that adding a reinforcement outside the shell raises an error."""
+ material = GenericMaterial(
+ density=7850, constitutive_law=Elastic2D(E=200000, nu=0.2)
+ )
+ shell = ShellGeometry(thickness=200, material=material)
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500, ftk=500, Es=200000, epsuk=3e-2
+ )
+ reinf = ShellReinforcement(
+ z=invalid_z,
+ n_bars=1,
+ cc_bars=100,
+ diameter_bar=12,
+ material=reinforcement,
+ phi=0,
+ )
+ with pytest.raises(ValueError):
+ shell.add_reinforcement([reinf])
+
+
+def test_repr_svg():
+ """Test the SVG representation of the shell geometry."""
+ # Create a shell geometry
+ material = GenericMaterial(
+ density=7850, constitutive_law=Elastic2D(E=200000, nu=0.2)
+ )
+ shell = ShellGeometry(thickness=200, material=material)
+
+ # Add a reinforcement
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500, ftk=500, Es=200000, epsuk=3e-2
+ )
+ reinf = ShellReinforcement(
+ z=-60,
+ n_bars=4,
+ cc_bars=500,
+ diameter_bar=16,
+ material=reinforcement,
+ phi=0,
+ )
+ shell.add_reinforcement(reinf)
+ # Get the SVG representation
+ svg = shell._repr_svg_()
+ assert ' 0)
+
+
+@pytest.mark.parametrize(
+ 'fc, eps_0, eps_u, nu, strain, expected',
+ [
+ (
+ -30.0,
+ -0.002,
+ -0.0035,
+ 0.0,
+ [-0.001, -0.001, 0.0],
+ [[22500, 0.0, 0.0], [0.0, 22500, 0.0], [0.0, 0.0, 0.5 * 22500]],
+ ),
+ (
+ -45.0,
+ -0.002,
+ -0.0035,
+ 0.2,
+ [-0.003, 0.002, 0],
+ [[15000, 0, 0], [0, 0, 0], [0, 0, 3750]],
+ ),
+ (
+ -30.0,
+ -0.002,
+ -0.0035,
+ 0.2,
+ [-0.001, 0.0, 0.0],
+ [
+ [2.31119792e04, 5.27343750e03, 0.0],
+ [5.27343750e03, 2.96223958e04, 0.0],
+ [0.0, 0.0, 1.05468750e04],
+ ],
+ ),
+ ],
+)
+def test_get_secant(fc, eps_0, eps_u, nu, strain, expected):
+ """Test the secant stiffness matrix of the concrete smeared cracking
+ material.
+ """
+ uniaxial_compression = ParabolaRectangle(fc=fc, eps_0=eps_0, eps_u=eps_u)
+ strength_reduction = GeneralVecchioCollins(c_1=0.8, c_2=100)
+ poisson_reduction = ConstantPoissonReduction(initial_nu=nu)
+ mat = ConcreteSmearedCracking(
+ uniaxial_compression=uniaxial_compression,
+ strength_reduction_lateral_cracking=strength_reduction,
+ poisson_reduction=poisson_reduction,
+ )
+ assert np.allclose(mat.get_secant(strain), expected)
+
+
@pytest.mark.parametrize(
'x, y, flag, strain, expected',
[
diff --git a/tests/test_sections/test_shell_section.py b/tests/test_sections/test_shell_section.py
new file mode 100644
index 00000000..89ecdfae
--- /dev/null
+++ b/tests/test_sections/test_shell_section.py
@@ -0,0 +1,565 @@
+"""Tests for the Shell Section."""
+
+import math
+
+import numpy as np
+import pytest
+
+from structuralcodes.geometry import (
+ RectangularGeometry,
+ add_reinforcement_line,
+)
+from structuralcodes.geometry._shell_geometry import (
+ ShellGeometry,
+ ShellReinforcement,
+)
+from structuralcodes.materials.basic import GenericMaterial
+from structuralcodes.materials.concrete import ConcreteEC2_2004
+from structuralcodes.materials.constitutive_laws import (
+ ConcreteSmearedCracking,
+ ConstantPoissonReduction,
+ Elastic2D,
+ ElasticPlastic,
+ GeneralVecchioCollins,
+ ParabolaRectangle,
+)
+from structuralcodes.materials.reinforcement import ReinforcementEC2_2004
+from structuralcodes.sections import GenericSection, ShellSection
+
+# Membrane and bending-strain parameter ranges
+eps_x = np.linspace(0.0, 1.0e-3, 2)
+eps_y = np.linspace(0.0, 1.0e-3, 2)
+eps_xy = np.linspace(0.0, 1.0e-3, 2)
+chi_x = np.linspace(0.0, 1.0e-6, 2)
+chi_y = np.linspace(0.0, 1.0e-6, 2)
+chi_xy = np.linspace(0.0, 1.0e-6, 2)
+
+
+@pytest.mark.parametrize('chi_xy', chi_xy)
+@pytest.mark.parametrize('chi_y', chi_y)
+@pytest.mark.parametrize('chi_x', chi_x)
+@pytest.mark.parametrize('eps_xy', eps_xy)
+@pytest.mark.parametrize('eps_y', eps_y)
+@pytest.mark.parametrize('eps_x', eps_x)
+@pytest.mark.parametrize(
+ 'Ec, nu, thickness',
+ [(30000, 0.20, 200), (30000, 0.20, 600)],
+)
+def test_integrate_strain_profile(
+ Ec,
+ nu,
+ thickness,
+ eps_x,
+ eps_y,
+ eps_xy,
+ chi_x,
+ chi_y,
+ chi_xy,
+):
+ """Elastic plate: strains → stress-resultants."""
+ constitutive_law = Elastic2D(Ec, nu)
+ material = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell = ShellSection(ShellGeometry(thickness, material), mesh_size=1e-3)
+
+ # Numerically integrated forces/moments
+ R = shell.section_calculator.integrate_strain_profile(
+ np.array([eps_x, eps_y, eps_xy, chi_x, chi_y, chi_xy])
+ )
+
+ # Analytical forces/moments
+ A_membrane = Ec * thickness / (1 - nu**2) # Coefficient in Hooke's law
+ A_bending = Ec * thickness**3 / 12 / (1 - nu**2) # Plate bending stiffness
+ A_membrane_shear = A_membrane * (1 - nu) / 2
+ A_bending_shear = A_bending * (1 - nu) / 2
+
+ Nx = A_membrane * (eps_x + nu * eps_y)
+ Ny = A_membrane * (eps_y + nu * eps_x)
+ Nxy = A_membrane_shear * eps_xy
+ Mx = A_bending * (chi_x + nu * chi_y)
+ My = A_bending * (chi_y + nu * chi_x)
+ Mxy = A_bending_shear * chi_xy
+
+ expected = np.array([Nx, Ny, Nxy, Mx, My, Mxy])
+
+ assert np.allclose(R, expected)
+
+
+@pytest.mark.parametrize('nu', ((0, 0.2)))
+def test_integrate_strain_profile_stiffness_matrix(nu):
+ """Elastic plate: stiffness matrix."""
+ E = 30000
+ t = 200
+ constitutive_law = Elastic2D(E, nu)
+ material = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell = ShellSection(ShellGeometry(t, material=material), mesh_size=1e-3)
+ K = shell.section_calculator.integrate_strain_profile(
+ np.array([1e-3] * 3 + [1e-6] * 3), integrate='modulus'
+ )
+
+ A11 = E * t / (1 - nu**2)
+ A12 = nu * A11
+ A33 = E * t / (2 * (1 + nu))
+ A44 = E * t**3 / 12 / (1 - nu**2)
+ A14 = nu * A44
+ A55 = E * t**3 / 12 / (2 * (1 + nu))
+
+ A = np.array(
+ [
+ [A11, A12, 0, 0, 0, 0],
+ [A12, A11, 0, 0, 0, 0],
+ [0, 0, A33, 0, 0, 0],
+ [0, 0, 0, A44, A14, 0],
+ [0, 0, 0, A14, A44, 0],
+ [0, 0, 0, 0, 0, A55],
+ ]
+ )
+
+ assert np.allclose(K, A)
+
+
+def test_wrong_integrator():
+ """Unknown keyword must raise ValueError."""
+ constitutive_law = Elastic2D(30000, 0.20)
+ material = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell = ShellSection(ShellGeometry(200, material=material))
+ with pytest.raises(ValueError):
+ shell.section_calculator.integrate_strain_profile(
+ np.zeros(6), integrate='tangent'
+ )
+
+
+# Loads for reverse strain-solution test
+nx = np.linspace(-1e5, 1e5, 2)
+ny = np.linspace(-1e5, 1e5, 2)
+nxy = np.linspace(-1e5, 1e5, 2)
+mx = np.linspace(-1e8, 1e8, 2)
+my = np.linspace(-1e8, 1e8, 2)
+mxy = np.linspace(-1e8, 1e8, 2)
+
+
+@pytest.mark.parametrize('mxy', mxy)
+@pytest.mark.parametrize('my', my)
+@pytest.mark.parametrize('mx', mx)
+@pytest.mark.parametrize('nxy', nxy)
+@pytest.mark.parametrize('ny', ny)
+@pytest.mark.parametrize('nx', nx)
+@pytest.mark.parametrize('Ec, nu, t', [(30000, 0.20, 200), (30000, 0.20, 600)])
+def test_elastic_strain_profile(Ec, nu, t, nx, ny, nxy, mx, my, mxy):
+ """Loads → strains for an isotropic plate."""
+ constitutive_law = Elastic2D(Ec, nu)
+ material = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell = ShellSection(ShellGeometry(t, material=material), mesh_size=1e-3)
+ eps = shell.section_calculator.calculate_strain_profile(
+ nx, ny, nxy, mx, my, mxy
+ )
+
+ sig_x = nx / t
+ sig_y = ny / t
+ tau_xy = nxy / t
+ eps_x = (sig_x - nu * sig_y) / Ec
+ eps_y = (sig_y - nu * sig_x) / Ec
+ eps_xy = 2 * (1 + nu) * tau_xy / Ec
+ D = Ec * t**3 / 12
+ chi_x = (mx - nu * my) / D
+ chi_y = (my - nu * mx) / D
+ chi_xy = 2 * (1 + nu) * mxy / D
+ expected = np.array([eps_x, eps_y, eps_xy, chi_x, chi_y, chi_xy])
+
+ assert np.allclose(eps, expected)
+
+
+def test_default_equals_explicit_mesh_size():
+ """Default mesh_size (0.01) equals explicit 0.01."""
+ t = 200
+ constitutive_law = Elastic2D(30000, 0.20)
+ material = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell_0 = ShellSection(ShellGeometry(t, material=material))
+ K0 = shell_0.section_calculator.integrate_strain_profile(
+ np.zeros(6), integrate='modulus'
+ )
+ shell_1 = ShellSection(ShellGeometry(t, material=material), mesh_size=0.01)
+ K1 = shell_1.section_calculator.integrate_strain_profile(
+ np.zeros(6), integrate='modulus'
+ )
+ assert np.allclose(K0, K1)
+
+
+@pytest.mark.parametrize('invalid', [-0.2, 0.0, 1.5])
+def test_invalid_mesh_size_raises(invalid):
+ """mesh_size outside (0,1] → ValueError."""
+ t = 200
+ constitutive_law = Elastic2D(30000, 0.20)
+ material = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell = ShellSection(ShellGeometry(t, material), mesh_size=invalid)
+ with pytest.raises(ValueError):
+ shell.section_calculator.integrate_strain_profile(
+ np.zeros(6), integrate='stress'
+ )
+
+
+@pytest.mark.parametrize(
+ 'nx,nxy,expected',
+ [
+ (
+ 0,
+ 500,
+ np.array([1.43163e-3, 1.919842e-3, 3.466593e-3, 0, 0, 0]),
+ ),
+ (
+ -500,
+ 500,
+ np.array([0.618194e-3, 1.476430e-3, 2.065709e-3, 0, 0, 0]),
+ ),
+ (
+ 500,
+ 500,
+ np.array([2.446863e-3, 2.283832e-3, 4.894237e-3, 0, 0, 0]),
+ ),
+ ],
+)
+def test_parabola_zero_initial_nu(nx, nxy, expected):
+ """Test parabola rectangle with nu = 0."""
+ uniaxial_compression = ParabolaRectangle(fc=45)
+ strength_reduction = GeneralVecchioCollins(c_1=0.8, c_2=100)
+ poisson_reduction = ConstantPoissonReduction(initial_nu=0)
+ parabola_rectangle = ConcreteSmearedCracking(
+ uniaxial_compression=uniaxial_compression,
+ strength_reduction_lateral_cracking=strength_reduction,
+ poisson_reduction=poisson_reduction,
+ )
+ elastic_plastic = ElasticPlastic(200000, 500)
+ concrete = GenericMaterial(
+ density=2500, constitutive_law=parabola_rectangle
+ )
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500,
+ Es=200000,
+ ftk=500,
+ epsuk=3e-2,
+ constitutive_law=elastic_plastic,
+ )
+ geo = ShellGeometry(350, concrete)
+ Asx1 = ShellReinforcement(-132, 1, 200, 16, reinforcement, 0)
+ Asx2 = ShellReinforcement(132, 1, 200, 16, reinforcement, 0)
+ Asy1 = ShellReinforcement(-118, 1, 200, 12, reinforcement, np.pi / 2)
+ Asy2 = ShellReinforcement(118, 1, 200, 12, reinforcement, np.pi / 2)
+ geo.add_reinforcement([Asx1, Asx2, Asy1, Asy2])
+ section = ShellSection(geo, mesh_size=0.5)
+ calculator = section.section_calculator
+ strain = calculator.calculate_strain_profile(
+ nx, 0, nxy, 0, 0, 0, max_iter=150
+ )
+
+ assert np.allclose(strain, expected)
+
+
+@pytest.mark.parametrize(
+ 'nx,nxy,expected',
+ [
+ (
+ 0,
+ 500,
+ np.array([1.431675e-3, 1.919788e-3, 3.466592e-3, 0, 0, 0]),
+ ),
+ (
+ -500,
+ 500,
+ np.array([0.618204e-3, 1.476423e-3, 2.065719e-3, 0, 0, 0]),
+ ),
+ (
+ 500,
+ 500,
+ np.array([2.446899e-3, 2.283762e-3, 4.894200e-3, 0, 0, 0]),
+ ),
+ ],
+)
+def test_parabola_initial_nu(nx, nxy, expected):
+ """Test parabola rectangle with nu = 0.2."""
+ uniaxial_compression = ParabolaRectangle(fc=45)
+ strength_reduction = GeneralVecchioCollins(c_1=0.8, c_2=100)
+ poisson_reduction = ConstantPoissonReduction(initial_nu=0.2)
+ parabola_rectangle = ConcreteSmearedCracking(
+ uniaxial_compression=uniaxial_compression,
+ strength_reduction_lateral_cracking=strength_reduction,
+ poisson_reduction=poisson_reduction,
+ )
+ elastic_plastic = ElasticPlastic(200000, 500)
+ concrete = GenericMaterial(
+ density=2500, constitutive_law=parabola_rectangle
+ )
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500,
+ Es=200000,
+ ftk=500,
+ epsuk=3e-2,
+ constitutive_law=elastic_plastic,
+ )
+ geo = ShellGeometry(350, concrete)
+ Asx1 = ShellReinforcement(-132, 1, 200, 16, reinforcement, 0)
+ Asx2 = ShellReinforcement(132, 1, 200, 16, reinforcement, 0)
+ Asy1 = ShellReinforcement(-118, 1, 200, 12, reinforcement, np.pi / 2)
+ Asy2 = ShellReinforcement(118, 1, 200, 12, reinforcement, np.pi / 2)
+ geo.add_reinforcement([Asx1, Asx2, Asy1, Asy2])
+ section = ShellSection(geo, mesh_size=0.5)
+ strain = section.section_calculator.calculate_strain_profile(
+ nx,
+ 0,
+ nxy,
+ 0,
+ 0,
+ 0,
+ max_iter=200,
+ )
+
+ assert np.allclose(strain, expected)
+
+
+def test_exceed_max_iterations():
+ """Test that the maximum number of iterations is exceeded."""
+ uniaxial_compression = ParabolaRectangle(fc=45)
+ strength_reduction = GeneralVecchioCollins(c_1=0.8, c_2=100)
+ poisson_reduction = ConstantPoissonReduction(initial_nu=0)
+ parabola_rectangle = ConcreteSmearedCracking(
+ uniaxial_compression=uniaxial_compression,
+ strength_reduction_lateral_cracking=strength_reduction,
+ poisson_reduction=poisson_reduction,
+ )
+ elastic_plastic = ElasticPlastic(200000, 500)
+ concrete = GenericMaterial(
+ density=2500, constitutive_law=parabola_rectangle
+ )
+ reinforcement = ReinforcementEC2_2004(
+ fyk=500,
+ Es=200000,
+ ftk=500,
+ epsuk=3e-2,
+ constitutive_law=elastic_plastic,
+ )
+ geo = ShellGeometry(350, concrete)
+ Asx1 = ShellReinforcement(-132, 1, 200, 16, reinforcement, 0)
+ Asx2 = ShellReinforcement(132, 1, 200, 16, reinforcement, 0)
+ Asy1 = ShellReinforcement(-118, 1, 200, 12, reinforcement, np.pi / 2)
+ Asy2 = ShellReinforcement(118, 1, 200, 12, reinforcement, np.pi / 2)
+ geo.add_reinforcement([Asx1, Asx2, Asy1, Asy2])
+ section = ShellSection(geo)
+
+ with pytest.raises(
+ StopIteration, match='Maximum number of iterations reached'
+ ):
+ section.section_calculator.calculate_strain_profile(
+ -1000,
+ 0,
+ 1000,
+ 0,
+ 0,
+ 0,
+ )
+
+
+@pytest.mark.parametrize(
+ 'axial_force, with_reinforcement',
+ (
+ (-1e6, False),
+ (-1e6, True),
+ (1e5, True),
+ ),
+)
+def test_compare_uniaxial_with_generic_section_reinforcement( # noqa: PLR0915
+ axial_force: float, with_reinforcement: bool
+):
+ """Compare the uniaxial response of the shell section with the generic
+ section, with reinforcement.
+ """
+ # Arrange
+ # Geometry
+ width = 1000
+ height = 350
+ cover = 35
+ diameter = 16
+ spacing = 200
+
+ # Material parameters
+ fc = 45 # Concrete compressive strength
+ nu = 0.0 # Poisson's ratio
+ fyk = 500 # Steel yield strength
+ Es = 200000 # Steel modulus of elasticity
+
+ # Create a GenericSection
+ reinforcement = ReinforcementEC2_2004(
+ fyk=fyk,
+ Es=Es,
+ epsuk=3e-2,
+ ftk=fyk,
+ constitutive_law='elasticperfectlyplastic',
+ )
+ parabola_rectangle = ParabolaRectangle(fc=fc)
+ concrete_for_generic = ConcreteEC2_2004(
+ fck=fc, constitutive_law=parabola_rectangle
+ )
+ generic_geo = RectangularGeometry(
+ height=height, width=width, material=concrete_for_generic
+ )
+
+ if with_reinforcement:
+ # Bottom reinforcement
+ generic_geo = add_reinforcement_line(
+ generic_geo,
+ (
+ -width / 2 + cover + diameter / 2,
+ -height / 2 + cover + diameter / 2,
+ ),
+ (
+ width / 2 - cover - diameter / 2,
+ -height / 2 + cover + diameter / 2,
+ ),
+ diameter=diameter,
+ material=reinforcement,
+ s=spacing,
+ )
+
+ # Top reinforcement
+ generic_geo = add_reinforcement_line(
+ generic_geo,
+ (
+ -width / 2 + cover + diameter / 2,
+ height / 2 - cover - diameter / 2,
+ ),
+ (
+ width / 2 - cover - diameter / 2,
+ height / 2 - cover - diameter / 2,
+ ),
+ diameter=diameter,
+ material=reinforcement,
+ s=spacing,
+ )
+
+ generic_sec = GenericSection(geometry=generic_geo, integrator='fiber')
+
+ # Create a ShellSection
+ uniaxial_compression = ParabolaRectangle(fc=fc)
+ strength_reduction = GeneralVecchioCollins(c_1=0.8, c_2=100)
+ poisson_reduction = ConstantPoissonReduction(initial_nu=nu)
+ smeared_cracking = ConcreteSmearedCracking(
+ uniaxial_compression=uniaxial_compression,
+ strength_reduction_lateral_cracking=strength_reduction,
+ poisson_reduction=poisson_reduction,
+ )
+ concrete_for_shell = GenericMaterial(
+ density=2500, constitutive_law=smeared_cracking
+ )
+ shell_geo = ShellGeometry(material=concrete_for_shell, thickness=height)
+
+ if with_reinforcement:
+ # Create reinforcement for the shell section
+ z_btm = -height / 2 + cover + diameter / 2 # -132
+ z_top = height / 2 - cover - diameter / 2 # 132
+
+ # Bottom reinforcement
+ shell_rein_btm = ShellReinforcement(
+ z_btm, 1, spacing, diameter, reinforcement, 0
+ )
+
+ # Top reinforcement
+ shell_rein_top = ShellReinforcement(
+ z_top, 1, spacing, diameter, reinforcement, 0
+ )
+
+ # Add reinforcement to the shell geometry
+ shell_geo.add_reinforcement([shell_rein_btm, shell_rein_top])
+
+ shell_sec = ShellSection(geometry=shell_geo, mesh_size=0.5)
+
+ # Calculate "exact" solution
+ if axial_force > 0:
+ exact_longitudinal_strain = axial_force / (
+ generic_sec.gross_properties.area_reinforcement * Es
+ )
+ else:
+ exact_longitudinal_strain = 0
+ num_iter = 0
+ area = width * height
+ while True:
+ num_iter += 1
+ if num_iter >= 40:
+ break
+ exact_longitudinal_strain = axial_force / (
+ area
+ * concrete_for_generic.constitutive_law.get_secant(
+ exact_longitudinal_strain
+ )
+ + generic_sec.gross_properties.area_reinforcement * Es
+ )
+
+ # Act
+ # Calculate strains
+ generic_strain = generic_sec.section_calculator.calculate_strain_profile(
+ n=axial_force, my=0.0, mz=0.0
+ )
+ shell_strain = shell_sec.section_calculator.calculate_strain_profile(
+ nx=axial_force / width,
+ ny=0.0,
+ nxy=0.0,
+ mx=0.0,
+ my=0.0,
+ mxy=0.0,
+ )
+
+ # Compare with exact solution
+ # Note that the tolerance in isclose is set rather loose because the
+ # iterative solver of the shell section converges slower than the generic
+ # section.
+ assert math.isclose(
+ shell_strain[0],
+ exact_longitudinal_strain,
+ rel_tol=1e-4,
+ )
+
+ # Compare GenericSection and ShellSection
+ assert math.isclose(
+ generic_strain[0],
+ shell_strain[0],
+ rel_tol=1e-4,
+ )
+
+
+@pytest.mark.parametrize('nu', [0.0, 0.2])
+@pytest.mark.parametrize(
+ 'strain',
+ [
+ (-2e-3, 0.0, 0.0),
+ (-2e-3, -2e-3, 0.0),
+ (0.0, 0.0, 0.75e-3),
+ (-2e-3, 0.0, 0.75e-3),
+ (-2e-3, -2e-3, 0.75e-3),
+ ],
+)
+def test_compare_constitutive_law_and_section(strain, nu):
+ """Compare the stress calculated with a constitutive law with the stress
+ resultant from the section.
+ """
+ # Arrange
+ fck = 45
+ thickness = 450
+ uniaxial_compression = ParabolaRectangle(fc=fck)
+ strength_reduction = GeneralVecchioCollins(c_1=0.8, c_2=100)
+ poisson_reduction = ConstantPoissonReduction(initial_nu=nu)
+ constitutive_law = ConcreteSmearedCracking(
+ uniaxial_compression=uniaxial_compression,
+ strength_reduction_lateral_cracking=strength_reduction,
+ poisson_reduction=poisson_reduction,
+ )
+ concrete = GenericMaterial(density=2500, constitutive_law=constitutive_law)
+ shell_geometry = ShellGeometry(thickness=thickness, material=concrete)
+ shell_section = ShellSection(geometry=shell_geometry)
+
+ # Act
+ stress = constitutive_law.get_stress(eps=strain)
+ stress_resultant = (
+ shell_section.section_calculator.integrate_strain_profile(
+ strain=[*strain, 0.0, 0.0, 0.0],
+ integrate='stress',
+ )
+ )
+
+ # Assert
+ assert np.allclose(stress * thickness, stress_resultant[:3])