"""Single 8-node hex under uniaxial tension.
A unit cube of isotropic linear elastic material, with one face
clamped axially and the opposite face loaded by a uniform axial
traction, must reproduce 1-D Hooke's law and Poisson contraction
to machine precision:
.. math::
\\sigma_{xx} = E\\, \\varepsilon_{xx}, \\qquad
\\varepsilon_{yy} = \\varepsilon_{zz} = -\\nu\\, \\varepsilon_{xx}.
This is the smallest possible check that the constitutive matrix
:math:`D` is wired correctly — any deviation greater than a few
ULP of the solve tolerance indicates a bug in the element
kernel's material law or the strain recovery.
References
----------
* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)
*Concepts and Applications of Finite Element Analysis*, 4th
ed., Wiley, §1.3 (Hooke's law), §6.2 (single-hex verification).
* Hughes, T. J. R. (2000) *The Finite Element Method*, Dover, §2.7.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any
import numpy as np
import pyvista as pv
import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.validation._benchmark import (
BenchmarkProblem,
PublishedValue,
)
[docs]
@dataclass
class SingleHexUniaxial(BenchmarkProblem):
"""Unit hex under uniaxial traction — Hooke's law + Poisson."""
name: str = "single_hex_uniaxial"
description: str = (
"Single 8-node hex under uniaxial tension — verifies 1-D "
"Hooke's law and Poisson contraction to machine precision."
)
analysis: str = "static"
L: float = 1.0
E: float = 2.0e11
nu: float = 0.3
rho: float = 7850.0
#: Total axial load on the x=L face; distributed across its 4 nodes.
P_total: float = 1.0e6
@property
def published_values(self) -> tuple[PublishedValue, ...]:
# Applied axial stress: total force / face area.
sigma_xx = self.P_total / (self.L * self.L)
eps_xx = sigma_xx / self.E
eps_yy = -self.nu * eps_xx
# Axial tip displacement on the loaded face: epsilon * L.
u_tip = eps_xx * self.L
# Transverse displacement on the loaded face: eps_yy * L (at y = L).
u_trans = eps_yy * self.L
return (
PublishedValue(
name="tip_axial_disp",
value=u_tip,
unit="m",
source="Hughes 2000 §2.7",
formula="u_x(L) = sigma/E * L = P L / (E A)",
# Machine-precision agreement is expected on a single
# hex — use a loose relative tolerance.
tolerance=1e-6,
),
PublishedValue(
name="transverse_disp",
value=u_trans,
unit="m",
source="Hughes 2000 §2.7",
formula="u_y(L) = -nu eps_xx L",
tolerance=1e-6,
),
PublishedValue(
name="sigma_xx_avg",
value=sigma_xx,
unit="Pa",
source="Hughes 2000 §2.7",
formula="sigma_xx = E eps_xx = P / A",
# σ = E ε on a single hex — residual is solver noise.
tolerance=1e-6,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
# Single element — the classic verification geometry.
xs = np.array([0.0, self.L])
ys = np.array([0.0, self.L])
zs = np.array([0.0, self.L])
grid = pv.StructuredGrid(
*np.meshgrid(xs, ys, zs, indexing="ij")
).cast_to_unstructured_grid()
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.HEX8,
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
)
pts = np.asarray(m.grid.points)
# x=0 face: UX=0 on all four nodes (axial clamp).
face_x0 = np.where(pts[:, 0] < 1e-9)[0]
m.fix(nodes=(face_x0 + 1).tolist(), dof="UX")
# Remove rigid-body translation in y and z by pinning one node
# (the x=0, y=0, z=0 corner).
corner = int(np.where(np.linalg.norm(pts, axis=1) < 1e-9)[0][0])
m.fix(nodes=[corner + 1], dof="UY")
m.fix(nodes=[corner + 1], dof="UZ")
# Remove z-rotation by anchoring (0, L, 0).
c0L0 = int(
np.where(
(np.abs(pts[:, 0]) < 1e-9)
& (np.abs(pts[:, 1] - self.L) < 1e-9)
& (np.abs(pts[:, 2]) < 1e-9)
)[0][0]
)
m.fix(nodes=[c0L0 + 1], dof="UZ")
# Uniaxial load on x=L face: distribute total force evenly.
face_xL = np.where(pts[:, 0] > self.L - 1e-9)[0]
f_per_node = self.P_total / len(face_xL)
for n in face_xL:
m.apply_force(int(n + 1), fx=f_per_node)
m._bench_face_xL = face_xL # type: ignore[attr-defined]
m._bench_pts = pts # type: ignore[attr-defined]
return m