"""Irons constant-strain patch test.
A patch of distorted elements subjected to boundary displacements
consistent with a **uniform strain field** must reproduce that
strain field exactly, to machine precision, in every element. A
conforming FE discretisation that fails this test cannot pass the
inf-sup / consistency conditions and therefore cannot converge at
the optimal rate.
Implementation: apply linear boundary displacements
:math:`u = \\varepsilon_0 \\cdot x` on every boundary node of a
2×2×2 hex patch with an **interior node displaced** from its
regular position (to stress the isoparametric mapping), and check
that every Gauss point in the interior recovers
:math:`\\varepsilon = \\varepsilon_0` to machine precision.
References
----------
* Irons, B. M. and Razzaque, A. (1972) "Experience with the
patch test for convergence of finite elements," in *The
Mathematical Foundations of the Finite Element Method*, pp.
557-587.
* Zienkiewicz, O. C., Taylor, R. L. and Zhu, J. Z. (2013)
*The Finite Element Method: Its Basis and Fundamentals*, 7th
ed., Butterworth-Heinemann, §10.3.
"""
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 IronsPatchTest(BenchmarkProblem):
"""2×2×2 hex patch under imposed uniform strain (Irons 1972)."""
name: str = "patch_test"
description: str = (
"Irons constant-strain patch test: prescribed linear "
"boundary displacement must reproduce a uniform strain to "
"machine precision."
)
analysis: str = "static"
L: float = 1.0
E: float = 2.0e11
nu: float = 0.3
rho: float = 7850.0
#: Target uniaxial strain imposed via boundary displacements.
#: The published expectation is that every interior DOF
#: reproduces ``u = eps0 * x`` at machine precision.
eps0: float = 1.0e-4
@property
def published_values(self) -> tuple[PublishedValue, ...]:
# Max deviation of computed strain from eps0 — expected zero
# to machine precision; set the tolerance at 1e-9 (the
# conservative FE-solve residual floor).
return (
PublishedValue(
name="max_strain_error",
value=0.0,
unit="",
source="Irons & Razzaque 1972",
formula="max_e |eps_xx^h(x_e) - eps0| == 0 (machine precision)",
tolerance=1.0e-9,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
# 2×2×2 hex patch with an interior node displaced so the
# isoparametric mapping is non-trivial. The displacement
# below preserves a convex element — the 8 corner nodes
# stay regular; the single interior node (node 14 in a
# 3×3×3 lattice) is offset by ~10% of the edge length.
xs = np.linspace(0.0, self.L, 3)
ys = np.linspace(0.0, self.L, 3)
zs = np.linspace(0.0, self.L, 3)
grid = pv.StructuredGrid(
*np.meshgrid(xs, ys, zs, indexing="ij")
).cast_to_unstructured_grid()
# Jitter the single interior node (x=y=z=L/2) by 10% of
# the edge length in a fixed direction.
pts = np.asarray(grid.points)
interior = np.argmin(
np.linalg.norm(pts - np.array([self.L / 2, self.L / 2, self.L / 2]), axis=1)
)
pts = pts.copy()
pts[interior] = pts[interior] + 0.1 * self.L * np.array([1.0, 0.5, 0.3])
grid.points = pts
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.HEX8,
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
)
# Prescribe linear boundary displacement ``u_x = eps0 * x``
# on every boundary node; leave ``u_y, u_z`` free (the
# Poisson-contraction-adjusted response is what the interior
# should recover). Simpler: fix UX on the full boundary to
# eps0 * x value, and pin UY / UZ at corner nodes only to
# remove rigid-body motion.
tol = 1e-9
on_x_edge = (np.abs(pts[:, 0]) < tol) | (np.abs(pts[:, 0] - self.L) < tol)
on_y_edge = (np.abs(pts[:, 1]) < tol) | (np.abs(pts[:, 1] - self.L) < tol)
on_z_edge = (np.abs(pts[:, 2]) < tol) | (np.abs(pts[:, 2] - self.L) < tol)
boundary = on_x_edge | on_y_edge | on_z_edge
# Pin UY at the z=0, y=0 edge; UZ at the y=0, z=0 edge; UX
# on each boundary node to ``eps0 * x``. We don't have a
# generic prescribed-displacement API in the native surface
# yet — apply_force with a stiffness trick is unnecessary if
# we simply fix UX on the two x-faces to the known value.
for n_idx in np.where(boundary)[0]:
# Fix UX at x=0 to zero, at x=L to eps0 * L, via the
# deprecated ``Model.d`` APDL-verb path since ``fix`` pins
# to zero only. The deprecation is surfaced as a warning;
# this path is acceptable since the framework itself is
# documentation + regression.
x_val = float(pts[n_idx, 0])
ux_prescribed = self.eps0 * x_val
if abs(pts[n_idx, 0]) < tol:
m.fix(nodes=[int(n_idx + 1)], dof="UX")
elif abs(pts[n_idx, 0] - self.L) < tol:
m.fix(nodes=[int(n_idx + 1)], dof="UX", value=ux_prescribed)
# Pin rigid-body y/z at one corner.
c000 = int(np.where(np.linalg.norm(pts, axis=1) < tol)[0][0])
m.fix(nodes=[c000 + 1], dof="UY")
m.fix(nodes=[c000 + 1], dof="UZ")
c0L0 = int(
np.where(
(np.abs(pts[:, 0]) < tol)
& (np.abs(pts[:, 1] - self.L) < tol)
& (np.abs(pts[:, 2]) < tol)
)[0][0]
)
m.fix(nodes=[c0L0 + 1], dof="UZ")
m._bench_pts = pts # type: ignore[attr-defined]
return m