"""Beam2 — 3D 2-node linear beam.
femorph-solver implements the slender-beam (Euler–Bernoulli) form of Beam2:
axial, torsional, and two independent bending planes combine into a
12×12 local stiffness matrix rotated into global coordinates. This is
exactly the limit Beam2 converges to when the Timoshenko shear
parameter Φ → 0, i.e. for slender sections where bending dominates over
shear deformation. A future revision can introduce Φ from an explicit
shear area or section integration.
Local DOF ordering at each node:
(u, v, w, θx, θy, θz)
so global assembly fills rows/cols 0..5 (node i) and 6..11 (node j) in
the canonical ``(UX, UY, UZ, ROTX, ROTY, ROTZ)`` DOF order.
Real constants
--------------
real[0]: AREA — cross-sectional area A
real[1]: IZZ — bending inertia about local z (resists y-deflection)
real[2]: IYY — bending inertia about local y (resists z-deflection)
real[3]: J — Saint-Venant torsion constant
Material
--------
EX : Young's modulus
PRXY : Poisson's ratio (→ G = E/(2(1+ν)) for torsion)
DENS : mass density (required for M_e)
Orientation
-----------
The local x-axis points from node I → node J. Without an explicit
orientation node, the local y-axis is chosen so that local z stays
"up" (world +Z) whenever the beam is not parallel to the global Z-axis,
falling back to world +Y for vertical beams. For axisymmetric sections
(``IYY == IZZ``) the orientation is immaterial.
Stiffness (slender-beam limit)
------------------------------
K_local = [ K_axial · · · ]
[ · K_bend_xy · · ]
[ · · K_bend_xz · ]
[ · · · K_torsion ]
with the familiar Hermite-cubic 4×4 bending blocks. A block-diagonal
direction-cosine matrix ``Λ`` rotates it: ``K_global = Λᵀ K_local Λ``.
Mass (consistent, Hermite cubic)
--------------------------------
Axial / torsional: ``ρ A L / 6 · [[2,1],[1,2]]`` (resp. ``ρ I_p L / 6``
with ``I_p ≈ I_y + I_z``). Bending uses the standard 4×4 block
``ρ A L / 420 · (…)`` in both planes.
Lumped mass puts ``ρ A L / 2`` on each translational DOF and zeros the
rotational ones — a simple diagonal approximation useful for explicit
dynamics but less accurate for modal work than the consistent form.
References
----------
* **Euler--Bernoulli beam** (slender-beam limit, shear neglected):
Timoshenko, S.P., *Vibration Problems in Engineering*, 4th ed.,
Wiley, 1974, §5.3. Beam kinematics and derivation of the
12-DOF stiffness / consistent-mass blocks: Cook, Malkus, Plesha,
Witt, *Concepts and Applications of Finite Element Analysis*,
4th ed., Wiley, 2002, §2.4 (axial), §2.5 (bending Hermite
cubics), §2.6 (torsion).
* Hermite cubic shape functions for transverse displacement and
slope continuity at the nodes: Zienkiewicz, O.C. and Taylor, R.L.
*The Finite Element Method: Its Basis and Fundamentals*, 7th
ed., Butterworth-Heinemann, 2013, §2.5.1 eqs. (2.26)–(2.27).
* Consistent mass matrix from Hermite cubics (coefficients
``ρAL/420`` on the 4×4 bending block and ``ρAL/6`` on the
axial / torsional block): Cook et al. Table 16.3-1.
* Timoshenko-beam form (non-zero shear, Φ = 12EI/(κAGL²) — future
revision) and comparison with the Euler-Bernoulli limit:
Bathe, K.J., *Finite Element Procedures*, 2nd ed., Prentice Hall,
2014, §5.4.
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._beam_fibers import FiberGrid
from femorph_solver.elements._registry import register
from femorph_solver.materials.plasticity import (
ElastoplasticMaterial,
IntegrationPointState,
stress_update_uniaxial,
)
# 2-point Gauss-Legendre over [0, L] — natural points and weights for
# integrating Hermite-cubic curvatures (linear in s) under quadratic
# integrands. Exact for the elastic K_bb assembly and the standard
# choice for fibre-section integration in distributed-plasticity beam-
# column elements (Spacone, Filippou, Taucer 1996).
_AXIAL_GAUSS_XI: tuple[float, float] = (-1.0 / np.sqrt(3.0), +1.0 / np.sqrt(3.0))
_AXIAL_GAUSS_WT: tuple[float, float] = (1.0, 1.0)
def _beam_frame(coords: np.ndarray) -> tuple[float, np.ndarray]:
"""Return ``(L, R)`` — beam length and 3×3 direction-cosine matrix."""
d = coords[1] - coords[0]
L = float(np.linalg.norm(d))
if L == 0.0:
raise ValueError("BEAM2: coincident nodes, element length is zero")
ex = d / L
world_z = np.array([0.0, 0.0, 1.0])
if abs(np.dot(ex, world_z)) > 0.99:
ref = np.array([0.0, 1.0, 0.0])
else:
ref = world_z
ey = np.cross(ref, ex)
ey /= np.linalg.norm(ey)
ez = np.cross(ex, ey)
R = np.stack([ex, ey, ez])
return L, R
def _transform_block(R: np.ndarray) -> np.ndarray:
"""Return the 12×12 block-diagonal transformation matrix ``Λ``."""
T = np.zeros((12, 12), dtype=np.float64)
for i in range(4):
T[3 * i : 3 * i + 3, 3 * i : 3 * i + 3] = R
return T
def _k_local(L: float, E: float, G: float, A: float, Iy: float, Iz: float, J: float) -> np.ndarray:
k = np.zeros((12, 12), dtype=np.float64)
# Axial: u_i (0), u_j (6)
ka = E * A / L
k[0, 0] = k[6, 6] = ka
k[0, 6] = k[6, 0] = -ka
# Torsion: θx_i (3), θx_j (9)
kt = G * J / L
k[3, 3] = k[9, 9] = kt
k[3, 9] = k[9, 3] = -kt
# Bending in local x-y plane (v, θz) — uses I_z
# DOFs: v_i (1), θz_i (5), v_j (7), θz_j (11)
a = E * Iz / L**3
b = E * Iz / L**2
c = E * Iz / L
dofs_z = [1, 5, 7, 11]
kb = np.array(
[
[12 * a, 6 * b, -12 * a, 6 * b],
[6 * b, 4 * c, -6 * b, 2 * c],
[-12 * a, -6 * b, 12 * a, -6 * b],
[6 * b, 2 * c, -6 * b, 4 * c],
]
)
for ii, I in enumerate(dofs_z):
for jj, J_ in enumerate(dofs_z):
k[I, J_] += kb[ii, jj]
# Bending in local x-z plane (w, θy) — uses I_y, sign flips from dw/dx = -θy
# DOFs: w_i (2), θy_i (4), w_j (8), θy_j (10)
a = E * Iy / L**3
b = E * Iy / L**2
c = E * Iy / L
dofs_y = [2, 4, 8, 10]
kb = np.array(
[
[12 * a, -6 * b, -12 * a, -6 * b],
[-6 * b, 4 * c, 6 * b, 2 * c],
[-12 * a, 6 * b, 12 * a, 6 * b],
[-6 * b, 2 * c, 6 * b, 4 * c],
]
)
for ii, I in enumerate(dofs_y):
for jj, J_ in enumerate(dofs_y):
k[I, J_] += kb[ii, jj]
return k
def _m_local_consistent(L: float, rho: float, A: float, Iy: float, Iz: float) -> np.ndarray:
m = np.zeros((12, 12), dtype=np.float64)
mass = rho * A * L
# Axial
ka = mass / 6.0
m[0, 0] = m[6, 6] = 2 * ka
m[0, 6] = m[6, 0] = ka
# Torsional rotary inertia: I_p ≈ I_y + I_z per unit length
kt = rho * (Iy + Iz) * L / 6.0
m[3, 3] = m[9, 9] = 2 * kt
m[3, 9] = m[9, 3] = kt
# Bending x-y (v, θz)
base = mass / 420.0
dofs_z = [1, 5, 7, 11]
mb = base * np.array(
[
[156, 22 * L, 54, -13 * L],
[22 * L, 4 * L**2, 13 * L, -3 * L**2],
[54, 13 * L, 156, -22 * L],
[-13 * L, -3 * L**2, -22 * L, 4 * L**2],
]
)
for ii, I in enumerate(dofs_z):
for jj, J_ in enumerate(dofs_z):
m[I, J_] += mb[ii, jj]
# Bending x-z (w, θy) — sign flips on rotation coupling
dofs_y = [2, 4, 8, 10]
mb = base * np.array(
[
[156, -22 * L, 54, 13 * L],
[-22 * L, 4 * L**2, -13 * L, -3 * L**2],
[54, -13 * L, 156, 22 * L],
[13 * L, -3 * L**2, 22 * L, 4 * L**2],
]
)
for ii, I in enumerate(dofs_y):
for jj, J_ in enumerate(dofs_y):
m[I, J_] += mb[ii, jj]
return m
def _m_local_lumped(L: float, rho: float, A: float) -> np.ndarray:
m = np.zeros((12, 12), dtype=np.float64)
half = 0.5 * rho * A * L
for i in (0, 1, 2, 6, 7, 8):
m[i, i] = half
return m
[docs]
@register
class Beam2(ElementBase):
name = "BEAM2"
n_nodes = 2
n_dof_per_node = 6
dof_labels = ("UX", "UY", "UZ", "ROTX", "ROTY", "ROTZ")
vtk_cell_type = 3 # VTK_LINE
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
r = np.asarray(real, dtype=np.float64)
if r.size < 4:
raise ValueError(f"Beam2 requires REAL[0..3] = AREA, IZZ, IYY, J; got {r.size} values")
E = float(material["EX"])
nu = float(material.get("PRXY", 0.3))
G = E / (2.0 * (1.0 + nu))
A = float(r[0])
Iz = float(r[1])
Iy = float(r[2])
J = float(r[3])
k_loc = _k_local(L, E, G, A, Iy, Iz, J)
T = _transform_block(R)
return T.T @ k_loc @ T
[docs]
@staticmethod
def strain_batch(
coords: np.ndarray,
u_e: np.ndarray,
material: dict[str, float] | None = None,
) -> np.ndarray | None:
"""Per-element axial-only 3D Voigt strain at the two end nodes.
Centroid (neutral-axis) recovery only - bending fibre stress
is a per-section follow-up that needs ``y_f`` / ``z_f``
positions from the registered ``BeamSectionProperties``. The
recovered strain corresponds to the pure-uniaxial centroid
stress state ``sigma = E * (du_axial / L) * (e^e)`` where
``e`` is the element-axis unit vector and
``du_axial = (u_J - u_I) . e``.
Mirrors :meth:`Truss2.strain_batch` so
:func:`compute_nodal_stress` recovers the right global
Cartesian stress projection via the standard isotropic
``C @ eps`` contraction. Voigt order is
``[exx, eyy, ezz, gxy, gyz, gxz]`` (engineering shears).
"""
coords = np.asarray(coords, dtype=np.float64)
u_e = np.asarray(u_e, dtype=np.float64)
if coords.ndim != 3 or coords.shape[1] != 2 or coords.shape[2] != 3:
return None
# u_e is (n_e, 12) for BEAM2 (6 DOFs per node).
if u_e.ndim != 2 or u_e.shape[1] != 12 or u_e.shape[0] != coords.shape[0]:
return None
n_e = coords.shape[0]
d = coords[:, 1, :] - coords[:, 0, :]
L = np.linalg.norm(d, axis=1)
if (L <= 0.0).any():
return None
e_hat = d / L[:, None]
# Axial component of the displacement difference (ignore the
# rotational DOFs at indices 3:6 and 9:12).
u_I = u_e[:, 0:3]
u_J = u_e[:, 6:9]
du = u_J - u_I
eps_a = np.einsum("ij,ij->i", du, e_hat) / L
# 3D Voigt strain matching uniaxial stress along e_hat.
nu = float(material.get("PRXY", 0.0)) if material else 0.0
eps_v = np.empty((n_e, 6), dtype=np.float64)
ex, ey, ez = e_hat[:, 0], e_hat[:, 1], e_hat[:, 2]
eps_v[:, 0] = eps_a * ((1.0 + nu) * ex * ex - nu)
eps_v[:, 1] = eps_a * ((1.0 + nu) * ey * ey - nu)
eps_v[:, 2] = eps_a * ((1.0 + nu) * ez * ez - nu)
eps_v[:, 3] = 2.0 * eps_a * (1.0 + nu) * ex * ey
eps_v[:, 4] = 2.0 * eps_a * (1.0 + nu) * ey * ez
eps_v[:, 5] = 2.0 * eps_a * (1.0 + nu) * ex * ez
return np.repeat(eps_v[:, None, :], 2, axis=1)
[docs]
@staticmethod
def distributed_load(
coords: np.ndarray,
real: np.ndarray, # noqa: ARG004 — accepted for kernel-hook signature uniformity
material: dict[str, float], # noqa: ARG004
*,
face: int,
value: float,
) -> np.ndarray:
"""Consistent nodal-force vector for a uniform distributed load.
Implements the Hermite-cubic equivalent nodal forces for a
uniform line load ``q = value`` (force per unit length) applied
along the element axis, perpendicular to the requested local
face. Standard textbook result (Cook §5.7 Table 5.7-1) on a
Hermite-cubic Euler-Bernoulli beam:
F_v_each_end = q · L / 2
M_each_end = ± q · L² / 12 (signs flip end-to-end)
Face conventions (BEAM188 / BEAM2 local frame, X = element axis):
* face 1 → +Y (load acts in local +Y direction)
* face 2 → +Z
* face 3 → -Y
* face 4 → -Z
Returns a 12-element local-frame load vector laid out as
``[ux_I, uy_I, uz_I, θx_I, θy_I, θz_I, ux_J, uy_J, uz_J,
θx_J, θy_J, θz_J]`` already transformed to the **global** frame
via the same rotation matrix ``R`` :func:`_k_local` uses, so the
caller can scatter it directly into ``F``.
References
----------
* Cook et al., *Concepts and Applications of Finite Element
Analysis*, 4th ed., Table 5.7-1 (consistent equivalent nodal
forces for uniform / triangular loads on Hermite beams).
* Logan, D.L., *A First Course in the Finite Element Method*,
5th ed., Cengage, 2012, §5.6 (work-equivalent loads).
"""
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
q = float(value)
f_local = np.zeros(12, dtype=np.float64)
# Translation amplitudes are q·L/2 at each end; rotation (moment)
# amplitudes are ±q·L²/12 with the sign flipping between I and J.
force_each = q * L / 2.0
moment_each = q * L * L / 12.0
if face == 1 or face == 3:
sign = 1.0 if face == 1 else -1.0
# Transverse load along local Y → couples to v (uy) and θz.
f_local[1] = sign * force_each
f_local[7] = sign * force_each
f_local[5] = sign * moment_each
f_local[11] = -sign * moment_each
elif face == 2 or face == 4:
sign = 1.0 if face == 2 else -1.0
# Transverse load along local Z → couples to w (uz) and θy.
# Sign of the moment about Y flips relative to Z because of
# the right-hand-rule convention used by ``_k_local``'s
# ``[w, θy]`` bending block (note the ``-6b`` cross terms).
f_local[2] = sign * force_each
f_local[8] = sign * force_each
f_local[4] = -sign * moment_each
f_local[10] = sign * moment_each
else: # pragma: no cover - guarded at the Model._sfbeam_impl boundary
raise NotImplementedError(f"BEAM2.distributed_load: face={face} not supported")
T = _transform_block(R)
return T.T @ f_local
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
r = np.asarray(real, dtype=np.float64)
if r.size < 4:
raise ValueError(f"Beam2 requires REAL[0..3] = AREA, IZZ, IYY, J; got {r.size} values")
rho = float(material.get("DENS", 0.0))
A = float(r[0])
Iz = float(r[1])
Iy = float(r[2])
if lumped:
m_loc = _m_local_lumped(L, rho, A)
else:
m_loc = _m_local_consistent(L, rho, A, Iy, Iz)
T = _transform_block(R)
return T.T @ m_loc @ T
[docs]
@staticmethod
def fiber_stress(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
u_e: np.ndarray,
*,
fiber_y: float,
fiber_z: float = 0.0,
station_s: float = 0.0,
) -> float:
"""Axial stress ``σ_xx`` at a section fiber along a Beam2 element.
Combines axial, strong-axis bending, and weak-axis bending::
σ_xx(s, y, z) = N(s)/A − M_z(s)·y/Iz + M_y(s)·z/Iy
where ``s`` is the natural station along the element axis
(``s = 0`` at node I, ``s = L`` at node J), ``(y, z)`` are
the fiber's location relative to the section centroid in the
beam's local frame, and ``(N, M_z, M_y)`` are the internal
forces recovered from the element's nodal displacement
vector ``u_e`` via the same shape functions used in
:meth:`Beam2.ke`:
* Linear axial: ``N(s) = E·A · (u_J − u_I)/L`` (constant)
* Hermite-cubic bending about local z (``v(s)``):
``M_z(s) = E·Iz · v''(s)`` — linear in s.
* Hermite-cubic bending about local y (``w(s)``):
with the ``dw/dx = −θ_y`` sign convention used by
``_k_local``, ``M_y(s) = −E·Iy · w''(s)``.
Parameters
----------
coords : (2, 3)
Element node coordinates in global frame.
material, real :
Same as :meth:`ke`. ``real`` carries ``(A, Iz, Iy, J)``
in slots 0..3.
u_e : (12,)
Per-element global-frame DOF vector laid out as
``[ux_I, uy_I, uz_I, θx_I, θy_I, θz_I, ux_J, …]``.
fiber_y, fiber_z :
Fiber location relative to the section centroid in the
beam's local frame. ``y`` is the depth direction (the
orientation node defines local +y), ``z`` is the
cross-bending direction. Defaults to ``z = 0``.
station_s :
Natural station ``s ∈ [0, L]`` (default 0 — node I).
References
----------
* Cook, R.D., Malkus, D.S., Plesha, M.E., Witt, R.J.,
*Concepts and Applications of FEA*, 4th ed., Wiley, 2002,
§2.5 (Hermite-cubic Euler-Bernoulli beam, internal moment
recovery from second derivative of the deflection field).
* Logan, D.L., *A First Course in the Finite Element
Method*, 5th ed., Cengage, 2012, §5.4 (axial-bending
superposition for σ_xx at a section fiber).
"""
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
r = np.asarray(real, dtype=np.float64)
if r.size < 4:
raise ValueError(
f"Beam2.fiber_stress requires REAL[0..3] = AREA, IZZ, IYY, J; got {r.size} values"
)
E = float(material["EX"])
A = float(r[0])
Iz = float(r[1])
Iy = float(r[2])
u = np.asarray(u_e, dtype=np.float64)
if u.shape != (12,):
raise ValueError(f"Beam2.fiber_stress expects u_e shape (12,), got {u.shape}")
# Pull u_e into the element-local frame so the DOF slots align
# with the K_local layout (see _transform_block).
T = _transform_block(R)
u_local = T @ u
# Linear axial: ε_xx = (u_J − u_I)/L → N = E·A·ε_xx
ux_I = float(u_local[0])
ux_J = float(u_local[6])
n_axial = E * A * (ux_J - ux_I) / L
# Hermite cubic bending in x-y plane (DOFs v, θz at slots 1, 5, 7, 11).
v_I = float(u_local[1])
tz_I = float(u_local[5])
v_J = float(u_local[7])
tz_J = float(u_local[11])
s = float(station_s)
if not (0.0 - 1e-12 <= s <= L + 1e-12):
raise ValueError(f"Beam2.fiber_stress: station_s={s} outside [0, L={L}] of element")
# Hermite second derivatives at s:
# H1''(s) = -6/L² + 12s/L³
# H2''(s) = -4/L + 6s/L²
# H3''(s) = +6/L² - 12s/L³
# H4''(s) = -2/L + 6s/L²
h1pp = -6.0 / L**2 + 12.0 * s / L**3
h2pp = -4.0 / L + 6.0 * s / L**2
h3pp = +6.0 / L**2 - 12.0 * s / L**3
h4pp = -2.0 / L + 6.0 * s / L**2
v_dd = h1pp * v_I + h2pp * tz_I + h3pp * v_J + h4pp * tz_J
m_z = E * Iz * v_dd
# Hermite cubic bending in x-z plane (DOFs w, θy at slots 2, 4, 8, 10).
# _k_local uses the dw/dx = -θy sign convention (the -6L cross-terms
# in the bending block). Mirror that here so M_y has the same
# sign convention as the element stiffness produces:
# w(s) = H1·w_I - H2·θy_I + H3·w_J - H4·θy_J
# w''(s) = H1''·w_I - H2''·θy_I + H3''·w_J - H4''·θy_J
# M_y(s) = -E·Iy · w''(s)
w_I = float(u_local[2])
ty_I = float(u_local[4])
w_J = float(u_local[8])
ty_J = float(u_local[10])
w_dd = h1pp * w_I - h2pp * ty_I + h3pp * w_J - h4pp * ty_J
m_y = -E * Iy * w_dd
# σ_xx superposition. Sign convention:
# M_z bends in x-y plane; positive M_z makes the +y fiber
# compress → contribution -y·M_z/Iz.
# M_y bends in x-z plane; positive M_y makes the +z fiber
# compress → contribution -z·M_y/Iy ... but with the local
# sign convention encoded above (dw/dx = -θy), M_y > 0
# actually puts +z in tension, so the contribution flips:
# +z·M_y/Iy.
sigma = n_axial / A - m_z * float(fiber_y) / Iz + m_y * float(fiber_z) / Iy
return float(sigma)
[docs]
@staticmethod
def initial_static_nonlinear_state(
coords: np.ndarray, # noqa: ARG004
material: dict[str, float], # noqa: ARG004
real: np.ndarray, # noqa: ARG004
plastic: ElastoplasticMaterial | None, # noqa: ARG004
fiber_grid: FiberGrid | None = None,
):
"""Initial plastic state — shape depends on plasticity layout.
* No fibre grid: a single :class:`IntegrationPointState` (axial-
only, matches the Phase 3 part 1 path).
* Fibre grid: a list-of-lists indexed ``[axial_station]
[fiber_idx]`` — each fibre at each axial Gauss point carries
its own ``(ε_p, α, ε_p_cum)``.
"""
if fiber_grid is None:
return IntegrationPointState()
n_axial = len(_AXIAL_GAUSS_XI)
n_fiber = fiber_grid.n_fiber
return [[IntegrationPointState() for _ in range(n_fiber)] for _ in range(n_axial)]
[docs]
@staticmethod
def stress_update_static_nonlinear(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
u_e_global: np.ndarray,
state,
plastic: ElastoplasticMaterial | None,
fiber_grid: FiberGrid | None = None,
) -> tuple[np.ndarray, np.ndarray, object]:
"""Phase-3 plasticity hook on a Beam2 — axial or fibre.
Two integration regimes are supported, keyed on ``fiber_grid``:
* **Axial-only** (``fiber_grid is None``): the section's
axial fibre is treated as a single integration point. The
consistent algorithmic ``D_t · A / L`` axial block sits in
the element tangent; bending and torsion remain linear-
elastic. This matches the Phase 3 part 1 path.
* **Fibre integration** (``fiber_grid`` is a :class:`FiberGrid`):
full distributed-plasticity beam-column at 2 axial Gauss
stations × ``n_fiber`` cross-section fibres. Each fibre
runs :func:`stress_update_uniaxial` on its uniaxial strain
``ε_f = ε_axial − y_f · κ_z(s) − z_f · κ_y(s)``. Section
resultants ``N``, ``M_z``, ``M_y`` come from fibre sums;
the consistent algorithmic ``D_section`` couples axial-
bending DOFs through ``Σ_f D_t · y_f · dA`` etc. Torsion
stays linear-elastic.
Parameters
----------
coords : (2, 3)
Element node coordinates in the global frame.
material, real :
Same as :meth:`ke`. ``real`` carries ``(A, Iz, Iy, J)``.
u_e_global : (12,)
Per-element global-frame DOF vector.
state : IntegrationPointState
Trial state at the start of the Newton iteration (committed
state from the previous load step until convergence).
plastic : ElastoplasticMaterial or None
The bilinear (BKIN / BISO) law deposited by ``TB,PLAS,…``.
``None`` means a fully elastic element — return-map is
skipped and the elastic axial stiffness is used.
Returns
-------
K_t_global : (12, 12)
Consistent algorithmic tangent in the global frame.
F_int_global : (12,)
Internal-force vector in the global frame.
new_state : IntegrationPointState
Trial state after this Newton step's return-map (axial).
"""
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
r = np.asarray(real, dtype=np.float64)
if r.size < 4:
raise ValueError(
f"Beam2.stress_update_static_nonlinear requires REAL[0..3] = "
f"AREA, IZZ, IYY, J; got {r.size} values"
)
E = float(material["EX"])
nu = float(material.get("PRXY", 0.3))
G = E / (2.0 * (1.0 + nu))
A = float(r[0])
Iz = float(r[1])
Iy = float(r[2])
J = float(r[3])
T = _transform_block(R)
u_local = T @ np.asarray(u_e_global, dtype=np.float64)
if fiber_grid is not None:
return _stress_update_fibre(
T=T,
u_local=u_local,
L=L,
E=E,
G=G,
A=A,
Iy=Iy,
Iz=Iz,
J=J,
state=state,
plastic=plastic,
fiber_grid=fiber_grid,
)
# Axial-only path (Phase 3 part 1): the section's axial fibre
# is treated as a single integration point. ``ε_xx`` is constant
# along the element under linear axial shape functions.
eps_axial = (u_local[6] - u_local[0]) / L
if plastic is not None:
sigma, new_state, tangent = stress_update_uniaxial(
eps_total=eps_axial,
state=state,
young_modulus=E,
plastic=plastic,
)
else:
sigma = E * eps_axial
tangent = E
new_state = state
# Local K_t — linear-elastic bending/torsion blocks plus the
# plasticity-tangent axial block (replaces E·A/L with D_t·A/L
# so the iteration matrix tracks the algorithmic tangent).
k_loc = _k_local(L, E, G, A, Iy, Iz, J)
k_axial = tangent * A / L
k_loc[0, 0] = k_loc[6, 6] = k_axial
k_loc[0, 6] = k_loc[6, 0] = -k_axial
# Internal-force vector in the local frame. Bending and
# torsion are linear, so ``K · u_local`` is exact for those
# rows. For the axial slots we substitute the explicit ±σ·A —
# the secant stress, which is what the residual ``F_ext − F_int``
# is meant to drive to zero (the tangent ``D_t`` in ``k_loc[0,6]``
# only governs the iteration matrix, not the residual).
f_int_local = k_loc @ u_local
f_int_local[0] = -sigma * A
f_int_local[6] = +sigma * A
K_t_global = T.T @ k_loc @ T
F_int_global = T.T @ f_int_local
return K_t_global, F_int_global, new_state
# ---------------------------------------------------------------------
# Fibre-section integration helpers (Phase 3 part 2)
# ---------------------------------------------------------------------
def _hermite_d2(s: float, L: float) -> tuple[float, float, float, float]:
"""Second derivatives of the four Hermite cubic shape functions
at station ``s ∈ [0, L]`` — linear in ``s``, used for curvature
recovery and B-matrices in bending-DOF integrals.
See :meth:`Beam2.fiber_stress` (which uses the identical four
polynomials)."""
return (
-6.0 / L**2 + 12.0 * s / L**3,
-4.0 / L + 6.0 * s / L**2,
+6.0 / L**2 - 12.0 * s / L**3,
-2.0 / L + 6.0 * s / L**2,
)
def _stress_update_fibre(
*,
T: np.ndarray,
u_local: np.ndarray,
L: float,
E: float,
G: float,
A: float,
Iy: float,
Iz: float,
J: float,
state,
plastic: ElastoplasticMaterial | None,
fiber_grid: FiberGrid,
) -> tuple[np.ndarray, np.ndarray, list]:
"""Distributed-plasticity fibre integration over a Beam2 element.
Two-point Gauss along the axis × ``n_fiber`` cross-section fibres
per station. Returns ``(K_t_global, F_int_global, new_state)``.
Sign / DOF conventions match :func:`_k_local`:
* Local DOFs ``(u, v, w, θx, θy, θz)`` per node.
* x-y bending uses ``(v_I=1, θz_I=5, v_J=7, θz_J=11)``; the
curvature ``κ_z = v''(s)`` is recovered with the standard
Hermite second derivatives.
* x-z bending uses ``(w_I=2, θy_I=4, w_J=8, θy_J=10)`` with
the ``dw/dx = −θ_y`` convention — so the slot signs are
``(N1, −N2, N3, −N4)`` in front of the second-derivative
stencil.
* Fibre strain: ``ε_f = ε_axial − y_f · κ_z(s) − z_f · κ_y(s)``.
* Section resultants: ``N = Σσ·dA``, ``M_z = −Σσ·y·dA``,
``M_y = +Σσ·z·dA``.
"""
n_axial = len(_AXIAL_GAUSS_XI)
if not (isinstance(state, list) and len(state) == n_axial):
raise ValueError(
f"Beam2 fibre stress-update expects state shaped ({n_axial}, n_fiber); "
f"got {type(state).__name__}"
)
fy = fiber_grid.y
fz = fiber_grid.z
fA = fiber_grid.dA
# Axial DOFs.
u_I_axial = float(u_local[0])
u_J_axial = float(u_local[6])
eps_axial = (u_J_axial - u_I_axial) / L
# Bending DOFs (local frame).
v_I, tz_I, v_J, tz_J = (
float(u_local[1]),
float(u_local[5]),
float(u_local[7]),
float(u_local[11]),
)
w_I, ty_I, w_J, ty_J = (
float(u_local[2]),
float(u_local[4]),
float(u_local[8]),
float(u_local[10]),
)
# Container for new fibre states (committed only when the outer
# Newton step converges).
new_state: list[list[IntegrationPointState]] = [
[None] * fiber_grid.n_fiber
for _ in range(n_axial) # type: ignore[list-item]
]
# Local element internal force / tangent.
k_loc = np.zeros((12, 12), dtype=np.float64)
f_int_local = np.zeros(12, dtype=np.float64)
# Torsion stays linear-elastic (Saint-Venant).
kt = G * J / L
k_loc[3, 3] = k_loc[9, 9] = kt
k_loc[3, 9] = k_loc[9, 3] = -kt
f_int_local[3] = kt * (u_local[3] - u_local[9])
f_int_local[9] = -f_int_local[3]
# B-matrix slot maps for scatter into the 12-DOF vector / 12×12 K.
axial_slots = (0, 6)
bend_z_slots = (1, 5, 7, 11)
bend_y_slots = (2, 4, 8, 10)
# Axial B (constant along s).
Ba = np.array([-1.0 / L, +1.0 / L])
for ig, (xi, w_xi) in enumerate(zip(_AXIAL_GAUSS_XI, _AXIAL_GAUSS_WT)):
# Map ξ ∈ [-1, 1] → s ∈ [0, L]; quadrature weight = w_xi · L/2.
s = L * 0.5 * (xi + 1.0)
w_g = w_xi * L * 0.5
h1pp, h2pp, h3pp, h4pp = _hermite_d2(s, L)
# Curvatures at this station.
kappa_z = h1pp * v_I + h2pp * tz_I + h3pp * v_J + h4pp * tz_J
kappa_y = h1pp * w_I - h2pp * ty_I + h3pp * w_J - h4pp * ty_J
# B-vectors at this station.
Bz = np.array([h1pp, h2pp, h3pp, h4pp])
By = np.array([h1pp, -h2pp, h3pp, -h4pp])
# Per-fibre stress update + accumulated section resultants /
# tangent integrals.
N_s = 0.0
Mz_s = 0.0
My_s = 0.0
DA_NN = 0.0
DA_NMz = 0.0
DA_NMy = 0.0
DA_MzMz = 0.0
DA_MzMy = 0.0
DA_MyMy = 0.0
for f, (yf, zf, dA) in enumerate(zip(fy, fz, fA)):
eps_f = eps_axial - yf * kappa_z - zf * kappa_y
if plastic is not None:
sigma, new_state_f, D_t = stress_update_uniaxial(
eps_total=eps_f,
state=state[ig][f],
young_modulus=E,
plastic=plastic,
)
else:
sigma = E * eps_f
D_t = E
new_state_f = state[ig][f]
new_state[ig][f] = new_state_f
N_s += sigma * dA
Mz_s += -sigma * yf * dA
My_s += +sigma * zf * dA
DA_NN += D_t * dA
DA_NMz += -D_t * yf * dA
DA_NMy += +D_t * zf * dA
DA_MzMz += D_t * yf * yf * dA
DA_MzMy += -D_t * yf * zf * dA
DA_MyMy += D_t * zf * zf * dA
# Internal-force scatter (multiply by w_g now).
for k_axial_idx, slot in enumerate(axial_slots):
f_int_local[slot] += Ba[k_axial_idx] * N_s * w_g
for k_bend_idx, slot in enumerate(bend_z_slots):
f_int_local[slot] += Bz[k_bend_idx] * Mz_s * w_g
for k_bend_idx, slot in enumerate(bend_y_slots):
f_int_local[slot] += By[k_bend_idx] * My_s * w_g
# Tangent K scatter — six block-pairs touching the section's
# 3-DOF (axial, M_z, M_y) tangent matrix.
for i_a, slot_a_i in enumerate(axial_slots):
for j_a, slot_a_j in enumerate(axial_slots):
k_loc[slot_a_i, slot_a_j] += Ba[i_a] * DA_NN * Ba[j_a] * w_g
for i_a, slot_a_i in enumerate(axial_slots):
for j_z, slot_z_j in enumerate(bend_z_slots):
contrib = Ba[i_a] * DA_NMz * Bz[j_z] * w_g
k_loc[slot_a_i, slot_z_j] += contrib
k_loc[slot_z_j, slot_a_i] += contrib
for i_a, slot_a_i in enumerate(axial_slots):
for j_y, slot_y_j in enumerate(bend_y_slots):
contrib = Ba[i_a] * DA_NMy * By[j_y] * w_g
k_loc[slot_a_i, slot_y_j] += contrib
k_loc[slot_y_j, slot_a_i] += contrib
for i_z, slot_z_i in enumerate(bend_z_slots):
for j_z, slot_z_j in enumerate(bend_z_slots):
k_loc[slot_z_i, slot_z_j] += Bz[i_z] * DA_MzMz * Bz[j_z] * w_g
for i_z, slot_z_i in enumerate(bend_z_slots):
for j_y, slot_y_j in enumerate(bend_y_slots):
contrib = Bz[i_z] * DA_MzMy * By[j_y] * w_g
k_loc[slot_z_i, slot_y_j] += contrib
k_loc[slot_y_j, slot_z_i] += contrib
for i_y, slot_y_i in enumerate(bend_y_slots):
for j_y, slot_y_j in enumerate(bend_y_slots):
k_loc[slot_y_i, slot_y_j] += By[i_y] * DA_MyMy * By[j_y] * w_g
K_t_global = T.T @ k_loc @ T
F_int_global = T.T @ f_int_local
return K_t_global, F_int_global, new_state