r"""
Irons constant-strain patch test
================================

The canonical FE consistency test: a 2×2×2 hex patch with an
interior node displaced from its regular position, subjected to
boundary displacements :math:`u = \varepsilon_0\, x` consistent
with a uniform axial strain field, must recover
:math:`\varepsilon = \varepsilon_0` everywhere to **machine
precision**.

Reference: Irons, B. M. and Razzaque, A. (1972) "Experience with
the patch test for convergence of finite elements," *The
Mathematical Foundations of the Finite Element Method*, pp.
557-587.

A discretisation that fails this test cannot pass the inf-sup /
consistency conditions and therefore cannot converge at the
optimal rate.  A pass at machine precision says the element
kernel is wiring :math:`B`, :math:`D`, and the isoparametric
mapping correctly under non-trivial geometry.

Vendor cross-references
-----------------------

======================================  =====================  ============================================
Source                                   Expected               Problem ID / location
======================================  =====================  ============================================
Irons-Razzaque 1972 (original)           identity (exact)       Math Foundations of FEM, pp. 557-587
Taylor-Beresford-Wilson 1976             identity (exact)       IJNME 10 (1976), pp. 1211-1219
Hughes, The Finite Element Method        identity (exact)       Dover (2000), §4.4.3
NAFEMS Background to Benchmarks          identity (exact)       BtB-4.2 (patch-test completeness)
Abaqus Verification Manual               identity (exact)       AVM 1.4.4 (C3D8 patch test)
======================================  =====================  ============================================
"""

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation.problems import IronsPatchTest

# %%
# Problem statement

problem = IronsPatchTest()
pv_ = problem.published_values[0]
print(f"{problem.name}: {problem.description}")
print(f"\n  source : {pv_.source}")
print(f"  formula: {pv_.formula}")
print(f"  target : {pv_.value} (tolerance {pv_.tolerance:.0e})")

# %%
# Run the check — single-shot ``validate()``

results = problem.validate()
r = results[0]
print(f"\n  computed max_strain_error: {r.computed:.3e}")
print(f"  passes tolerance          : {r.passed}")

# %%
# Figure — the distorted patch
# ----------------------------
# The validation's ``build_model`` offsets the single interior
# node by 10 % of the edge length so the isoparametric mapping
# is non-trivial; rendering the resulting mesh is a nice sanity
# check that the geometry actually has the promised distortion.

m = problem.build_model()
plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(m.grid, show_edges=True, opacity=0.7, color="#b8cde3")
plotter.add_points(
    np.asarray(m.grid.points),
    render_points_as_spheres=True,
    point_size=12,
    color="red",
)
plotter.view_isometric()
plotter.camera.zoom(1.1)
plotter.show()

# %%
# Acceptance
# ----------
# The 1e-9 tolerance is conservative — we measure ~1e-19 on
# this build.  A regression above 1e-12 would indicate a
# precision-dropping change in the assembly / solve pipeline.

assert r.passed, f"patch test failed: {r.computed:.3e} > {pv_.tolerance:.1e}"