diff --git a/include/fesa/elements/element.h b/include/fesa/elements/element.h
index d8ff473..c0a02d7 100644
--- a/include/fesa/elements/element.h
+++ b/include/fesa/elements/element.h
@@ -11,6 +11,7 @@
 #include "fesa/math/matrix.h"
 #include "fesa/math/vector.h"
 #include "fesa/results/result_records.h"
+#include "fesa/results/solid_result_records.h"
 
 namespace fesa {
 
@@ -52,7 +53,8 @@ struct ShellElementResultRows {
 
 /// @brief Selects the physical recovery shape of one runtime element.
 using ElementResultPayload =
-    std::variant<BeamElementResultRows, ShellElementResultRows>;
+    std::variant<BeamElementResultRows, ShellElementResultRows,
+                 SolidElementResultRows>;
 
 /// @brief Carries stable source identity with one typed recovery payload.
 struct ElementResultBundle {
diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h
new file mode 100644
index 0000000..29e6866
--- /dev/null
+++ b/include/fesa/elements/hex8_element.h
@@ -0,0 +1,83 @@
+#ifndef FESA_ELEMENTS_HEX8_ELEMENT_H_
+#define FESA_ELEMENTS_HEX8_ELEMENT_H_
+
+#include <array>
+#include <optional>
+
+#include "fesa/elements/element.h"
+#include "fesa/materials/isotropic_linear_elastic_material.h"
+#include "fesa/results/solid_result_records.h"
+
+namespace fesa {
+
+/// @brief Owns all local quantities from one validated linear or TL trial.
+struct Hex8Response {
+  Vector internal_force{24U};
+  Matrix material_tangent{24U, 24U};
+  Matrix geometric_tangent{24U, 24U};
+  Matrix tangent{24U, 24U};
+  double strain_energy{0.0};
+  std::array<SolidIntegrationPointRow, 8> points{};
+};
+
+/// @brief Implements full-integration C3D8 linear and TL StVK mechanics.
+/// @note Reference geometry is owned and immutable; no trial is committed.
+class Hex8Element final : public Element {
+ public:
+  using Coordinates = std::array<std::array<double, 3>, 8>;
+
+  /// @brief Validates identity, isotropic material and whole reference cell.
+  /// @param coordinates Reference global coordinates in source node order.
+  /// @param density Optional reference mass per volume, positive if present.
+  /// @note All inputs are copied; no Domain or section lifetime is required.
+  static Result<Hex8Element> Create(
+      EntityIndex element_index, ElementDofLayout layout,
+      Coordinates coordinates, const LinearElasticMaterial& material,
+      std::optional<double> density = std::nullopt,
+      SourceLocation location = {});
+
+  const ElementDofLayout& DofLayout() const noexcept override;
+  Result<ElementStiffnessContribution> ComputeStiffness() const override;
+  Result<ElementResultBundle> Recover(
+      const Vector& displacement) const override;
+
+  /// @brief Returns f=K0*u and infinitesimal tensor strain at eight points.
+  Result<Hex8Response> EvaluateLinear(const Vector& displacement) const;
+  /// @brief Returns one certified StVK trial and its analytic energy Hessian.
+  Result<Hex8Response> EvaluateNonlinear(const Vector& displacement) const;
+  /// @brief Returns the constant reference consistent mass; density required.
+  Result<ElementStiffnessContribution> ComputeMass() const;
+
+ private:
+  /// @brief Caches shape, reference gradients and positive volume weights.
+  struct ReferencePoint {
+    std::array<double, 8> shape{};
+    std::array<std::array<double, 3>, 8> gradients{};
+    std::array<double, 3> natural{};
+    std::array<double, 3> position{};
+    double weight{0.0};
+  };
+  /// @brief Constructs an owned candidate before reference validation.
+  Hex8Element(EntityIndex element_index, ElementDofLayout layout,
+              Coordinates coordinates, double lambda, double mu,
+              std::optional<double> density, SourceLocation location);
+  /// @brief Evaluates without mutating the cached reference or previous values.
+  Result<Hex8Response> Evaluate(const Vector& displacement,
+                                bool nonlinear) const;
+  /// @brief Attaches stable source context to a model failure.
+  Status Failure(const char* message) const;
+
+  EntityIndex element_index_;
+  ElementDofLayout layout_;
+  Coordinates coordinates_;
+  double lambda_;
+  double mu_;
+  std::optional<double> density_;
+  SourceLocation location_;
+  std::array<ReferencePoint, 8> points_{};
+  Matrix stiffness_{24U, 24U};
+};
+
+}  // namespace fesa
+
+#endif  // FESA_ELEMENTS_HEX8_ELEMENT_H_
diff --git a/include/fesa/results/solid_result_records.h b/include/fesa/results/solid_result_records.h
new file mode 100644
index 0000000..e5028ac
--- /dev/null
+++ b/include/fesa/results/solid_result_records.h
@@ -0,0 +1,31 @@
+#ifndef FESA_RESULTS_SOLID_RESULT_RECORDS_H_
+#define FESA_RESULTS_SOLID_RESULT_RECORDS_H_
+
+#include <array>
+#include <vector>
+
+#include "fesa/core/source_identity.h"
+
+namespace fesa {
+
+/// @brief Stores one unaveraged global solid integration-point result.
+/// @note Tensor order is [11,22,33,12,13,23], with tensor shear strain.
+struct SolidIntegrationPointRow {
+  EntityIndex element;
+  int integration_point;
+  std::array<double, 3> natural_coordinates;
+  std::array<double, 3> reference_coordinates;
+  std::array<double, 6> cauchy_stress;
+  // Infinitesimal strain for linear response; reference Green strain for TL.
+  std::array<double, 6> strain;
+};
+
+/// @brief Keeps eight physical solid rows and their reference-volume energy.
+struct SolidElementResultRows {
+  std::vector<SolidIntegrationPointRow> rows;
+  double physical_strain_energy{0.0};
+};
+
+}  // namespace fesa
+
+#endif  // FESA_RESULTS_SOLID_RESULT_RECORDS_H_
diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp
new file mode 100644
index 0000000..6a4dfe9
--- /dev/null
+++ b/src/fesa/elements/hex8_element.cpp
@@ -0,0 +1,325 @@
+#include "fesa/elements/hex8_element.h"
+
+#include <cmath>
+#include <set>
+#include <utility>
+
+#include "hex8_geometry.h"
+
+namespace fesa {
+namespace {
+using hex8_internal::Tensor;
+using hex8_internal::Triple;
+/// @brief Checks every tensor component before recovery or inverse use.
+bool Finite(const Tensor& value) {
+  for (const auto& row : value)
+    for (double entry : row)
+      if (!std::isfinite(entry)) return false;
+  return true;
+}
+/// @brief Packs actual tensor shear, never engineering shear, in schema order.
+std::array<double, 6> Pack(const Tensor& value) {
+  return {value[0][0], value[1][1], value[2][2],
+          value[0][1], value[0][2], value[1][2]};
+}
+/// @brief Computes a Euclidean contraction in deterministic component order.
+double Dot(const Triple& a, const Triple& b) {
+  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+/// @brief Gives the finite-control failure without accepting a geometry guess.
+const char* CertificateMessage(hex8_internal::Certificate certificate) {
+  switch (certificate) {
+    case hex8_internal::Certificate::kNonpositive:
+      return "Whole-cell Jacobian has a nonpositive witness.";
+    case hex8_internal::Certificate::kUncertain:
+      return "Whole-cell Jacobian positivity is unresolved at depth 10 or cell "
+             "budget 65536.";
+    case hex8_internal::Certificate::kNonfinite:
+      return "Whole-cell Jacobian interval arithmetic is nonfinite.";
+    case hex8_internal::Certificate::kPositive:
+      break;
+  }
+  return "Invalid Jacobian certificate.";
+}
+}  // namespace
+
+Hex8Element::Hex8Element(EntityIndex index, ElementDofLayout layout,
+                         Coordinates coordinates, double lambda, double mu,
+                         std::optional<double> density, SourceLocation location)
+    : element_index_{index},
+      layout_{std::move(layout)},
+      coordinates_{coordinates},
+      lambda_{lambda},
+      mu_{mu},
+      density_{density},
+      location_{std::move(location)} {}
+
+Status Hex8Element::Failure(const char* message) const {
+  return Status::Failure(FailureCategory::kModel,
+                         {{Severity::kError, "invalid-hex8", location_, "C3D8",
+                           layout_.source_id.instance_name + ":" +
+                               layout_.source_id.source_label_text,
+                           message}});
+}
+
+Result<Hex8Element> Hex8Element::Create(EntityIndex index,
+                                        ElementDofLayout layout,
+                                        Coordinates coordinates,
+                                        const LinearElasticMaterial& material,
+                                        std::optional<double> density,
+                                        SourceLocation location) {
+  const double e = material.YoungsModulus();
+  const double nu = material.PoissonsRatio();
+  Hex8Element candidate{index, std::move(layout), coordinates,        0.0,
+                        0.0,   density,           std::move(location)};
+  const std::vector<DofComponent> components{
+      DofComponent::kUx, DofComponent::kUy, DofComponent::kUz};
+  const std::set<EntityIndex> nodes(candidate.layout_.node_indices.begin(),
+                                    candidate.layout_.node_indices.end());
+  if (candidate.layout_.node_indices.size() != 8U || nodes.size() != 8U ||
+      candidate.layout_.components_per_node != components)
+    return Result<Hex8Element>::Failure(candidate.Failure(
+        "Hex8 requires eight distinct nodes with exactly UX, UY, UZ."));
+  if (!std::isfinite(e) || e <= 0.0 || !std::isfinite(nu) || nu <= -1.0 ||
+      nu >= 0.5)
+    return Result<Hex8Element>::Failure(
+        candidate.Failure("Hex8 requires finite E>0 and -1<nu<0.5."));
+  candidate.mu_ = e / (2.0 * (1.0 + nu));
+  candidate.lambda_ = e * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
+  if (!std::isfinite(candidate.mu_) || candidate.mu_ <= 0.0 ||
+      !std::isfinite(candidate.lambda_) ||
+      !std::isfinite(candidate.lambda_ + 2.0 * candidate.mu_))
+    return Result<Hex8Element>::Failure(
+        candidate.Failure("Hex8 Lame moduli must remain finite with "
+                          "representable positive shear modulus."));
+  if (density && (!std::isfinite(*density) || *density <= 0.0))
+    return Result<Hex8Element>::Failure(candidate.Failure(
+        "Hex8 reference density must be finite and positive."));
+  const auto certificate = hex8_internal::Certify(coordinates);
+  if (certificate != hex8_internal::Certificate::kPositive)
+    return Result<Hex8Element>::Failure(
+        candidate.Failure(CertificateMessage(certificate)));
+  const double gauss = 1.0 / std::sqrt(3.0);
+  for (std::size_t g = 0; g < 8; ++g) {
+    auto& point = candidate.points_[g];
+    point.natural = {g % 2U == 0U ? -gauss : gauss,
+                     (g / 2U) % 2U == 0U ? -gauss : gauss,
+                     g < 4U ? -gauss : gauss};
+    std::array<Triple, 8> derivatives{};
+    hex8_internal::Shape(point.natural, point.shape, derivatives);
+    const auto jacobian = hex8_internal::Jacobian(coordinates, derivatives);
+    point.weight = hex8_internal::Determinant(jacobian);
+    if (!Finite(jacobian) || !std::isfinite(point.weight) ||
+        point.weight <= 0.0)
+      return Result<Hex8Element>::Failure(candidate.Failure(
+          "Reference quadrature Jacobian must be finite and positive."));
+    const auto inverse = hex8_internal::Inverse(jacobian, point.weight);
+    if (!Finite(inverse))
+      return Result<Hex8Element>::Failure(
+          candidate.Failure("Reference Jacobian inverse is nonfinite."));
+    point.position = coordinates[0];
+    for (std::size_t a = 0; a < 8; ++a)
+      for (std::size_t i = 0; i < 3; ++i) {
+        point.position[i] +=
+            point.shape[a] * (coordinates[a][i] - coordinates[0][i]);
+        for (std::size_t d = 0; d < 3; ++d)
+          point.gradients[a][i] += inverse[d][i] * derivatives[a][d];
+        if (!std::isfinite(point.gradients[a][i]) ||
+            !std::isfinite(point.position[i]))
+          return Result<Hex8Element>::Failure(candidate.Failure(
+              "Reference gradient or point coordinate is nonfinite."));
+      }
+    for (std::size_t a = 0; a < 8; ++a)
+      for (std::size_t b = 0; b < 8; ++b) {
+        const auto& ga = point.gradients[a];
+        const auto& gb = point.gradients[b];
+        const double dot = Dot(ga, gb);
+        for (std::size_t i = 0; i < 3; ++i)
+          for (std::size_t k = 0; k < 3; ++k) {
+            candidate.stiffness_(3U * a + i, 3U * b + k) +=
+                point.weight * (candidate.lambda_ * ga[i] * gb[k] +
+                                candidate.mu_ * ga[k] * gb[i] +
+                                (i == k ? candidate.mu_ * dot : 0.0));
+          }
+      }
+  }
+  for (std::size_t i = 0; i < 24; ++i)
+    for (std::size_t j = 0; j < 24; ++j)
+      if (!std::isfinite(candidate.stiffness_(i, j)))
+        return Result<Hex8Element>::Failure(
+            candidate.Failure("Hex8 reference stiffness is nonfinite."));
+  for (std::size_t i = 0; i < 24; ++i)
+    if (candidate.stiffness_(i, i) <= 0.0)
+      return Result<Hex8Element>::Failure(
+          candidate.Failure("Hex8 reference stiffness has an unrepresentable "
+                            "positive diagonal."));
+  return Result<Hex8Element>::Success(std::move(candidate));
+}
+
+const ElementDofLayout& Hex8Element::DofLayout() const noexcept {
+  return layout_;
+}
+Result<ElementStiffnessContribution> Hex8Element::ComputeStiffness() const {
+  return Result<ElementStiffnessContribution>::Success({layout_, stiffness_});
+}
+Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
+  if (!density_)
+    return Result<ElementStiffnessContribution>::Failure(
+        Failure("Hex8 mass requires reference density."));
+  Matrix mass{24U, 24U};
+  for (const auto& point : points_)
+    for (std::size_t a = 0; a < 8; ++a)
+      for (std::size_t b = 0; b < 8; ++b) {
+        const double value =
+            *density_ * point.shape[a] * point.shape[b] * point.weight;
+        if (!std::isfinite(value))
+          return Result<ElementStiffnessContribution>::Failure(
+              Failure("Hex8 reference mass is nonfinite."));
+        for (std::size_t i = 0; i < 3; ++i)
+          mass(3U * a + i, 3U * b + i) += value;
+      }
+  for (std::size_t i = 0; i < 24; ++i)
+    for (std::size_t j = 0; j < 24; ++j)
+      if (!std::isfinite(mass(i, j)))
+        return Result<ElementStiffnessContribution>::Failure(
+            Failure("Hex8 accumulated reference mass is nonfinite."));
+  for (std::size_t i = 0; i < 24; ++i)
+    if (mass(i, i) <= 0.0)
+      return Result<ElementStiffnessContribution>::Failure(Failure(
+          "Hex8 reference mass has an unrepresentable positive diagonal."));
+  return Result<ElementStiffnessContribution>::Success(
+      {layout_, std::move(mass)});
+}
+Result<ElementResultBundle> Hex8Element::Recover(const Vector& u) const {
+  auto response = EvaluateLinear(u);
+  if (!response.HasValue())
+    return Result<ElementResultBundle>::Failure(response.GetStatus());
+  SolidElementResultRows rows;
+  rows.rows.assign(response.Value().points.begin(),
+                   response.Value().points.end());
+  rows.physical_strain_energy = response.Value().strain_energy;
+  return Result<ElementResultBundle>::Success(
+      {layout_.source_id, std::move(rows)});
+}
+Result<Hex8Response> Hex8Element::EvaluateLinear(const Vector& u) const {
+  return Evaluate(u, false);
+}
+Result<Hex8Response> Hex8Element::EvaluateNonlinear(const Vector& u) const {
+  return Evaluate(u, true);
+}
+Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
+                                           bool nonlinear) const {
+  if (u.Size() != 24U)
+    return Result<Hex8Response>::Failure(
+        Failure("Hex8 trial requires 24 translations."));
+  const auto certificate = hex8_internal::Certify(coordinates_, &u);
+  if (certificate != hex8_internal::Certificate::kPositive)
+    return Result<Hex8Response>::Failure(
+        Failure(CertificateMessage(certificate)));
+  Hex8Response response;
+  if (!nonlinear) {
+    response.material_tangent = stiffness_;
+    response.tangent = stiffness_;
+    response.internal_force = stiffness_.Multiply(u);
+    response.strain_energy = 0.5 * u.Dot(response.internal_force);
+  }
+  for (std::size_t g = 0; g < 8; ++g) {
+    const auto& point = points_[g];
+    Tensor gradient{};
+    for (std::size_t a = 0; a < 8; ++a)
+      for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+          gradient[i][j] += (u[3U * a + i] - u[i]) * point.gradients[a][j];
+    Tensor f = gradient;
+    for (std::size_t i = 0; i < 3; ++i) f[i][i] += 1.0;
+    Tensor strain{};
+    for (std::size_t i = 0; i < 3; ++i)
+      for (std::size_t j = 0; j < 3; ++j) {
+        if (nonlinear) {
+          for (std::size_t k = 0; k < 3; ++k)
+            strain[i][j] += 0.5 * f[k][i] * f[k][j];
+          if (i == j) strain[i][j] -= 0.5;
+        } else
+          strain[i][j] = 0.5 * (gradient[i][j] + gradient[j][i]);
+      }
+    const double trace = strain[0][0] + strain[1][1] + strain[2][2];
+    Tensor s{};
+    for (std::size_t i = 0; i < 3; ++i)
+      for (std::size_t j = 0; j < 3; ++j)
+        s[i][j] = 2.0 * mu_ * strain[i][j] + (i == j ? lambda_ * trace : 0.0);
+    Tensor stress = s;
+    if (nonlinear) {
+      const double det_f = hex8_internal::Determinant(f);
+      if (!std::isfinite(det_f) || det_f <= 0.0)
+        return Result<Hex8Response>::Failure(Failure(
+            "Hex8 trial deformation determinant must be finite and positive."));
+      Tensor p{};
+      Tensor fft{};
+      double energy = 0.5 * lambda_ * trace * trace;
+      for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j) {
+          energy += mu_ * strain[i][j] * strain[i][j];
+          for (std::size_t k = 0; k < 3; ++k) {
+            p[i][j] += f[i][k] * s[k][j];
+            fft[i][j] += f[i][k] * f[j][k];
+          }
+        }
+      response.strain_energy += energy * point.weight;
+      stress = {};
+      for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j)
+          for (std::size_t k = 0; k < 3; ++k)
+            stress[i][j] += p[i][k] * f[j][k] / det_f;
+      for (std::size_t a = 0; a < 8; ++a) {
+        const auto& ga = point.gradients[a];
+        for (std::size_t i = 0; i < 3; ++i)
+          response.internal_force[3U * a + i] += Dot(p[i], ga) * point.weight;
+        for (std::size_t b = 0; b < 8; ++b) {
+          const auto& gb = point.gradients[b];
+          double geometric = 0.0;
+          for (std::size_t j = 0; j < 3; ++j)
+            for (std::size_t l = 0; l < 3; ++l)
+              geometric += ga[j] * s[l][j] * gb[l];
+          for (std::size_t i = 0; i < 3; ++i)
+            for (std::size_t k = 0; k < 3; ++k) {
+              const auto row = 3U * a + i;
+              const auto column = 3U * b + k;
+              response.material_tangent(row, column) +=
+                  point.weight * (lambda_ * Dot(f[i], ga) * Dot(f[k], gb) +
+                                  mu_ * Dot(f[i], gb) * Dot(f[k], ga) +
+                                  mu_ * fft[i][k] * Dot(ga, gb));
+              if (i == k)
+                response.geometric_tangent(row, column) +=
+                    geometric * point.weight;
+            }
+        }
+      }
+    }
+    if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(stress))
+      return Result<Hex8Response>::Failure(
+          Failure("Hex8 trial tensor is nonfinite."));
+    response.points[g] = {element_index_, static_cast<int>(g + 1U),
+                          point.natural,  point.position,
+                          Pack(stress),   Pack(strain)};
+  }
+  if (!std::isfinite(response.strain_energy))
+    return Result<Hex8Response>::Failure(
+        Failure("Hex8 trial energy is nonfinite."));
+  for (std::size_t i = 0; i < 24; ++i) {
+    if (!std::isfinite(response.internal_force[i]))
+      return Result<Hex8Response>::Failure(
+          Failure("Hex8 trial force is nonfinite."));
+    for (std::size_t j = 0; j < 24; ++j) {
+      if (nonlinear)
+        response.tangent(i, j) =
+            response.material_tangent(i, j) + response.geometric_tangent(i, j);
+      if (!std::isfinite(response.tangent(i, j)) ||
+          !std::isfinite(response.material_tangent(i, j)) ||
+          !std::isfinite(response.geometric_tangent(i, j)))
+        return Result<Hex8Response>::Failure(
+            Failure("Hex8 trial tangent is nonfinite."));
+    }
+  }
+  return Result<Hex8Response>::Success(std::move(response));
+}
+}  // namespace fesa
diff --git a/src/fesa/elements/hex8_geometry.cpp b/src/fesa/elements/hex8_geometry.cpp
new file mode 100644
index 0000000..e758daf
--- /dev/null
+++ b/src/fesa/elements/hex8_geometry.cpp
@@ -0,0 +1,239 @@
+#include "hex8_geometry.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <utility>
+#include <vector>
+
+namespace fesa::hex8_internal {
+namespace {
+/// @brief Encloses exact-real arithmetic on the supplied binary64 inputs.
+struct Interval {
+  double lower{0.0};
+  double upper{0.0};
+};
+using Coefficients = std::array<Interval, 27>;
+
+/// @brief Rounds each primitive outward by one representable binary64 value.
+/// @note /fp:precise (repository default) forbids reassociation; nextafter also
+/// encloses subnormal rounding. Overflow/nonfinite bounds fail certification.
+Interval Outward(double lower, double upper) {
+  if (!std::isfinite(lower) || !std::isfinite(upper)) {
+    const double invalid = std::numeric_limits<double>::quiet_NaN();
+    return {invalid, invalid};
+  }
+  return {std::nextafter(lower, -std::numeric_limits<double>::infinity()),
+          std::nextafter(upper, std::numeric_limits<double>::infinity())};
+}
+/// @brief Recognizes an exact zero, for exact identity operations only.
+bool Zero(Interval a) { return a.lower == 0.0 && a.upper == 0.0; }
+/// @brief Encloses addition, including the current X+u construction.
+Interval Add(Interval a, Interval b) {
+  if (Zero(a)) return b;
+  if (Zero(b)) return a;
+  return Outward(a.lower + b.lower, a.upper + b.upper);
+}
+/// @brief Encloses subtraction without assuming correlated interval operands.
+Interval Subtract(Interval a, Interval b) {
+  if (Zero(b)) return a;
+  if (a.lower == a.upper && b.lower == b.upper && a.lower == b.lower) return {};
+  return Outward(a.lower - b.upper, a.upper - b.lower);
+}
+/// @brief Encloses multiplication across all endpoint sign combinations.
+Interval Multiply(Interval a, Interval b) {
+  if (Zero(a) || Zero(b)) return {};
+  const std::array<double, 4> products{a.lower * b.lower, a.lower * b.upper,
+                                       a.upper * b.lower, a.upper * b.upper};
+  for (double product : products)
+    if (!std::isfinite(product)) {
+      const double invalid = std::numeric_limits<double>::quiet_NaN();
+      return {invalid, invalid};
+    }
+  return Outward(*std::min_element(products.begin(), products.end()),
+                 *std::max_element(products.begin(), products.end()));
+}
+/// @brief Multiplies by an exact binary64 constant with outward rounding.
+Interval Scale(Interval a, double value) { return Multiply(a, {value, value}); }
+/// @brief Detects an invalid enclosure before any positivity decision.
+bool Finite(Interval a) {
+  return std::isfinite(a.lower) && std::isfinite(a.upper) && a.lower <= a.upper;
+}
+/// @brief Preserves xi-fastest tensor ordering.
+std::size_t Index(std::size_t i, std::size_t j, std::size_t k) {
+  return i + 3U * j + 9U * k;
+}
+/// @brief Evaluates determinant samples with every primitive enclosed.
+Interval Sample(const std::array<std::array<Interval, 3>, 8>& x,
+                const Triple& natural) {
+  std::array<std::array<Interval, 3>, 3> jacobian{};
+  for (std::size_t a = 0; a < 8; ++a) {
+    for (std::size_t d = 0; d < 3; ++d) {
+      Interval derivative{0.125 * kNodeSigns[a][d], 0.125 * kNodeSigns[a][d]};
+      for (std::size_t axis = 0; axis < 3; ++axis) {
+        if (axis != d) {
+          const auto factor =
+              Add({1.0, 1.0},
+                  Scale({natural[axis], natural[axis]}, kNodeSigns[a][axis]));
+          derivative = Multiply(derivative, factor);
+        }
+      }
+      for (std::size_t i = 0; i < 3; ++i)
+        jacobian[i][d] = Add(jacobian[i][d],
+                             Multiply(Subtract(x[a][i], x[0][i]), derivative));
+    }
+  }
+  const auto& j = jacobian;
+  return Add(Subtract(Multiply(j[0][0], Subtract(Multiply(j[1][1], j[2][2]),
+                                                 Multiply(j[1][2], j[2][1]))),
+                      Multiply(j[0][1], Subtract(Multiply(j[1][0], j[2][2]),
+                                                 Multiply(j[1][2], j[2][0])))),
+             Multiply(j[0][2], Subtract(Multiply(j[1][0], j[2][1]),
+                                        Multiply(j[1][1], j[2][0]))));
+}
+/// @brief Applies the degree-two nodal-to-Bernstein transform on one axis.
+void Transform(Coefficients& coefficients, std::size_t axis) {
+  const auto original = coefficients;
+  const std::size_t stride = axis == 0U ? 1U : (axis == 1U ? 3U : 9U);
+  for (std::size_t k = 0; k < 3; ++k)
+    for (std::size_t j = 0; j < 3; ++j)
+      for (std::size_t i = 0; i < 3; ++i) {
+        const std::array<std::size_t, 3> coordinate{i, j, k};
+        if (coordinate[axis] != 0U) continue;
+        const auto start = Index(i, j, k);
+        coefficients[start + stride] = Subtract(
+            Scale(original[start + stride], 2.0),
+            Scale(Add(original[start], original[start + 2U * stride]), 0.5));
+      }
+}
+/// @brief Uses interval de Casteljau midpoint subdivision along one axis.
+std::pair<Coefficients, Coefficients> Split(const Coefficients& original,
+                                            std::size_t axis) {
+  Coefficients left = original;
+  Coefficients right = original;
+  const std::size_t stride = axis == 0U ? 1U : (axis == 1U ? 3U : 9U);
+  for (std::size_t k = 0; k < 3; ++k)
+    for (std::size_t j = 0; j < 3; ++j)
+      for (std::size_t i = 0; i < 3; ++i) {
+        const std::array<std::size_t, 3> coordinate{i, j, k};
+        if (coordinate[axis] != 0U) continue;
+        const auto start = Index(i, j, k);
+        const auto a =
+            Scale(Add(original[start], original[start + stride]), 0.5);
+        const auto b = Scale(
+            Add(original[start + stride], original[start + 2U * stride]), 0.5);
+        const auto middle = Scale(Add(a, b), 0.5);
+        left[start + stride] = a;
+        left[start + 2U * stride] = middle;
+        right[start] = middle;
+        right[start + stride] = b;
+      }
+  return {left, right};
+}
+}  // namespace
+
+Certificate Certify(const Hex8Element::Coordinates& coordinates,
+                    const Vector* displacement) {
+  std::array<std::array<Interval, 3>, 8> x{};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i) {
+      const auto value = coordinates[a][i];
+      if (!std::isfinite(value)) return Certificate::kNonfinite;
+      x[a][i] = {value, value};
+      if (displacement) {
+        const auto u = (*displacement)[3U * a + i];
+        if (!std::isfinite(u)) return Certificate::kNonfinite;
+        x[a][i] = Add(x[a][i], {u, u});
+      }
+      if (!Finite(x[a][i])) return Certificate::kNonfinite;
+    }
+  Coefficients coefficients{};
+  for (std::size_t k = 0; k < 3; ++k)
+    for (std::size_t j = 0; j < 3; ++j)
+      for (std::size_t i = 0; i < 3; ++i) {
+        const auto sample = Sample(
+            x, {static_cast<double>(i) - 1.0, static_cast<double>(j) - 1.0,
+                static_cast<double>(k) - 1.0});
+        if (!Finite(sample)) return Certificate::kNonfinite;
+        if (sample.upper <= 0.0) return Certificate::kNonpositive;
+        coefficients[Index(i, j, k)] = sample;
+      }
+  for (std::size_t axis = 0; axis < 3; ++axis) Transform(coefficients, axis);
+  struct Cell {
+    Coefficients coefficients;
+    unsigned depth;
+  };
+  std::vector<Cell> pending{{coefficients, 0U}};
+  std::size_t processed = 0U;
+  while (!pending.empty()) {
+    if (++processed > 65536U) return Certificate::kUncertain;
+    Cell cell = std::move(pending.back());
+    pending.pop_back();
+    bool positive = true;
+    for (const auto c : cell.coefficients) {
+      if (!Finite(c)) return Certificate::kNonfinite;
+      positive = positive && c.lower > 0.0;
+    }
+    if (positive) continue;
+    // Corner coefficients are exact polynomial values enclosed by intervals.
+    for (std::size_t k = 0; k < 3; k += 2)
+      for (std::size_t j = 0; j < 3; j += 2)
+        for (std::size_t i = 0; i < 3; i += 2)
+          if (cell.coefficients[Index(i, j, k)].upper <= 0.0)
+            return Certificate::kNonpositive;
+    if (cell.depth >= 10U) return Certificate::kUncertain;
+    const auto xi = Split(cell.coefficients, 0U);
+    for (const auto& xchild : {xi.first, xi.second}) {
+      const auto eta = Split(xchild, 1U);
+      for (const auto& ychild : {eta.first, eta.second}) {
+        const auto zeta = Split(ychild, 2U);
+        pending.push_back({zeta.second, cell.depth + 1U});
+        pending.push_back({zeta.first, cell.depth + 1U});
+      }
+    }
+  }
+  return Certificate::kPositive;
+}
+
+void Shape(const Triple& natural, std::array<double, 8>& values,
+           std::array<Triple, 8>& derivatives) {
+  for (std::size_t a = 0; a < 8; ++a) {
+    values[a] = 0.125;
+    for (std::size_t d = 0; d < 3; ++d)
+      values[a] *= 1.0 + kNodeSigns[a][d] * natural[d];
+    for (std::size_t d = 0; d < 3; ++d) {
+      derivatives[a][d] = 0.125 * kNodeSigns[a][d];
+      for (std::size_t axis = 0; axis < 3; ++axis)
+        if (axis != d)
+          derivatives[a][d] *= 1.0 + kNodeSigns[a][axis] * natural[axis];
+    }
+  }
+}
+Tensor Jacobian(const Hex8Element::Coordinates& coordinates,
+                const std::array<Triple, 8>& derivatives) {
+  Tensor j{};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i)
+      for (std::size_t d = 0; d < 3; ++d)
+        j[i][d] += (coordinates[a][i] - coordinates[0][i]) * derivatives[a][d];
+  return j;
+}
+double Determinant(const Tensor& j) {
+  return j[0][0] * (j[1][1] * j[2][2] - j[1][2] * j[2][1]) -
+         j[0][1] * (j[1][0] * j[2][2] - j[1][2] * j[2][0]) +
+         j[0][2] * (j[1][0] * j[2][1] - j[1][1] * j[2][0]);
+}
+Tensor Inverse(const Tensor& j, double determinant) {
+  Tensor inverse{};
+  for (std::size_t i = 0; i < 3; ++i)
+    for (std::size_t d = 0; d < 3; ++d) {
+      const auto r1 = (d + 1U) % 3U;
+      const auto r2 = (d + 2U) % 3U;
+      const auto c1 = (i + 1U) % 3U;
+      const auto c2 = (i + 2U) % 3U;
+      inverse[i][d] =
+          (j[r1][c1] * j[r2][c2] - j[r1][c2] * j[r2][c1]) / determinant;
+    }
+  return inverse;
+}
+}  // namespace fesa::hex8_internal
diff --git a/src/fesa/elements/hex8_geometry.h b/src/fesa/elements/hex8_geometry.h
new file mode 100644
index 0000000..54a8a44
--- /dev/null
+++ b/src/fesa/elements/hex8_geometry.h
@@ -0,0 +1,32 @@
+#ifndef FESA_ELEMENTS_HEX8_GEOMETRY_H_
+#define FESA_ELEMENTS_HEX8_GEOMETRY_H_
+
+#include <array>
+
+#include "fesa/elements/hex8_element.h"
+
+namespace fesa::hex8_internal {
+using Triple = std::array<double, 3>;
+using Tensor = std::array<Triple, 3>;
+inline constexpr int kNodeSigns[8][3] = {{-1, -1, -1}, {1, -1, -1}, {1, 1, -1},
+                                         {-1, 1, -1},  {-1, -1, 1}, {1, -1, 1},
+                                         {1, 1, 1},    {-1, 1, 1}};
+/// @brief Distinguishes a witness from an exhausted conservative certificate.
+enum class Certificate { kPositive, kNonpositive, kUncertain, kNonfinite };
+/// @brief Certifies the closed natural cube with interval Bernstein bounds.
+/// @note Current addition rounding is enclosed when displacement is supplied.
+Certificate Certify(const Hex8Element::Coordinates& coordinates,
+                    const Vector* displacement = nullptr);
+/// @brief Evaluates trilinear shapes and natural derivatives in source order.
+void Shape(const Triple& natural, std::array<double, 8>& values,
+           std::array<Triple, 8>& derivatives);
+/// @brief Computes a three-column map Jacobian using translated coordinates.
+Tensor Jacobian(const Hex8Element::Coordinates& coordinates,
+                const std::array<Triple, 8>& derivatives);
+/// @brief Computes the signed determinant; no absolute value or clamp.
+double Determinant(const Tensor& matrix);
+/// @brief Computes the inverse for an already validated positive determinant.
+Tensor Inverse(const Tensor& matrix, double determinant);
+}  // namespace fesa::hex8_internal
+
+#endif  // FESA_ELEMENTS_HEX8_GEOMETRY_H_
diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp
index 854c79d..b23041a 100644
--- a/src/fesa/results/result_recovery.cpp
+++ b/src/fesa/results/result_recovery.cpp
@@ -356,10 +356,16 @@ Status ResultRecovery::Recover(const AnalysisModel& model,
             return results_internal::BeamResultRecovery(
                 {domain.SourcePath(), 0U}, recovered.Value().source_id,
                 element_index, expected_nodes, payload, candidate);
-          } else {
+          } else if constexpr (std::is_same_v<Payload, ShellElementResultRows>) {
             return results_internal::ShellResultRecovery(
                 {domain.SourcePath(), 0U}, recovered.Value().source_id,
                 element_index, payload, candidate);
+          } else {
+            return RecoveryFailure(
+                "unsupported-solid-result-recovery",
+                {domain.SourcePath(), 0U},
+                recovered.Value().source_id.source_label_text,
+                "Solid results require the solid procedure recovery owner.");
           }
         },
         recovered.Value().payload);
diff --git a/tests/unit/elements/hex8_element_test.cpp b/tests/unit/elements/hex8_element_test.cpp
new file mode 100644
index 0000000..5206693
--- /dev/null
+++ b/tests/unit/elements/hex8_element_test.cpp
@@ -0,0 +1,617 @@
+#include "fesa/elements/hex8_element.h"
+
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <optional>
+#include <string>
+
+namespace {
+using Coordinates = fesa::Hex8Element::Coordinates;
+constexpr int kSigns[8][3] = {{-1, -1, -1}, {1, -1, -1}, {1, 1, -1},
+                              {-1, 1, -1},  {-1, -1, 1}, {1, -1, 1},
+                              {1, 1, 1},    {-1, 1, 1}};
+Coordinates Brick() {
+  Coordinates x{};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i)
+      x[a][i] = static_cast<double>(kSigns[a][i]) * static_cast<double>(i + 1U);
+  return x;
+}
+fesa::Result<fesa::Hex8Element> Make(Coordinates x = Brick(), double e = 120.0,
+                                     double nu = 0.25,
+                                     std::optional<double> rho = 2.0) {
+  fesa::ElementDofLayout layout{
+      {"PART-1", 42, "42"},
+      {0, 1, 2, 3, 4, 5, 6, 7},
+      {fesa::DofComponent::kUx, fesa::DofComponent::kUy,
+       fesa::DofComponent::kUz}};
+  fesa::LinearElasticMaterial material{"elastic", e, nu, {"brick.inp", 12U}};
+  return fesa::Hex8Element::Create(7U, layout, x, material, rho,
+                                   {"brick.inp", 4U});
+}
+
+// Independent exact polynomial integrals on the [-h,h] brick. Derivatives
+// remove one linear factor; integral of unlike derivatives is sign-only.
+double GradientProduct(std::size_t a, std::size_t b, std::size_t i,
+                       std::size_t j) {
+  double integral = 48.0 * kSigns[a][i] * kSigns[b][j] /
+                    (64.0 * static_cast<double>((i + 1U) * (j + 1U)));
+  for (std::size_t d = 0; d < 3; ++d)
+    if (d != i && d != j) integral *= 1.0 + kSigns[a][d] * kSigns[b][d] / 3.0;
+  return integral;
+}
+TEST(Hex8Element, ExactBrickStiffnessAndConsistentMass) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  auto stiffness = element.Value().ComputeStiffness();
+  auto mass = element.Value().ComputeMass();
+  ASSERT_TRUE(stiffness.HasValue());
+  ASSERT_TRUE(mass.HasValue());
+  const auto& k = stiffness.Value().values;
+  const auto& m = mass.Value().values;
+  ASSERT_EQ(k.Rows(), 24U);
+  ASSERT_EQ(m.Columns(), 24U);
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t b = 0; b < 8; ++b) {
+      double scalar = 2.0 * 48.0 / 216.0;
+      for (std::size_t d = 0; d < 3; ++d)
+        scalar *= kSigns[a][d] == kSigns[b][d] ? 2.0 : 1.0;
+      for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j) {
+          double exact = 48.0 * GradientProduct(a, b, i, j) +
+                         48.0 * GradientProduct(a, b, j, i);
+          if (i == j)
+            for (std::size_t d = 0; d < 3; ++d)
+              exact += 48.0 * GradientProduct(a, b, d, d);
+          EXPECT_NEAR(k(3 * a + i, 3 * b + j), exact, 2e-12);
+          EXPECT_NEAR(m(3 * a + i, 3 * b + j), i == j ? scalar : 0.0, 2e-14);
+        }
+    }
+}
+TEST(Hex8Element, AffineTensorShearAndEightDistinctRecoveryPoints) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector u{24U};
+  const auto x = Brick();
+  for (std::size_t a = 0; a < 8; ++a) u[3 * a] = 0.02 * x[a][1];
+  auto r = element.Value().EvaluateLinear(u);
+  ASSERT_TRUE(r.HasValue());
+  EXPECT_NEAR(r.Value().strain_energy, 48.0 * 48.0 * 0.0002, 1e-12);
+  for (std::size_t g = 0; g < 8; ++g) {
+    const auto& p = r.Value().points[g];
+    EXPECT_EQ(p.element, 7U);
+    EXPECT_EQ(p.integration_point, static_cast<int>(g + 1));
+    EXPECT_NEAR(p.strain[3], 0.01, 1e-16);
+    EXPECT_NEAR(p.cauchy_stress[3], 0.96, 1e-14);
+    const double gauss = 1.0 / std::sqrt(3.0);
+    const std::array<double, 3> natural{g % 2U == 0U ? -gauss : gauss,
+                                        (g / 2U) % 2U == 0U ? -gauss : gauss,
+                                        g < 4U ? -gauss : gauss};
+    for (std::size_t i = 0; i < 3; ++i) {
+      EXPECT_DOUBLE_EQ(p.natural_coordinates[i], natural[i]);
+      EXPECT_NEAR(p.reference_coordinates[i],
+                  natural[i] * static_cast<double>(i + 1U), 2e-15);
+    }
+  }
+  EXPECT_LT(r.Value().points[2].natural_coordinates[0], 0.0);
+  EXPECT_GT(r.Value().points[2].natural_coordinates[1], 0.0);
+  auto recovered = element.Value().Recover(u);
+  ASSERT_TRUE(recovered.HasValue());
+  EXPECT_EQ(std::get<fesa::SolidElementResultRows>(recovered.Value().payload)
+                .rows.size(),
+            8U);
+}
+TEST(Hex8Element, FiniteRigidRotationHasZeroGreenEnergyAndForce) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector u{24U};
+  const auto x = Brick();
+  for (std::size_t a = 0; a < 8; ++a) {
+    u[3 * a] = -x[a][1] - x[a][0];
+    u[3 * a + 1] = x[a][0] - x[a][1];
+  }
+  auto r = element.Value().EvaluateNonlinear(u);
+  ASSERT_TRUE(r.HasValue());
+  EXPECT_NEAR(r.Value().strain_energy, 0.0, 1e-25);
+  EXPECT_LT(r.Value().internal_force.Norm(), 1e-11);
+}
+TEST(Hex8Element, RejectsCollapsedAndInvalidMaterial) {
+  auto x = Brick();
+  for (auto& a : x) a[2] = 0.0;
+  EXPECT_FALSE(Make(x).HasValue());
+  EXPECT_FALSE(Make(Brick(), 120.0, 0.5).HasValue());
+  EXPECT_FALSE(Make(Brick(), 120.0, -1.0).HasValue());
+}
+
+using Triple = std::array<double, 3>;
+using Tensor = std::array<Triple, 3>;
+Triple Rotate(const Triple& x) {
+  return {-x[1], 0.6 * x[0] - 0.8 * x[2], 0.8 * x[0] + 0.6 * x[2]};
+}
+fesa::Vector Affine(const Coordinates& x, const Tensor& f) {
+  fesa::Vector u{24U};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i) {
+      u[3 * a + i] = -x[a][i];
+      for (std::size_t j = 0; j < 3; ++j) u[3 * a + i] += f[i][j] * x[a][j];
+    }
+  return u;
+}
+Coordinates Warped() {
+  Coordinates x = Brick();
+  for (std::size_t a = 0; a < 8; ++a) {
+    x[a][0] += 0.12 * kSigns[a][1] * kSigns[a][2];
+    x[a][1] += 0.09 * kSigns[a][0] * kSigns[a][2];
+    x[a][2] += 0.15 * kSigns[a][0] * kSigns[a][1];
+  }
+  return x;
+}
+Coordinates QuadraticMap(double a, double b, double c, double d, double e,
+                         double f) {
+  Coordinates x{};
+  for (std::size_t n = 0; n < 8; ++n) {
+    const double xi = kSigns[n][0], eta = kSigns[n][1], zeta = kSigns[n][2];
+    x[n] = {xi, eta * (a + b * xi) + c * xi * zeta,
+            zeta * (d + e * xi) + f * xi * eta};
+  }
+  return x;
+}
+double MaxEntry(const fesa::Matrix& m) {
+  double value = 0.0;
+  for (std::size_t i = 0; i < m.Rows(); ++i)
+    for (std::size_t j = 0; j < m.Columns(); ++j)
+      value = (std::max)(value, std::abs(m(i, j)));
+  return value;
+}
+TEST(Hex8Element, CertifiesPositiveQuadraticAfterSubdivision) {
+  // det J=xi^2-.5xi+.1125, min=.05. Root Bernstein middle=-.8875.
+  auto r = Make(QuadraticMap(.375, -1.0, 1.0, .3, -8.0 / 15.0, -7.0 / 15.0));
+  ASSERT_TRUE(r.HasValue());
+  auto mass = r.Value().ComputeMass();
+  ASSERT_TRUE(mass.HasValue());
+  double total = 0.0;
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t b = 0; b < 8; ++b)
+      total += mass.Value().values(3 * a, 3 * b);
+  EXPECT_NEAR(total, 2.0 * 8.0 * (.1125 + 1.0 / 3.0), 2e-14);
+}
+TEST(Hex8Element, RejectsInteriorFoldDespiteAll27SamplesAndGaussBeingPositive) {
+  // det J=xi^2-.5xi+.0575, min=-.005 at xi=.25. Positive at -1,0,1
+  // and +/-1/sqrt(3); midpoint subdivision detects the hidden folded region.
+  for (double xi :
+       {-1.0, 0.0, 1.0, -1.0 / std::sqrt(3.0), 1.0 / std::sqrt(3.0)})
+    EXPECT_GT(xi * xi - .5 * xi + .0575, 0.0);
+  EXPECT_LT(.25 * .25 - .5 * .25 + .0575, 0.0);
+  const auto folded =
+      QuadraticMap(.23, -1.0, 1.0, .25, -25.0 / 23.0, 2.0 / 23.0);
+  auto invalid = Make(folded);
+  ASSERT_FALSE(invalid.HasValue());
+  EXPECT_NE(invalid.GetStatus().Diagnostics()[0].message.find("nonpositive"),
+            std::string::npos);
+  auto valid = Make();
+  ASSERT_TRUE(valid.HasValue());
+  fesa::Vector u{24U};
+  const auto brick = Brick();
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i)
+      u[3 * a + i] = folded[a][i] - brick[a][i];
+  EXPECT_FALSE(valid.Value().EvaluateNonlinear(u).HasValue());
+  EXPECT_FALSE(valid.Value().EvaluateLinear(u).HasValue());
+}
+TEST(Hex8Element, UnresolvedZeroInteriorJacobianFailsClosed) {
+  // det J=(xi-.25)^2; all27 initial samples are positive but positivity
+  // on the closed cube cannot be certified at the interior collapsed sheet.
+  auto invalid = Make(QuadraticMap(.25, -1.0, 0.0, .25, -1.0, 0.0));
+  ASSERT_FALSE(invalid.HasValue());
+  EXPECT_NE(invalid.GetStatus().Diagnostics()[0].message.find("unresolved"),
+            std::string::npos);
+}
+TEST(Hex8Element, TinyUncertainCurrentAdditionCannotBeAccepted) {
+  Coordinates x = Brick();
+  // The reference is exactly representable with width128. Current exact-real
+  // additions leave width2^-47, invisible when rounded at offset2^58.
+  for (auto& a : x) a[0] = std::ldexp(1.0, 58) + 64.0 * a[0];
+  auto valid = Make(x);
+  ASSERT_TRUE(valid.HasValue());
+  fesa::Vector u{24U};
+  for (std::size_t a = 0; a < 8; ++a) {
+    u[3 * a] = -64.0 * kSigns[a][0];
+    if (kSigns[a][0] > 0) u[3 * a] = std::nextafter(u[3 * a], 0.0);
+  }
+  EXPECT_FALSE(valid.Value().EvaluateNonlinear(u).HasValue());
+}
+TEST(Hex8Element, SixRigidModesAndZeroStateTangent) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  const auto k = element.Value().ComputeStiffness().Value().values;
+  const auto x = Brick();
+  for (std::size_t mode = 0; mode < 6; ++mode) {
+    fesa::Vector u{24U};
+    for (std::size_t a = 0; a < 8; ++a) {
+      if (mode < 3U)
+        u[3 * a + mode] = 1.0;
+      else {
+        const auto axis = mode - 3U;
+        u[3 * a + (axis + 1U) % 3U] = -x[a][(axis + 2U) % 3U];
+        u[3 * a + (axis + 2U) % 3U] = x[a][(axis + 1U) % 3U];
+      }
+    }
+    EXPECT_LT(k.Multiply(u).Norm(), 2e-12);
+  }
+  auto trial = element.Value().EvaluateNonlinear(fesa::Vector{24U});
+  ASSERT_TRUE(trial.HasValue());
+  for (std::size_t i = 0; i < 24; ++i)
+    for (std::size_t j = 0; j < 24; ++j) {
+      EXPECT_NEAR(trial.Value().tangent(i, j), k(i, j), 2e-12);
+      EXPECT_DOUBLE_EQ(trial.Value().geometric_tangent(i, j), 0.0);
+      EXPECT_NEAR(k(i, j), k(j, i), 2e-12);
+    }
+}
+TEST(Hex8Element,
+     ConsistentMassHasAll24PositiveDirectionsAndReferenceInvariance) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  auto before = element.Value().ComputeMass().Value().values;
+  // Tensor-product Walsh vectors diagonalize the independent brick scalar
+  // mass; eigenvalues are rho*V/(8*3^number_of_odd_axes), each repeated3.
+  for (std::size_t mask = 0; mask < 8; ++mask)
+    for (std::size_t i = 0; i < 3; ++i) {
+      fesa::Vector v{24U};
+      unsigned odd = 0U;
+      for (std::size_t d = 0; d < 3; ++d)
+        if ((mask & (1U << d)) != 0U) ++odd;
+      for (std::size_t a = 0; a < 8; ++a) {
+        v[3 * a + i] = 1.0;
+        for (std::size_t d = 0; d < 3; ++d)
+          if ((mask & (1U << d)) != 0U) v[3 * a + i] *= kSigns[a][d];
+      }
+      const double eigenvalue = 12.0 / std::pow(3.0, odd);
+      const auto mv = before.Multiply(v);
+      for (std::size_t a = 0; a < 24; ++a)
+        EXPECT_NEAR(mv[a], eigenvalue * v[a], 3e-14);
+      EXPECT_GT(v.Dot(mv), 0.0);
+    }
+  const Tensor f{{{1.1, .03, 0.0}, {0.0, .97, 0.0}, {0.0, 0.0, 1.02}}};
+  ASSERT_TRUE(element.Value().EvaluateNonlinear(Affine(Brick(), f)).HasValue());
+  auto after = element.Value().ComputeMass().Value().values;
+  for (std::size_t i = 0; i < 24; ++i)
+    for (std::size_t j = 0; j < 24; ++j)
+      EXPECT_DOUBLE_EQ(before(i, j), after(i, j));
+}
+TEST(Hex8Element, GeometryMaterialDensityScalingAndCoordinateCovariance) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  auto x = Brick();
+  auto rotated = x;
+  for (auto& a : rotated) a = Rotate(a);
+  auto rotated_element = Make(rotated);
+  ASSERT_TRUE(rotated_element.HasValue());
+  for (auto& a : x)
+    for (auto& value : a) value *= 2.5;
+  auto scaled = Make(x, 360.0, .25, 10.0);
+  ASSERT_TRUE(scaled.HasValue());
+  const auto k = element.Value().ComputeStiffness().Value().values;
+  const auto m = element.Value().ComputeMass().Value().values;
+  const auto kr = rotated_element.Value().ComputeStiffness().Value().values;
+  const auto ks = scaled.Value().ComputeStiffness().Value().values;
+  const auto ms = scaled.Value().ComputeMass().Value().values;
+  Tensor r{};
+  for (std::size_t j = 0; j < 3; ++j) {
+    Triple basis{};
+    basis[j] = 1.0;
+    const auto column = Rotate(basis);
+    for (std::size_t i = 0; i < 3; ++i) r[i][j] = column[i];
+  }
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t b = 0; b < 8; ++b)
+      for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j) {
+          double oracle = 0.0;
+          for (std::size_t l = 0; l < 3; ++l)
+            for (std::size_t n = 0; n < 3; ++n)
+              oracle += r[i][l] * k(3 * a + l, 3 * b + n) * r[j][n];
+          EXPECT_NEAR(kr(3 * a + i, 3 * b + j), oracle, 3e-12);
+          EXPECT_NEAR(ks(3 * a + i, 3 * b + j), 7.5 * k(3 * a + i, 3 * b + j),
+                      3e-11);
+          EXPECT_NEAR(ms(3 * a + i, 3 * b + j),
+                      5.0 * 15.625 * m(3 * a + i, 3 * b + j), 3e-12);
+        }
+}
+TEST(Hex8Element, WarpedAffinePatchHasExactEnergyAndGlobalCauchyGreenTensors) {
+  const auto x = Warped();
+  auto element = Make(x);
+  ASSERT_TRUE(element.HasValue());
+  const Tensor f{{{1.1, .07, -.04}, {0.0, .96, .02}, {0.0, 0.0, 1.03}}};
+  auto trial = element.Value().EvaluateNonlinear(Affine(x, f));
+  ASSERT_TRUE(trial.HasValue());
+  Tensor green{};
+  Tensor second{};
+  Tensor sigma{};
+  double trace = 0.0, contraction = 0.0;
+  for (std::size_t i = 0; i < 3; ++i)
+    for (std::size_t j = 0; j < 3; ++j) {
+      double c = 0.0;
+      for (std::size_t d = 0; d < 3; ++d) c += f[d][i] * f[d][j];
+      green[i][j] = (c - (i == j ? 1.0 : 0.0)) / 2.0;
+      if (i == j) trace += green[i][j];
+      contraction += green[i][j] * green[i][j];
+    }
+  const double det = 1.1 * .96 * 1.03;
+  for (std::size_t i = 0; i < 3; ++i)
+    for (std::size_t j = 0; j < 3; ++j)
+      second[i][j] = 96.0 * green[i][j] + (i == j ? 48.0 * trace : 0.0);
+  for (std::size_t i = 0; i < 3; ++i)
+    for (std::size_t j = 0; j < 3; ++j)
+      for (std::size_t l = 0; l < 3; ++l)
+        for (std::size_t n = 0; n < 3; ++n)
+          sigma[i][j] += f[i][l] * second[l][n] * f[j][n] / det;
+  // Exact polynomial volume of x=xi+.12eta*zeta,
+  // y=2eta+.09xi*zeta, z=3zeta+.15xi*eta.
+  const double volume =
+      48.0 - (8.0 / 3.0) * (.09 * .15 + 2 * .12 * .15 + 3 * .12 * .09);
+  EXPECT_NEAR(trial.Value().strain_energy,
+              volume * (24.0 * trace * trace + 48.0 * contraction), 2e-12);
+  const std::array<std::array<std::size_t, 2>, 6> component{
+      {{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}};
+  for (const auto& p : trial.Value().points)
+    for (std::size_t c = 0; c < 6; ++c) {
+      const auto i = component[c][0], j = component[c][1];
+      EXPECT_NEAR(p.strain[c], green[i][j], 5e-16);
+      EXPECT_NEAR(p.cauchy_stress[c], sigma[i][j], 8e-14);
+    }
+  // Uniform affine stress has no resultant force or moment on one cell.
+  Triple force{}, moment{};
+  for (std::size_t a = 0; a < 8; ++a) {
+    Triple current{};
+    for (std::size_t i = 0; i < 3; ++i)
+      for (std::size_t j = 0; j < 3; ++j) current[i] += f[i][j] * x[a][j];
+    for (std::size_t i = 0; i < 3; ++i) {
+      force[i] += trial.Value().internal_force[3 * a + i];
+      moment[i] += current[(i + 1U) % 3U] *
+                       trial.Value().internal_force[3 * a + (i + 2U) % 3U] -
+                   current[(i + 2U) % 3U] *
+                       trial.Value().internal_force[3 * a + (i + 1U) % 3U];
+    }
+  }
+  for (std::size_t i = 0; i < 3; ++i) {
+    EXPECT_NEAR(force[i], 0.0, 5e-13);
+    EXPECT_NEAR(moment[i], 0.0, 2e-12);
+  }
+}
+TEST(Hex8Element,
+     StrainedRigidRotationPreservesEnergyAndCovariantForceTangent) {
+  const auto x = Warped();
+  auto element = Make(x);
+  ASSERT_TRUE(element.HasValue());
+  const Tensor f{{{1.1, .07, -.04}, {0.0, .96, .02}, {0.0, 0.0, 1.03}}};
+  const auto u = Affine(x, f);
+  fesa::Vector rotated_u{24U};
+  for (std::size_t a = 0; a < 8; ++a) {
+    Triple current{};
+    for (std::size_t i = 0; i < 3; ++i) current[i] = x[a][i] + u[3 * a + i];
+    const auto rotated = Rotate(current);
+    for (std::size_t i = 0; i < 3; ++i)
+      rotated_u[3 * a + i] = rotated[i] - x[a][i];
+  }
+  const auto baseline = element.Value().EvaluateNonlinear(u);
+  const auto rotated = element.Value().EvaluateNonlinear(rotated_u);
+  ASSERT_TRUE(baseline.HasValue());
+  ASSERT_TRUE(rotated.HasValue());
+  EXPECT_NEAR(rotated.Value().strain_energy, baseline.Value().strain_energy,
+              3e-12);
+  Tensor r{};
+  for (std::size_t j = 0; j < 3; ++j) {
+    Triple basis{};
+    basis[j] = 1.0;
+    const auto column = Rotate(basis);
+    for (std::size_t i = 0; i < 3; ++i) r[i][j] = column[i];
+  }
+  for (std::size_t a = 0; a < 8; ++a) {
+    const Triple original{baseline.Value().internal_force[3 * a],
+                          baseline.Value().internal_force[3 * a + 1],
+                          baseline.Value().internal_force[3 * a + 2]};
+    const auto expected = Rotate(original);
+    for (std::size_t i = 0; i < 3; ++i)
+      EXPECT_NEAR(rotated.Value().internal_force[3 * a + i], expected[i],
+                  3e-12);
+    for (std::size_t b = 0; b < 8; ++b)
+      for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t j = 0; j < 3; ++j) {
+          double expected_k = 0.0;
+          for (std::size_t l = 0; l < 3; ++l)
+            for (std::size_t n = 0; n < 3; ++n)
+              expected_k += r[i][l] *
+                            baseline.Value().tangent(3 * a + l, 3 * b + n) *
+                            r[j][n];
+          EXPECT_NEAR(rotated.Value().tangent(3 * a + i, 3 * b + j), expected_k,
+                      7e-12);
+        }
+  }
+}
+TEST(Hex8Element, EnergyForceAndForceTangentDirectionalStepSweep) {
+  auto element = Make(Warped());
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector u{24U}, direction{24U};
+  for (std::size_t i = 0; i < 24; ++i) {
+    u[i] = .05 * std::sin(.7 * static_cast<double>(i));
+    direction[i] = std::cos(1.1 * static_cast<double>(i));
+  }
+  const auto trial = element.Value().EvaluateNonlinear(u);
+  ASSERT_TRUE(trial.HasValue());
+  const double derivative = direction.Dot(trial.Value().internal_force);
+  const auto tangent_direction = trial.Value().tangent.Multiply(direction);
+  std::array<double, 6> energy_error{}, force_error{};
+  const std::array<double, 6> steps{1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7};
+  for (std::size_t s = 0; s < steps.size(); ++s) {
+    auto plus = u, minus = u;
+    plus.Axpy(steps[s], direction);
+    minus.Axpy(-steps[s], direction);
+    const auto rp = element.Value().EvaluateNonlinear(plus),
+               rm = element.Value().EvaluateNonlinear(minus);
+    ASSERT_TRUE(rp.HasValue());
+    ASSERT_TRUE(rm.HasValue());
+    energy_error[s] =
+        std::abs((rp.Value().strain_energy - rm.Value().strain_energy) /
+                     (2 * steps[s]) -
+                 derivative) /
+        std::abs(derivative);
+    auto difference = rp.Value().internal_force;
+    difference.Axpy(-1.0, rm.Value().internal_force);
+    difference.Scale(1.0 / (2 * steps[s]));
+    difference.Axpy(-1.0, tangent_direction);
+    force_error[s] = difference.Norm() / tangent_direction.Norm();
+    std::cout << "h=" << steps[s] << " energy->f=" << energy_error[s]
+              << " f->Kt=" << force_error[s] << '\n';
+  }
+  EXPECT_LT(energy_error[1], .02 * energy_error[0]);
+  EXPECT_LT(force_error[1], .02 * force_error[0]);
+  EXPECT_LT(*std::min_element(energy_error.begin(), energy_error.end()), 1e-7);
+  EXPECT_LT(*std::min_element(force_error.begin(), force_error.end()), 1e-7);
+  EXPECT_GT(energy_error.back(),
+            *std::min_element(energy_error.begin(), energy_error.end()));
+  for (std::size_t i = 0; i < 24; ++i)
+    for (std::size_t j = 0; j < 24; ++j) {
+      EXPECT_NEAR(trial.Value().tangent(i, j), trial.Value().tangent(j, i),
+                  MaxEntry(trial.Value().tangent) * 2e-14);
+      EXPECT_DOUBLE_EQ(trial.Value().tangent(i, j),
+                       trial.Value().material_tangent(i, j) +
+                           trial.Value().geometric_tangent(i, j));
+    }
+}
+TEST(Hex8Element,
+     InvalidInputsAndFailedTrialPreservePriorResponseAndReference) {
+  EXPECT_FALSE(Make(Brick(), 0.0).HasValue());
+  EXPECT_FALSE(
+      Make(Brick(), std::numeric_limits<double>::infinity()).HasValue());
+  EXPECT_FALSE(Make(Brick(), 120.0, std::numeric_limits<double>::quiet_NaN())
+                   .HasValue());
+  EXPECT_FALSE(Make(Brick(), 120.0, .25, 0.0).HasValue());
+  EXPECT_FALSE(Make(Brick(), 120.0, .25, -1.0).HasValue());
+  EXPECT_FALSE(
+      Make(Brick(), 120.0, .25, std::numeric_limits<double>::infinity())
+          .HasValue());
+  auto nonfinite = Brick();
+  nonfinite[2][1] = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_FALSE(Make(nonfinite).HasValue());
+  auto no_density = Make(Brick(), 120.0, .25, std::nullopt);
+  ASSERT_TRUE(no_density.HasValue());
+  EXPECT_FALSE(no_density.Value().ComputeMass().HasValue());
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  const auto u = Affine(
+      Brick(), Tensor{{{1.01, .01, 0.0}, {0.0, 1.02, 0.0}, {0.0, 0.0, .98}}});
+  const auto accepted = element.Value().EvaluateNonlinear(u);
+  ASSERT_TRUE(accepted.HasValue());
+  auto failed = u;
+  failed[4] = std::numeric_limits<double>::infinity();
+  EXPECT_FALSE(element.Value().EvaluateNonlinear(failed).HasValue());
+  EXPECT_FALSE(element.Value().EvaluateLinear(fesa::Vector{23U}).HasValue());
+  auto inverted = Brick();
+  for (auto& a : inverted) a[0] = -a[0];
+  EXPECT_FALSE(Make(inverted).HasValue());
+  const auto again = element.Value().EvaluateNonlinear(u);
+  ASSERT_TRUE(again.HasValue());
+  EXPECT_DOUBLE_EQ(accepted.Value().strain_energy, again.Value().strain_energy);
+  for (std::size_t i = 0; i < 24; ++i) {
+    EXPECT_DOUBLE_EQ(accepted.Value().internal_force[i],
+                     again.Value().internal_force[i]);
+    for (std::size_t j = 0; j < 24; ++j)
+      EXPECT_DOUBLE_EQ(accepted.Value().tangent(i, j),
+                       again.Value().tangent(i, j));
+  }
+  const auto failure = element.Value().EvaluateNonlinear(failed);
+  ASSERT_FALSE(failure.HasValue());
+  EXPECT_EQ(failure.GetStatus().Category(), fesa::FailureCategory::kModel);
+  const auto& diagnostic = failure.GetStatus().Diagnostics()[0];
+  EXPECT_EQ(diagnostic.location.line, 4U);
+  EXPECT_EQ(diagnostic.keyword, "C3D8");
+  EXPECT_EQ(diagnostic.entity_identity, "PART-1:42");
+}
+
+TEST(Hex8Element, InvalidLayoutsAreRejectedBeforeNumericalWork) {
+  fesa::LinearElasticMaterial material{
+      "elastic", 120.0, .25, {"brick.inp", 12U}};
+  fesa::ElementDofLayout layout{
+      {"PART-1", 42, "42"},
+      {0, 1, 2, 3, 4, 5, 6, 7},
+      {fesa::DofComponent::kUx, fesa::DofComponent::kUy,
+       fesa::DofComponent::kUz}};
+  auto duplicate = layout;
+  duplicate.node_indices[7] = 0U;
+  EXPECT_FALSE(
+      fesa::Hex8Element::Create(7U, duplicate, Brick(), material).HasValue());
+  auto missing = layout;
+  missing.node_indices.pop_back();
+  EXPECT_FALSE(
+      fesa::Hex8Element::Create(7U, missing, Brick(), material).HasValue());
+  auto wrong_order = layout;
+  std::swap(wrong_order.components_per_node[0],
+            wrong_order.components_per_node[1]);
+  EXPECT_FALSE(
+      fesa::Hex8Element::Create(7U, wrong_order, Brick(), material).HasValue());
+  auto rotation = layout;
+  rotation.components_per_node.push_back(fesa::DofComponent::kUrx);
+  EXPECT_FALSE(
+      fesa::Hex8Element::Create(7U, rotation, Brick(), material).HasValue());
+}
+TEST(Hex8Element, LinearStiffnessHasOnlyTheSixRigidNullDirections) {
+  auto element = Make();
+  ASSERT_TRUE(element.HasValue());
+  const auto k = element.Value().ComputeStiffness().Value().values;
+  // Independent minimal rigid-body fixing: node0 xyz, node1 yz, node3 z.
+  // Cholesky positivity of the principal18x18 proves no extra null mode;
+  // full K nullity is at most6, combined with six independent rigid modes.
+  std::array<std::size_t, 18> free{};
+  std::size_t count = 0U;
+  for (std::size_t i = 0; i < 24; ++i)
+    if (i != 0U && i != 1U && i != 2U && i != 4U && i != 5U && i != 11U)
+      free[count++] = i;
+  std::array<std::array<double, 18>, 18> lower{};
+  for (std::size_t i = 0; i < 18; ++i)
+    for (std::size_t j = 0; j <= i; ++j) {
+      double value = k(free[i], free[j]);
+      for (std::size_t n = 0; n < j; ++n) value -= lower[i][n] * lower[j][n];
+      if (i == j) {
+        ASSERT_GT(value, 0.0);
+        lower[i][j] = std::sqrt(value);
+      } else
+        lower[i][j] = value / lower[j][j];
+    }
+}
+TEST(Hex8Element, NonlinearSmallLoadApproachesLinearResponse) {
+  auto element = Make(Warped());
+  ASSERT_TRUE(element.HasValue());
+  fesa::Vector direction{24U};
+  for (std::size_t i = 0; i < 24; ++i)
+    direction[i] = std::sin(.7 * static_cast<double>(i));
+  const auto k = element.Value().ComputeStiffness().Value().values;
+  std::array<double, 3> error{};
+  const std::array<double, 3> factors{1e-2, 1e-3, 1e-4};
+  for (std::size_t n = 0; n < factors.size(); ++n) {
+    auto u = direction;
+    u.Scale(factors[n]);
+    const auto nonlinear = element.Value().EvaluateNonlinear(u);
+    const auto linear = element.Value().EvaluateLinear(u);
+    ASSERT_TRUE(nonlinear.HasValue());
+    ASSERT_TRUE(linear.HasValue());
+    auto difference = nonlinear.Value().internal_force;
+    difference.Axpy(-1.0, linear.Value().internal_force);
+    error[n] = difference.Norm() / linear.Value().internal_force.Norm();
+    EXPECT_NEAR(linear.Value().strain_energy, .5 * u.Dot(k.Multiply(u)), 1e-14);
+  }
+  EXPECT_LT(error[1], .12 * error[0]);
+  EXPECT_LT(error[2], .12 * error[1]);
+}
+TEST(Hex8Element, RejectsUnrepresentableShearAndMassInsteadOfReturningZeros) {
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  EXPECT_FALSE(Make(Brick(), tiny, .25).HasValue());
+  auto element = Make(Brick(), 120.0, .25, tiny);
+  ASSERT_TRUE(element.HasValue());
+  EXPECT_FALSE(element.Value().ComputeMass().HasValue());
+}
+}  // namespace