"""PLANE182 — 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 MAPDL (and VTK_QUAD) — the four corners
traversed counter-clockwise: ``(-,-), (+,-), (+,+), (-,+)``. Two
translational DOFs per node (UX, UY); element lives in the global
(X, Y) plane.
KEYOPT(3) selects the constitutive regime:
* 0 — plane stress (default, ``σ_zz = 0``)
* 1 — axisymmetric (not yet implemented)
* 2 — plane strain (``ε_zz = 0``)
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.
MAPDL compatibility — specification source
------------------------------------------
* Ansys, Inc., *Ansys Mechanical APDL Element Reference*,
Release 2022R2, section "PLANE182 — 2-D 4-Node Structural
Solid".
Short factual summary (paraphrased): 4-node 2-D solid; 2
translational DOFs per node; ``KEYOPT(3)`` selects plane stress
/ axisymmetric / plane strain; ``REAL[0]`` supplies out-of-plane
thickness in plane-stress mode. Ansys Element Reference is the
compat spec only; theory is Zienkiewicz-Taylor / Cook above.
"""
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 PLANE182(ElementBase):
name = "PLANE182"
n_nodes = 4
n_dof_per_node = 2 # UX, UY
vtk_cell_type = 9 # VTK_QUAD
# KEYOPT(3): 0 = plane stress (default), 2 = plane strain.
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"PLANE182 expects (4, 2) or (4, 3) coords, got {coords.shape}")
E = float(material["EX"])
nu = float(material["PRXY"])
mode = material.get("_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"PLANE182: 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"PLANE182 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