"""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])
def _ke_enhanced(coords: np.ndarray, C: np.ndarray, t: float) -> np.ndarray:
"""Q6 / Wilson-Taylor-Beresford enhanced-strain stiffness.
Adds two displacement bubbles ``b₁ = 1 − ξ²`` and ``b₂ = 1 − η²``
(vanishing on element edges) with four internal DOFs ``α``.
Patch-test compatibility on distorted elements requires two
coordinated modifications (Taylor-Beresford-Wilson 1976 / the
Pian-Tong simplification later popularised by Simo-Rifai 1990):
* Bubble derivatives use the *centroid* Jacobian ``J(0, 0)`` so the
enhanced strain matrix ``G`` is linear in ``(ξ, η)``.
* The enhanced-integration determinant is held constant at
``det J(0, 0)`` rather than ``det J(ξ, η)``. Without this,
``K_αu`` does not annihilate the constant-strain modes on a
distorted element and the patch test fails.
The compatible block ``K_uu`` still uses ``det J(ξ, η)`` (otherwise
the standard part loses accuracy on distorted geometries).
The internal ``α`` is statically condensed at element level::
K_eff = K_uu − K_uα · K_αα⁻¹ · K_αu
Returns the condensed 8×8 effective stiffness matrix.
References
----------
* Wilson, E.L., Taylor, R.L., Doherty, W., Ghaboussi, J.,
"Incompatible Displacement Models," in *Numerical and Computer
Methods in Structural Mechanics*, Academic Press, 1973.
* Taylor, R.L., Beresford, P.J., Wilson, E.L., "A Non-conforming
Element for Stress Analysis," *IJNME* 10 (1976), 1211–1219 —
the centroid-Jacobian fix that recovers patch-test compatibility.
* Simo, J.C., Rifai, M.S., "A Class of Mixed Assumed Strain
Methods and the Method of Incompatible Modes," *IJNME* 29
(1990), 1595–1638.
* Cook, Malkus, Plesha, Witt, *Concepts and Applications of FEA*,
4th ed., Wiley 2002, §6.7.
"""
# Centroid Jacobian — drives both the bubble derivatives and the
# constant-determinant integration of the enhanced block.
_, dN_dxi_0 = _shape_and_derivs(0.0, 0.0)
J0 = dN_dxi_0.T @ coords
det_J0 = float(np.linalg.det(J0))
if det_J0 <= 0:
raise ValueError(
f"QUAD4_PLANE: non-positive centroid Jacobian {det_J0:.3e} — check node ordering"
)
J0_inv = np.linalg.inv(J0)
K_uu = np.zeros((8, 8), dtype=np.float64)
K_ua = np.zeros((8, 4), dtype=np.float64)
K_aa = np.zeros((4, 4), dtype=np.float64)
for gp, w in zip(_GAUSS_POINTS, _GAUSS_WEIGHTS):
xi, eta = float(gp[0]), float(gp[1])
_, dN_dxi = _shape_and_derivs(xi, eta)
J = dN_dxi.T @ coords
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
B = _b_matrix(dN_dx)
# Bubble derivatives in natural coords:
# b₁ = 1 − ξ² → ∂b₁/∂ξ = −2ξ, ∂b₁/∂η = 0
# b₂ = 1 − η² → ∂b₂/∂ξ = 0, ∂b₂/∂η = −2η
# Push to physical coords via the *centroid* Jacobian.
db_dxi = np.array(
[[-2.0 * xi, 0.0], [0.0, -2.0 * eta]],
dtype=np.float64,
) # (2 bubbles, 2 nat coords)
db_dx = (J0_inv @ db_dxi.T).T # (2 bubbles, 2 phys coords)
bx1, by1 = db_dx[0]
bx2, by2 = db_dx[1]
# Strain-from-α matrix (3 strain components × 4 enhanced DOFs).
# α layout: [α_{x,1}, α_{x,2}, α_{y,1}, α_{y,2}] (bubble × component).
G = np.array(
[
[bx1, bx2, 0.0, 0.0],
[0.0, 0.0, by1, by2],
[by1, by2, bx1, bx2],
],
dtype=np.float64,
)
BTC = B.T @ C
K_uu += BTC @ B * det_J * w * t
# Enhanced blocks integrated with the *centroid* determinant.
K_ua += BTC @ G * det_J0 * w * t
K_aa += G.T @ C @ G * det_J0 * w * t
return K_uu - K_ua @ np.linalg.solve(K_aa, K_ua.T)
[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)
tech = material.get("_QUAD4_PLANE_TECH", "full")
if tech == "enhanced":
return _ke_enhanced(coords, C, t)
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.get("DENS", 0.0))
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 strain_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.strain_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.strain_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.strain_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