"""Quad4Plane — 2D 4-node structural solid.
Quad element with bilinear shape functions on the reference square
``ξ, η ∈ [-1, 1]``:
N_i(ξ, η) = (1/4)(1 + ξ_i ξ)(1 + η_i η)
Node ordering matches VTK_QUAD — the four corners traversed
counter-clockwise: ``(-,-), (+,-), (+,+), (-,+)``. Two
translational DOFs per node (UX, UY); element lives in the global
(X, Y) plane.
``material["_QUAD4_PLANE_MODE"]`` selects the constitutive regime:
* ``"stress"`` (default) — plane stress (``σ_zz = 0``)
* ``"strain"`` — plane strain (``ε_zz = 0``)
* axisymmetric mode is roadmap and not yet implemented.
Stiffness / mass
----------------
2×2 Gauss integration for both K and M::
K_e = ∫_A Bᵀ C B t dA ≈ Σ_g Bᵀ(ξ_g) C B(ξ_g) det J(ξ_g) t_g
M_e = ρ t ∫_A Nᵀ N dA ≈ Σ_g ρ t Nᵀ(ξ_g) N(ξ_g) det J(ξ_g)
where ``t`` is the out-of-plane thickness (real constant, default 1.0
for plane strain) and ``C`` is the 3×3 plane elastic matrix in
engineering shears.
Real constants
--------------
real[0] : THK — out-of-plane thickness (default 1.0)
Material
--------
EX, PRXY : elastic moduli (both mandatory)
DENS : density (required for M_e)
Lumped mass
-----------
Row-sum (HRZ) lumping distributes each row of the consistent mass onto
its diagonal.
References
----------
* Shape functions — bilinear Lagrange on the 4-node quadrilateral
(Q4 / tensor-product Lagrange on the reference square):
Zienkiewicz, O.C. and Taylor, R.L. *The Finite Element Method:
Its Basis and Fundamentals*, 7th ed., Butterworth-Heinemann,
2013, §6.3.2. Cook, Malkus, Plesha, Witt, *Concepts and
Applications of Finite Element Analysis*, 4th ed., Wiley, 2002,
§6.3.
* 2×2 Gauss-Legendre quadrature (product rule): Zienkiewicz &
Taylor §5.10 Table 5.3.
* Plane-stress / plane-strain elastic matrix: Cook et al. §3.5
(plane stress) and §3.6 (plane strain).
* Row-sum / HRZ lumped mass: Hinton, Rock & Zienkiewicz,
"A note on mass lumping and related processes in the finite
element method," *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
_XI_SIGNS = np.array(
[
[-1, -1],
[+1, -1],
[+1, +1],
[-1, +1],
],
dtype=np.float64,
)
_GP = 1.0 / np.sqrt(3.0)
_GAUSS_POINTS = np.array(
[
[-_GP, -_GP],
[+_GP, -_GP],
[+_GP, +_GP],
[-_GP, +_GP],
],
dtype=np.float64,
)
_GAUSS_WEIGHTS = np.ones(4, dtype=np.float64)
def _shape_and_derivs(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]:
s = _XI_SIGNS
one_p_xi = 1.0 + s[:, 0] * xi
one_p_eta = 1.0 + s[:, 1] * eta
N = 0.25 * one_p_xi * one_p_eta
dN = np.empty((4, 2), dtype=np.float64)
dN[:, 0] = 0.25 * s[:, 0] * one_p_eta
dN[:, 1] = 0.25 * s[:, 1] * one_p_xi
return N, dN
def _b_matrix(dN_dx: np.ndarray) -> np.ndarray:
"""Return the 3×8 strain-displacement matrix from ``dN/dx (4, 2)``.
Engineering-shear ordering: ``[ε_xx, ε_yy, γ_xy]``.
"""
B = np.zeros((3, 8), dtype=np.float64)
for i in range(4):
bx, by = dN_dx[i]
col = 2 * i
B[0, col + 0] = bx
B[1, col + 1] = by
B[2, col + 0] = by
B[2, col + 1] = bx
return B
def _shape_matrix(N: np.ndarray) -> np.ndarray:
Phi = np.zeros((2, 8), dtype=np.float64)
for i in range(4):
Phi[0, 2 * i + 0] = N[i]
Phi[1, 2 * i + 1] = N[i]
return Phi
def _plane_stress_C(E: float, nu: float) -> np.ndarray:
c = E / (1.0 - nu * nu)
C = np.array(
[
[c, c * nu, 0.0],
[c * nu, c, 0.0],
[0.0, 0.0, c * (1.0 - nu) / 2.0],
],
dtype=np.float64,
)
return C
def _plane_strain_C(E: float, nu: float) -> np.ndarray:
c = E / ((1.0 + nu) * (1.0 - 2.0 * nu))
C = np.array(
[
[c * (1.0 - nu), c * nu, 0.0],
[c * nu, c * (1.0 - nu), 0.0],
[0.0, 0.0, c * (1.0 - 2.0 * nu) / 2.0],
],
dtype=np.float64,
)
return C
def _thickness(real: np.ndarray) -> float:
real = np.asarray(real, dtype=np.float64)
if real.size == 0:
return 1.0
return float(real[0])
[docs]
@register
class Quad4Plane(ElementBase):
name = "QUAD4_PLANE"
n_nodes = 4
n_dof_per_node = 2 # UX, UY
vtk_cell_type = 9 # VTK_QUAD
# Plane-mode default; overridable via ``material["_QUAD4_PLANE_MODE"]``.
plane_mode: str = "stress"
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
coords = np.asarray(coords, dtype=np.float64)
if coords.shape == (4, 3):
coords = coords[:, :2]
if coords.shape != (4, 2):
raise ValueError(f"Quad4Plane expects (4, 2) or (4, 3) coords, got {coords.shape}")
E = float(material["EX"])
nu = float(material["PRXY"])
mode = material.get("_QUAD4_PLANE_MODE", "stress")
C = _plane_stress_C(E, nu) if mode == "stress" else _plane_strain_C(E, nu)
t = _thickness(real)
K = np.zeros((8, 8), dtype=np.float64)
for gp, w in zip(_GAUSS_POINTS, _GAUSS_WEIGHTS):
_, dN_dxi = _shape_and_derivs(*gp)
J = dN_dxi.T @ coords # 2×2
det_J = float(np.linalg.det(J))
if det_J <= 0:
raise ValueError(
f"QUAD4_PLANE: non-positive Jacobian {det_J:.3e} — check node ordering"
)
dN_dx = np.linalg.solve(J, dN_dxi.T).T # (4, 2)
B = _b_matrix(dN_dx)
K += B.T @ C @ B * det_J * w * t
return K
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
coords = np.asarray(coords, dtype=np.float64)
if coords.shape == (4, 3):
coords = coords[:, :2]
if coords.shape != (4, 2):
raise ValueError(f"Quad4Plane expects (4, 2) or (4, 3) coords, got {coords.shape}")
rho = float(material["DENS"])
t = _thickness(real)
M = np.zeros((8, 8), dtype=np.float64)
for gp, w in zip(_GAUSS_POINTS, _GAUSS_WEIGHTS):
N, dN_dxi = _shape_and_derivs(*gp)
J = dN_dxi.T @ coords
det_J = float(np.linalg.det(J))
Phi = _shape_matrix(N)
M += rho * t * (Phi.T @ Phi) * det_J * w
if lumped:
M = np.diag(M.sum(axis=1))
return M
[docs]
@staticmethod
def eel_batch(
coords: np.ndarray,
u_e: np.ndarray,
material: dict[str, float] | None = None,
) -> np.ndarray | None:
"""Per-element engineering-shear Voigt strain at element nodes.
``coords``: ``(n_elem, 4, 3)`` or ``(n_elem, 4, 2)``.
``u_e``: ``(n_elem, 8)`` — DOF layout ``[ux0, uy0, ux1, …]``.
Returns ``(n_elem, 4, 6)`` Voigt strain ordered
``[εxx, εyy, εzz, γxy, γyz, γxz]`` — matches the 3D-solid Voigt
convention so :func:`compute_nodal_stress` can drive the same
``σ = C·ε`` contraction with the standard isotropic ``(6, 6)``
``C``.
For plane stress (default) we fill ``εzz = -ν/(1-ν)·(εxx+εyy)``
so ``σ_zz`` recovered via the 6×6 ``C`` is identically zero (the
defining condition of plane stress). For plane strain we fill
``εzz = 0`` (the defining condition of plane strain); ``σ_zz``
is then ``ν·(σ_xx + σ_yy)`` per Hooke's law.
``material`` must be provided when the kernel is operating in
plane stress (the default) so ``ν`` is available; for plane
strain the ``material`` argument is unused.
Strain is evaluated at each corner node (ξ_i, η_i) via the
bilinear B-matrix — the same basis used in ``ke``. This is
the standard Q4 nodal recovery; the recovered strain is exact
at the centroid and accurate to ``O(h²)`` at corner nodes for
smooth solutions.
"""
coords = np.asarray(coords, dtype=np.float64)
if coords.ndim != 3 or coords.shape[1] != 4:
raise ValueError(
f"QUAD4_PLANE.eel_batch expects (n_elem, 4, 2|3) coords, got {coords.shape}"
)
if coords.shape[2] == 3:
coords = coords[:, :, :2]
u_e = np.asarray(u_e, dtype=np.float64)
if u_e.ndim != 2 or u_e.shape[1] != 8 or u_e.shape[0] != coords.shape[0]:
raise ValueError(
f"QUAD4_PLANE.eel_batch expects (n_elem, 8) u_e matching "
f"coords[0], got coords {coords.shape}, u_e {u_e.shape}"
)
n_elem = coords.shape[0]
eps = np.zeros((n_elem, 4, 6), dtype=np.float64)
for i, (xi, eta) in enumerate(_XI_SIGNS):
_, dN_dxi = _shape_and_derivs(float(xi), float(eta))
# Per-element Jacobian: J[k] = dN_dxi.T @ coords[k] -> (n_elem, 2, 2)
J = np.einsum("ij,kjm->kim", dN_dxi.T, coords)
J_inv = np.linalg.inv(J) # (n_elem, 2, 2)
# dN_dx[k, n, d] = sum_j J_inv[k, d, j] * dN_dxi[n, j]
# i.e. shape (n_elem, 4 nodes, 2 dims) ready for B-matrix
# assembly.
dN_dx = np.einsum("kij,nj->kni", J_inv, dN_dxi)
# Hand-build B for the batch: shape (n_elem, 3, 8).
B = np.zeros((n_elem, 3, 8), dtype=np.float64)
for nd in range(4):
bx = dN_dx[:, nd, 0]
by = dN_dx[:, nd, 1]
B[:, 0, 2 * nd + 0] = bx
B[:, 1, 2 * nd + 1] = by
B[:, 2, 2 * nd + 0] = by
B[:, 2, 2 * nd + 1] = bx
# eps_plane[k] = B[k] @ u_e[k] -> (n_elem, 3)
eps_plane = np.einsum("kij,kj->ki", B, u_e)
# Voigt slot: [εxx, εyy, εzz, γxy, γyz, γxz].
eps[:, i, 0] = eps_plane[:, 0] # εxx
eps[:, i, 1] = eps_plane[:, 1] # εyy
eps[:, i, 3] = eps_plane[:, 2] # γxy
# εzz fill-in depending on plane mode.
mode = material.get("_QUAD4_PLANE_MODE", "stress") if material is not None else "stress"
if mode == "stress":
if material is None:
raise ValueError(
"QUAD4_PLANE.eel_batch needs material for plane-stress "
"εzz fill-in; pass {'EX', 'PRXY'}"
)
nu = float(material["PRXY"])
eps[:, i, 2] = -nu / (1.0 - nu) * (eps_plane[:, 0] + eps_plane[:, 1])
# else plane strain: εzz = 0 (already initialized).
return eps