"""Quad4Shell — 4-node flat Mindlin–Reissner shell.
Kinematics (first-order shear-deformation theory):
u_point = u + z · θ_y
v_point = v − z · θ_x
w_point = w
giving a split into membrane / bending / transverse-shear energy after
through-thickness integration. Node ordering matches VTK_QUAD:
(-,-), (+,-), (+,+), (-,+) in natural ``(ξ, η) ∈ [-1, 1]``.
The element is built in a local 2-D frame (``ex`` along edge 1→2,
``ez`` normal to the best-fit plane, ``ey = ez × ex``). Coordinates are
projected into this frame, the plane stiffness is assembled, then a
block-diagonal 3×3 direction-cosine matrix rotates every node's
six-DOF block back to global coordinates.
Integration
-----------
* Membrane + bending: 2×2 Gauss (full)
* Transverse shear: 1×1 Gauss (selective reduced) — suppresses shear
locking for thin shells without needing MITC4 tying.
Drilling DOF (θ_z, local) has no true stiffness for a flat shell;
femorph-solver adds a small stabilization ``α · G · t · A_elem`` per corner so
the 24×24 local matrix stays non-singular when every DOF is free.
Real constants
--------------
real[0]: THICKNESS — shell thickness (mandatory)
Material
--------
EX, PRXY : plane-stress elastic moduli (mandatory)
DENS : mass density (required for M_e)
References
----------
* **Mindlin--Reissner first-order shear-deformation theory**
(transverse shear included so rotations are independent of
in-plane slope): Mindlin, R.D., "Influence of rotatory inertia
and shear on flexural motions of isotropic elastic plates,"
*J. Appl. Mech.* 18 (1951), pp. 31–38. Reissner, E., "The
effect of transverse shear deformation on the bending of
elastic plates," *J. Appl. Mech.* 12 (1945), pp. A69–A77.
* Bilinear isoparametric 4-node shell element with 6 DOF / node:
Zienkiewicz, O.C. and Taylor, R.L., *The Finite Element Method*,
Vol. 2: *Solid Mechanics*, 5th ed., Butterworth-Heinemann,
2000, Ch. 8 (plate & shell elements).
* **Selective reduced integration for transverse shear** (1×1 on
shear, 2×2 on membrane + bending — suppresses shear locking
for thin shells): Malkus, D.S. and Hughes, T.J.R., "Mixed
finite element methods — reduced and selective integration
techniques: a unification of concepts," *Comput. Methods Appl.
Mech. Engrg.* 15 (1978), pp. 63–81.
* Alternative MITC4 (Mixed-Interpolated Tensorial Components)
shear-locking treatment — roadmap for a future variant:
Bathe, K.J. and Dvorkin, E.N., "A four-node plate bending
element based on Mindlin/Reissner plate theory and a mixed
interpolation," *IJNME* 21 (1985), pp. 367–383.
* Drilling-DOF stabilisation (``α·G·t·A`` per corner to keep the
local 24×24 non-singular when θ_z is free): Allman, D.J.,
"A compatible triangular element including vertex rotations
for plane elasticity analysis," *Comput. Struct.* 19 (1984),
pp. 1–8, and Hughes, T.J.R. and Brezzi, F., "On drilling
degrees of freedom," *Comput. Methods Appl. Mech. Engrg.* 72
(1989), pp. 105–121.
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
_GP = 1.0 / np.sqrt(3.0)
_GAUSS_2x2 = np.array(
[
[-_GP, -_GP],
[+_GP, -_GP],
[+_GP, +_GP],
[-_GP, +_GP],
],
dtype=np.float64,
)
_W_2x2 = np.ones(4, dtype=np.float64)
_GAUSS_1x1 = np.array([[0.0, 0.0]], dtype=np.float64)
_W_1x1 = np.array([4.0], dtype=np.float64)
_SHEAR_FACTOR = 5.0 / 6.0 # Reissner–Mindlin shear correction
_DRILL_ALPHA = 1.0e-3
def _local_frame(coords: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(coords_2d (4,2), R (3,3))`` for a 4-node quad in 3D."""
v1 = coords[1] - coords[0]
v2 = coords[3] - coords[0]
ex = v1 / np.linalg.norm(v1)
normal = np.cross(v1, v2)
n_mag = np.linalg.norm(normal)
if n_mag == 0.0:
raise ValueError("QUAD4_SHELL: degenerate element (zero area)")
ez = normal / n_mag
ey = np.cross(ez, ex)
R = np.stack([ex, ey, ez])
relative = coords - coords[0]
coords_2d = (relative @ R.T)[:, :2]
return coords_2d, R
def _shape_and_derivs(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]:
s = np.array([[-1, -1], [+1, -1], [+1, +1], [-1, +1]], dtype=np.float64)
N = 0.25 * (1.0 + s[:, 0] * xi) * (1.0 + s[:, 1] * eta)
dN = np.empty((4, 2), dtype=np.float64)
dN[:, 0] = 0.25 * s[:, 0] * (1.0 + s[:, 1] * eta)
dN[:, 1] = 0.25 * s[:, 1] * (1.0 + s[:, 0] * xi)
return N, dN
def _d_plane_stress(E: float, nu: float) -> np.ndarray:
factor = E / (1.0 - nu * nu)
return factor * np.array(
[
[1.0, nu, 0.0],
[nu, 1.0, 0.0],
[0.0, 0.0, 0.5 * (1.0 - nu)],
]
)
def _assemble_b(dN_dx: np.ndarray, N: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Return (B_m (3,24), B_b (3,24), B_s (2,24)) in local 6-DOF-per-node ordering."""
Bm = np.zeros((3, 24), dtype=np.float64)
Bb = np.zeros((3, 24), dtype=np.float64)
Bs = np.zeros((2, 24), dtype=np.float64)
for i in range(4):
bx, by = dN_dx[i]
base = 6 * i
# Membrane (u, v)
Bm[0, base + 0] = bx
Bm[1, base + 1] = by
Bm[2, base + 0] = by
Bm[2, base + 1] = bx
# Bending (θx, θy): κxx = ∂θy/∂x, κyy = -∂θx/∂y, κxy = ∂θy/∂y - ∂θx/∂x
Bb[0, base + 4] = bx
Bb[1, base + 3] = -by
Bb[2, base + 4] = by
Bb[2, base + 3] = -bx
# Transverse shear (w, θx, θy): γxz = ∂w/∂x + θy, γyz = ∂w/∂y - θx
Bs[0, base + 2] = bx
Bs[0, base + 4] = N[i]
Bs[1, base + 2] = by
Bs[1, base + 3] = -N[i]
return Bm, Bb, Bs
def _k_local(coords_2d: np.ndarray, E: float, nu: float, t: float) -> np.ndarray:
D_p = _d_plane_stress(E, nu)
D_m = t * D_p
D_b = (t**3 / 12.0) * D_p
G = E / (2.0 * (1.0 + nu))
D_s = _SHEAR_FACTOR * G * t * np.eye(2)
K = np.zeros((24, 24), dtype=np.float64)
area = 0.0
# Membrane + bending: full 2x2 Gauss
for (xi, eta), w in zip(_GAUSS_2x2, _W_2x2):
N, dN_dxi = _shape_and_derivs(xi, eta)
J = dN_dxi.T @ coords_2d
det_J = float(np.linalg.det(J))
dN_dx = np.linalg.solve(J, dN_dxi.T).T
Bm, Bb, _ = _assemble_b(dN_dx, N)
K += (Bm.T @ D_m @ Bm + Bb.T @ D_b @ Bb) * det_J * w
area += det_J * w
# Transverse shear: reduced 1x1 Gauss
for (xi, eta), w in zip(_GAUSS_1x1, _W_1x1):
N, dN_dxi = _shape_and_derivs(xi, eta)
J = dN_dxi.T @ coords_2d
det_J = float(np.linalg.det(J))
dN_dx = np.linalg.solve(J, dN_dxi.T).T
_, _, Bs = _assemble_b(dN_dx, N)
K += Bs.T @ D_s @ Bs * det_J * w
# Drilling stabilization
k_drill = _DRILL_ALPHA * G * t * area
for i in range(4):
K[6 * i + 5, 6 * i + 5] += k_drill
return K
def _m_local(coords_2d: np.ndarray, rho: float, t: float) -> np.ndarray:
M = np.zeros((24, 24), dtype=np.float64)
trans_factor = rho * t
rot_factor = rho * t**3 / 12.0
for (xi, eta), w in zip(_GAUSS_2x2, _W_2x2):
N, dN_dxi = _shape_and_derivs(xi, eta)
J = dN_dxi.T @ coords_2d
det_J = float(np.linalg.det(J))
NN = np.outer(N, N) * det_J * w
for i in range(4):
for j in range(4):
# Translations (u, v, w)
for d in (0, 1, 2):
M[6 * i + d, 6 * j + d] += trans_factor * NN[i, j]
# Rotations (θx, θy) — rotary inertia
for d in (3, 4):
M[6 * i + d, 6 * j + d] += rot_factor * NN[i, j]
return M
def _transform_block(R: np.ndarray) -> np.ndarray:
T = np.zeros((24, 24), dtype=np.float64)
for i in range(8):
T[3 * i : 3 * i + 3, 3 * i : 3 * i + 3] = R
return T
[docs]
@register
class Quad4Shell(ElementBase):
name = "QUAD4_SHELL"
n_nodes = 4
n_dof_per_node = 6
dof_labels = ("UX", "UY", "UZ", "ROTX", "ROTY", "ROTZ")
vtk_cell_type = 9 # VTK_QUAD
[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):
raise ValueError(f"Quad4Shell expects (4, 3) coords, got {coords.shape}")
r = np.asarray(real, dtype=np.float64)
if r.size < 1:
raise ValueError("Quad4Shell requires REAL[0]=thickness; got empty real set")
t = float(r[0])
E = float(material["EX"])
nu = float(material["PRXY"])
coords_2d, R = _local_frame(coords)
K_loc = _k_local(coords_2d, E, nu, t)
T = _transform_block(R)
return T.T @ K_loc @ T
[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):
raise ValueError(f"Quad4Shell expects (4, 3) coords, got {coords.shape}")
r = np.asarray(real, dtype=np.float64)
if r.size < 1:
raise ValueError("Quad4Shell requires REAL[0]=thickness; got empty real set")
t = float(r[0])
rho = float(material["DENS"])
coords_2d, R = _local_frame(coords)
M_loc = _m_local(coords_2d, rho, t)
if lumped:
M_loc = np.diag(M_loc.sum(axis=1))
T = _transform_block(R)
return T.T @ M_loc @ T