"""Truss2 — 3D 2-node spar (truss).
Linear shape functions in the natural coordinate ``s ∈ [-1, 1]`` for
all three translations. Only axial strain contributes to internal
work::
B = (1/L) [-l, -m, -n, l, m, n] (1 × 6)
K_e = (E·A/L) · [[ C, -C],
[-C, C]] where C_ij = d_i·d_j (3 × 3)
with ``(l, m, n)`` the unit vector along the element axis.
Consistent mass::
M_e = ρ·A·L / 6 · [[2 I₃, I₃],
[ I₃, 2 I₃]]
Lumped mass: ``M_e = ρ·A·L/2 · I₆``.
Real constants
--------------
real[0]: AREA — cross-sectional area (mandatory)
real[1]: ADDMAS — added mass per unit length (optional; unused here)
real[2]: ISTRN — initial strain (optional; unused here — stress path only)
Material properties
-------------------
EX : Young's modulus (mandatory)
DENS : mass density (required for M_e)
References
----------
* Shape functions — linear Lagrange on a 1-D 2-node element
(``N₁ = ½(1 − s)``, ``N₂ = ½(1 + s)`` on ``s ∈ [−1, 1]``):
Zienkiewicz, O.C. and Taylor, R.L. *The Finite Element Method:
Its Basis and Fundamentals*, 7th ed., Butterworth-Heinemann,
2013, §2.3. Cook, Malkus, Plesha, Witt, *Concepts and
Applications of Finite Element Analysis*, 4th ed., Wiley, 2002,
§2.4.
* Axial-only strain (spar / truss element, no bending / shear
contribution, only longitudinal stiffness): Cook et al. §2.3.
* Consistent mass ``ρAL/6 · [[2, 1], [1, 2]]`` on the two axial
DOFs (and its 6×6 extension to 3-D for the three translational
pairs): Cook et al. Table 16.3-1; Zienkiewicz & Taylor §16.2.1.
* Row-sum lumped mass ``ρAL/2`` per translational DOF: Hinton,
Rock & Zienkiewicz, *Earthquake Engng. Struct. Dyn.* 4 (1976),
pp. 245–249.
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
def _axial_frame(coords: np.ndarray) -> tuple[float, np.ndarray]:
"""Return element length and unit direction vector (I → J)."""
d = coords[1] - coords[0]
L = float(np.linalg.norm(d))
if L == 0.0:
raise ValueError("TRUSS2: coincident nodes, element length is zero")
return L, d / L
[docs]
@register
class Truss2(ElementBase):
name = "TRUSS2"
n_nodes = 2
n_dof_per_node = 3 # UX, UY, UZ
vtk_cell_type = 3 # VTK_LINE
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
L, d = _axial_frame(np.asarray(coords, dtype=np.float64))
E = float(material["EX"])
A = float(np.asarray(real)[0])
C = np.outer(d, d) # 3 × 3
k = (E * A / L) * np.block([[C, -C], [-C, C]])
return np.ascontiguousarray(k)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
L, _ = _axial_frame(np.asarray(coords, dtype=np.float64))
rho = float(material.get("DENS", 0.0))
A = float(np.asarray(real)[0])
m_total = rho * A * L
if lumped:
return (m_total / 2.0) * np.eye(6)
I3 = np.eye(3)
return (m_total / 6.0) * np.block([[2.0 * I3, I3], [I3, 2.0 * I3]])
[docs]
@staticmethod
def stress_update_static_nonlinear(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
u_e_global: np.ndarray,
state,
plastic,
):
"""Phase-2 plasticity hook — uniaxial return-map on a TRUSS2.
Returns
-------
K_t_global : (6, 6)
Consistent algorithmic tangent in the global frame
(block ``e_hat ⊗ e_hat`` scaled by ``D_t·A/L``).
F_int_global : (6,)
Internal-force vector ``±σ·A·e_hat`` at each node.
new_state :
Trial state after this Newton step's return-map.
"""
from femorph_solver.materials.plasticity import stress_update_uniaxial
L, e_hat = _axial_frame(np.asarray(coords, dtype=np.float64))
E = float(material["EX"])
A = float(np.asarray(real, dtype=np.float64)[0])
u = np.asarray(u_e_global, dtype=np.float64)
u_i = u[:3]
u_j = u[3:]
eps_total = float(np.dot(u_j - u_i, e_hat) / L)
if plastic is not None:
sigma, new_state, tangent = stress_update_uniaxial(
eps_total=eps_total,
state=state,
young_modulus=E,
plastic=plastic,
)
else:
sigma = E * eps_total
tangent = E
new_state = state
f_axial = sigma * A
F_int_global = np.concatenate([-f_axial * e_hat, +f_axial * e_hat])
k_axial = tangent * A / L
e_block = np.outer(e_hat, e_hat)
K_t_global = k_axial * np.block([[e_block, -e_block], [-e_block, e_block]])
return np.ascontiguousarray(K_t_global), F_int_global, new_state
[docs]
@staticmethod
def strain_batch(
coords: np.ndarray,
u_e: np.ndarray,
material: dict[str, float] | None = None,
) -> np.ndarray | None:
"""Per-element 3D Voigt strain at the two end nodes.
The truss carries only an axial strain ``eps_a = (u2 - u1) . e / L``
- the perpendicular cross-section stretches under Poisson
contraction when the bar is in uniaxial stress. Returning the
full 3D continuum strain that corresponds to the pure-uniaxial
stress state ``sigma = E*eps_a * (e^e)`` lets
:func:`compute_nodal_stress` recover the right stress via the
standard isotropic ``C @ eps`` contraction without having to
special-case 1D kernels.
For a uniaxial stress state along unit axis ``e`` with strain
``eps_a`` (computed from the displacement field), the 3D strain
tensor is::
eps = eps_a * ((1 + nu) * (e^e) - nu * I)
Voigt order is ``[exx, eyy, ezz, gxy, gyz, gxz]`` (engineering
shears), matching :meth:`MaterialTable.eval_c_isotropic`.
Returns the per-element-node strain (same value at both nodes
since the truss is constant-strain).
"""
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
if u_e.ndim != 2 or u_e.shape[1] != 6 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)
bad = L <= 0.0
if bad.any():
return None
e_hat = d / L[:, None]
# Axial strain at each element: (u2 - u1) . e_hat / L.
du = u_e[:, 3:6] - u_e[:, 0:3]
eps_a = np.einsum("ij,ij->i", du, e_hat) / L
# Build per-element Voigt strain
# eps_v = eps_a * ((1+nu) (e^e) - nu I) (Voigt order [xx,yy,zz,xy,yz,xz])
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)
# Engineering shears: g_ij = 2 * eps_ij; eps_ij = eps_a * (1+nu) * ei*ej.
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
# Replicate the constant-strain tensor at both nodes.
return np.repeat(eps_v[:, None, :], 2, axis=1)
[docs]
@staticmethod
def thermal_load(
coords: np.ndarray,
real: np.ndarray,
material: dict[str, float],
*,
delta_t: float,
) -> np.ndarray:
"""Equivalent nodal force vector for a uniform thermal pre-strain.
For an axial bar with thermal-expansion coefficient ``α =
material["ALPX"]`` and uniform temperature change ``ΔT``, the
thermal pre-strain is ``ε_th = α·ΔT``. The work-equivalent
nodal force vector follows from ``f_th = ∫ Bᵀ E ε_th dV``::
f_local_axial = E · A · α · ΔT · [-1, +1] (along axis)
i.e. an outward force ``E·A·α·ΔT`` at each end pulling the
constrained bar apart on heating (or pushing it together on
cooling). The return value is the 6-component **global-frame**
force vector at ``[ux_I, uy_I, uz_I, ux_J, uy_J, uz_J]``,
already rotated by the same direction-cosine matrix
:func:`ke` uses, so the caller can scatter it directly into ``F``.
Returns ``np.zeros(6)`` when the material has no ``ALPX`` set —
a missing thermal-expansion coefficient just means thermally
inert, not an error.
References
----------
* Cook, R.D., Malkus, D.S., Plesha, M.E., Witt, R.J., *Concepts
and Applications of Finite Element Analysis*, 4th ed., Wiley,
2002, §3.6.4 (thermal pre-strain), §6.10 (initial-strain
equivalent nodal forces).
"""
alpha = float(material.get("ALPX", 0.0))
if alpha == 0.0 or float(delta_t) == 0.0:
return np.zeros(6, dtype=np.float64)
L, d = _axial_frame(np.asarray(coords, dtype=np.float64))
del L # length cancels in the f_th formula (B has 1/L; ∫ … dV gives back L)
E = float(material["EX"])
A = float(np.asarray(real)[0])
f_axial = E * A * alpha * float(delta_t)
# Local-frame [u_I_axial, u_J_axial] = f_axial · [-1, +1]; rotate
# to global via the unit direction vector. In the 6-DOF global
# frame, node-I axial component sits in (UX_I, UY_I, UZ_I) along
# ``d`` (and node-J in the opposite sense).
f_global = np.empty(6, dtype=np.float64)
f_global[0:3] = -f_axial * d
f_global[3:6] = +f_axial * d
return f_global