#include "fesa/elements/hex8_element.h"

#include <gtest/gtest.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <iostream>
#include <limits>
#include <optional>
#include <string>
#include <vector>

#include "../../../src/fesa/math/high_low_arithmetic.h"

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 = -(2.0 / 3.0) * 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);
          // Mean bulk term is the product of two independent first moments.
          exact += (80.0 / 48.0) * (48.0 * kSigns[a][i] / (8.0 * (i + 1))) *
                   (48.0 * kSigns[b][j] / (8.0 * (j + 1)));
          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, UniformTranslationHasExactlyZeroLinearForceAndEnergy) {
  auto element = Make(Brick(), 2.1e11, 0.3);
  ASSERT_TRUE(element.HasValue());
  fesa::Vector u{24U};
  for (std::size_t a = 0; a < 8; ++a) {
    u[3U * a] = 0.125;
    u[3U * a + 1U] = -0.25;
    u[3U * a + 2U] = 0.5;
  }
  const auto response = element.Value().EvaluateLinear(u);
  ASSERT_TRUE(response.HasValue());
  EXPECT_EQ(response.Value().strain_energy, 0.0);
  for (std::size_t i = 0; i < 24; ++i)
    EXPECT_EQ(response.Value().internal_force[i], 0.0);
  for (const auto& point : response.Value().points)
    for (std::size_t i = 0; i < 6; ++i) {
      EXPECT_EQ(point.strain[i], 0.0);
      EXPECT_EQ(point.cauchy_stress[i], 0.0);
    }
}
TEST(Hex8Element, LinearResponseIsUnchangedByExactlyRepresentableTranslation) {
  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[3U * a] = x[a][0] / 64.0 + x[a][1] / 128.0;
    u[3U * a + 1U] = x[a][2] / 256.0;
  }
  const auto original = element.Value().EvaluateLinear(u);
  ASSERT_TRUE(original.HasValue());
  for (std::size_t a = 0; a < 8; ++a) {
    u[3U * a] += 1024.0;
    u[3U * a + 1U] -= 512.0;
    u[3U * a + 2U] += 256.0;
  }
  const auto translated = element.Value().EvaluateLinear(u);
  ASSERT_TRUE(translated.HasValue());
  EXPECT_EQ(translated.Value().strain_energy, original.Value().strain_energy);
  for (std::size_t i = 0; i < 24; ++i)
    EXPECT_EQ(translated.Value().internal_force[i],
              original.Value().internal_force[i]);
  for (std::size_t g = 0; g < 8; ++g) {
    EXPECT_EQ(translated.Value().points[g].strain,
              original.Value().points[g].strain);
    EXPECT_EQ(translated.Value().points[g].cauchy_stress,
              original.Value().points[g].cauchy_stress);
  }
}
TEST(Hex8Element, TinySignedAffineStretchRetainsStvkStressForceAndEnergy) {
  constexpr double youngs_modulus = 2.1e11;
  constexpr double poisson_ratio = 0.3;
  const double lambda = youngs_modulus * poisson_ratio /
                        ((1.0 + poisson_ratio) * (1.0 - 2.0 * poisson_ratio));
  const double mu = youngs_modulus / (2.0 * (1.0 + poisson_ratio));
  const double axial_modulus = lambda + 2.0 * mu;
  auto element = Make(Brick(), youngs_modulus, poisson_ratio);
  ASSERT_TRUE(element.HasValue());
  const auto x = Brick();
  for (const double stretch : {1e-10, -1e-10, 1e-8, -1e-8}) {
    SCOPED_TRACE(stretch);
    fesa::Vector u{24U};
    for (std::size_t a = 0; a < 8; ++a)
      u[3U * a] = stretch * x[a][0];
    const auto response = element.Value().EvaluateNonlinear(u);
    ASSERT_TRUE(response.HasValue());
    // Independent homogeneous diagonal F and reference face-traction integrals.
    const double green = stretch + 0.5 * stretch * stretch;
    const double sxx = axial_modulus * green;
    const double syy = lambda * green;
    const double expected_energy = 24.0 * axial_modulus * green * green;
    EXPECT_NEAR(response.Value().strain_energy, expected_energy,
                1e-12 * expected_energy);
    for (const auto& point : response.Value().points) {
      EXPECT_NEAR(point.strain[0], green, 1e-12 * std::abs(green));
      EXPECT_NEAR(point.cauchy_stress[0], (1.0 + stretch) * sxx,
                  1e-12 * std::abs(sxx));
      EXPECT_NEAR(point.cauchy_stress[1], syy / (1.0 + stretch),
                  1e-12 * std::abs(syy));
      EXPECT_NEAR(point.cauchy_stress[2], syy / (1.0 + stretch),
                  1e-12 * std::abs(syy));
    }
    for (std::size_t a = 0; a < 8; ++a) {
      const std::array<double, 3> expected{
          6.0 * kSigns[a][0] * (1.0 + stretch) * sxx,
          3.0 * kSigns[a][1] * syy, 2.0 * kSigns[a][2] * syy};
      for (std::size_t i = 0; i < 3; ++i)
        EXPECT_NEAR(response.Value().internal_force[3U * a + i], expected[i],
                    1e-12 * std::abs(expected[i]));
    }
  }
}
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 mr = rotated_element.Value().ComputeMass().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);
          // Isotropic translation mass is unchanged by a rigid coordinate
          // rotation; this equals blockwise R*M*R^T in global components.
          EXPECT_NEAR(mr(3 * a + i, 3 * b + j), m(3 * a + i, 3 * b + j), 3e-14);
          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, 8> energy_error{}, force_error{};
  const std::array<double, 8> steps{1e-2, 1e-3, 1e-4, 1e-5,
                                    1e-6, 1e-7, 1e-8, 1e-9};
  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, TranslatedWarpedLinearResponseRetainsEnergyAndTangentVariation) {
  auto element = Make(Warped());
  ASSERT_TRUE(element.HasValue());
  fesa::Vector u{24U}, direction{24U};
  const std::array<double, 3> translation{1024.0, -512.0, 256.0};
  for (std::size_t i = 0; i < 24; ++i) {
    u[i] = translation[i % 3U] + .05 * std::sin(.7 * static_cast<double>(i));
    direction[i] = std::cos(1.1 * static_cast<double>(i));
  }
  const auto trial = element.Value().EvaluateLinear(u);
  ASSERT_TRUE(trial.HasValue());
  const double derivative = direction.Dot(trial.Value().internal_force);
  const auto tangent_direction = trial.Value().tangent.Multiply(direction);
  ASSERT_NE(derivative, 0.0);
  ASSERT_GT(tangent_direction.Norm(), 0.0);
  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().EvaluateLinear(plus);
    const auto rm = element.Value().EvaluateLinear(minus);
    ASSERT_TRUE(rp.HasValue());
    ASSERT_TRUE(rm.HasValue());
    energy_error[s] =
        std::abs((rp.Value().strain_energy - rm.Value().strain_energy) /
                     (2.0 * 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.0 * steps[s]));
    difference.Axpy(-1.0, tangent_direction);
    force_error[s] = difference.Norm() / tangent_direction.Norm();
    std::cout << "translated-linear h=" << steps[s]
              << " energy->f=" << energy_error[s]
              << " f->K0=" << force_error[s] << '\n';
  }
  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);
}
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());
}
TEST(Hex8Element, DirectBulkSurvivesNearNegativePoissonBound) {
  const double nu = std::nextafter(-1.0, 0.0);
  auto element = Make(Brick(), 1.0, nu);
  ASSERT_TRUE(element.HasValue());
  // Direct kappa is positive even where a lambda-based difference loses it.
  EXPECT_GT(1.0 / (3.0 * (1.0 - 2.0 * nu)), 0.0);
  fesa::Vector zero(24);
  EXPECT_TRUE(element.Value().EvaluateLinear(zero).HasValue());
}
TEST(Hex8Element, OrdinaryNegativePoissonRatioRetainsPositiveVolumetricEnergy) {
  auto element = Make(Brick(), 120.0, -.5);
  ASSERT_TRUE(element.HasValue());
  const auto u = Affine(
      Brick(), Tensor{{{1.01, 0.0, 0.0}, {0.0, 1.01, 0.0}, {0.0, 0.0, 1.01}}});
  const auto response = element.Value().EvaluateLinear(u);
  ASSERT_TRUE(response.HasValue());
  const double bulk = 120.0 / (3.0 * (1.0 - 2.0 * (-.5)));
  const double trace = .03;
  EXPECT_NEAR(response.Value().strain_energy, 48.0 * .5 * bulk * trace * trace,
              2e-14);
  for (const auto& point : response.Value().points)
    for (std::size_t i = 0; i < 3; ++i)
      EXPECT_NEAR(point.cauchy_stress[i], bulk * trace, 2e-14);
}

TEST(Hex8Element, FiniteInputsThatOverflowDerivedArithmeticFailClosed) {
  const double huge = std::numeric_limits<double>::max();
  auto bad_material = Make(Brick(), huge, .49);
  ASSERT_FALSE(bad_material.HasValue());
  EXPECT_NE(bad_material.GetStatus().Diagnostics()[0].message.find("finite"),
            std::string::npos);
  auto x = Brick();
  for (auto& a : x)
    for (auto& value : a) value *= 1e200;
  auto bad_geometry = Make(x);
  ASSERT_FALSE(bad_geometry.HasValue());
  EXPECT_NE(bad_geometry.GetStatus().Diagnostics()[0].message.find("nonfinite"),
            std::string::npos);
  auto element = Make(Brick(), 120.0, .25, huge);
  ASSERT_TRUE(element.HasValue());
  auto bad_mass = element.Value().ComputeMass();
  ASSERT_FALSE(bad_mass.HasValue());
  EXPECT_NE(bad_mass.GetStatus().Diagnostics()[0].message.find("nonfinite"),
            std::string::npos);
  // Current detJ is finite (~6e300), but the StVK energy and pushed-forward
  // stress arithmetic overflow despite every input displacement being finite.
  const auto u =
      Affine(Brick(),
             Tensor{{{1e100, 0.0, 0.0}, {0.0, 1e100, 0.0}, {0.0, 0.0, 1e100}}});
  auto bad_trial = element.Value().EvaluateNonlinear(u);
  ASSERT_FALSE(bad_trial.HasValue());
  EXPECT_NE(bad_trial.GetStatus().Diagnostics()[0].message.find("trial"),
            std::string::npos);
  EXPECT_NE(bad_trial.GetStatus().Diagnostics()[0].message.find("nonfinite"),
            std::string::npos);
}

// Independent Gauss-Legendre abscissae from polynomial root iteration. This
// test does not call the kernel's shapes, Jacobian, inverse or quadrature.
std::vector<std::array<double, 2>> LegendreRule(unsigned count) {
  std::vector<std::array<double, 2>> points;
  for (unsigned root = 0; root < count; ++root) {
    double x = std::cos(std::acos(-1.0) * (static_cast<double>(root) + .75) /
                        (count + .5));
    double derivative = 0.0;
    for (unsigned iteration = 0; iteration < 30; ++iteration) {
      double previous = 1.0, current = x;
      for (unsigned n = 2; n <= count; ++n) {
        const double next =
            ((2.0 * n - 1.0) * x * current - (n - 1.0) * previous) / n;
        previous = current;
        current = next;
      }
      derivative = count * (x * current - previous) / (x * x - 1.0);
      const double correction = current / derivative;
      x -= correction;
      if (std::abs(correction) < 2e-16) break;
    }
    points.push_back({x, 2.0 / ((1.0 - x * x) * derivative * derivative)});
  }
  return points;
}
Triple Cross(const Triple& a, const Triple& b) {
  return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
          a[0] * b[1] - a[1] * b[0]};
}
std::array<double, 2> WarpedQuadratureEnergies(const fesa::Vector& u,
                                               unsigned count,
                                               bool nonlinear = false) {
  const auto rule = LegendreRule(count);
  std::array<double, 2> result{};
  double volume = 0.0, trace_integral = 0.0;
  for (const auto& xi : rule)
    for (const auto& eta : rule)
      for (const auto& zeta : rule) {
        const Triple natural{xi[0], eta[0], zeta[0]};
        // Derivatives of the analytic warped map, rather than recovered nodal
        // J.
        const Triple c0{1.0, .09 * zeta[0], .15 * eta[0]};
        const Triple c1{.12 * zeta[0], 2.0, .15 * xi[0]};
        const Triple c2{.12 * eta[0], .09 * xi[0], 3.0};
        Tensor reciprocal{Cross(c1, c2), Cross(c2, c0), Cross(c0, c1)};
        const double det = c0[0] * reciprocal[0][0] + c0[1] * reciprocal[0][1] +
                           c0[2] * reciprocal[0][2];
        for (auto& row : reciprocal)
          for (auto& value : row) value /= det;
        Tensor natural_gradient{};
        Triple velocity{};
        for (std::size_t a = 0; a < 8; ++a) {
          double shape = .125;
          for (std::size_t j = 0; j < 3; ++j)
            shape *= 1.0 + kSigns[a][j] * natural[j];
          for (std::size_t i = 0; i < 3; ++i)
            velocity[i] += shape * u[3 * a + i];
          for (std::size_t j = 0; j < 3; ++j) {
            double derivative = .125 * kSigns[a][j];
            for (std::size_t other = 0; other < 3; ++other)
              if (other != j)
                derivative *= 1.0 + kSigns[a][other] * natural[other];
            for (std::size_t i = 0; i < 3; ++i)
              natural_gradient[i][j] += derivative * u[3 * a + i];
          }
        }
        Tensor gradient{};
        for (std::size_t i = 0; i < 3; ++i)
          for (std::size_t j = 0; j < 3; ++j)
            for (std::size_t d = 0; d < 3; ++d)
              gradient[i][j] += natural_gradient[i][d] * reciprocal[d][j];
        Tensor strains{};
        for (std::size_t i = 0; i < 3; ++i)
          for (std::size_t j = 0; j < 3; ++j) {
            strains[i][j] = .5 * (gradient[i][j] + gradient[j][i]);
            if (nonlinear)
              for (std::size_t k = 0; k < 3; ++k)
                strains[i][j] += .5 * gradient[k][i] * gradient[k][j];
          }
        double trace = 0.0, strain_squared = 0.0, velocity_squared = 0.0;
        for (std::size_t i = 0; i < 3; ++i) {
          trace += strains[i][i];
          velocity_squared += velocity[i] * velocity[i];
          for (std::size_t j = 0; j < 3; ++j) {
            const double strain = strains[i][j];
            strain_squared += strain * strain;
          }
        }
        const double weight = xi[1] * eta[1] * zeta[1] * det;
        result[0] += weight * 48.0 * (strain_squared - trace * trace / 3.0);
        volume += weight;
        trace_integral += weight * trace;
        result[1] += weight * velocity_squared;  // .5*rho*v^2 with rho=2.
      }
  result[0] += 40.0 * trace_integral * trace_integral / volume;
  return result;
}
TEST(Hex8Element, WarpedNonaffineQuadratureDifferenceIsQuantifiedAndRefines) {
  auto element = Make(Warped());
  ASSERT_TRUE(element.HasValue());
  fesa::Vector u{24U};
  for (std::size_t i = 0; i < 24; ++i)
    u[i] = .03 * std::sin(.7 * static_cast<double>(i));
  const auto k = element.Value().ComputeStiffness().Value().values;
  const auto m = element.Value().ComputeMass().Value().values;
  const auto q2 = WarpedQuadratureEnergies(u, 2),
             q3 = WarpedQuadratureEnergies(u, 3);
  const auto q5 = WarpedQuadratureEnergies(u, 5),
             q7 = WarpedQuadratureEnergies(u, 7);
  const double budget = 512.0 * std::numeric_limits<double>::epsilon();
  EXPECT_NEAR(.5 * u.Dot(k.Multiply(u)), q2[0], budget * q2[0]);
  EXPECT_NEAR(.5 * u.Dot(m.Multiply(u)), q2[1], budget * q2[1]);
  const double stiffness_discrepancy = std::abs(q2[0] - q7[0]) / q7[0];
  const double mass_discrepancy = std::abs(q2[1] - q7[1]) / q7[1];
  // Mandated2x2 is retained: higher-order integration gives a distinct warped
  // discretization. Rational stiffness integrands converge with refinement.
  EXPECT_GT(stiffness_discrepancy, 1e-8);
  EXPECT_GT(mass_discrepancy, 1e-8);
  EXPECT_LT(std::abs(q5[0] - q7[0]), .05 * std::abs(q3[0] - q7[0]));
  // N_a*N_b*detJ has per-axis degree<=4, so order3 integrates mass exactly.
  EXPECT_NEAR(q3[1], q7[1], budget * q7[1]);
  EXPECT_NEAR(q5[1], q7[1], budget * q7[1]);
  std::cout << "warped K2-vs-K7 energy relative=" << stiffness_discrepancy
            << " M2-vs-M7 kinetic relative=" << mass_discrepancy
            << " K3-vs-K7=" << std::abs(q3[0] - q7[0]) / q7[0]
            << " K5-vs-K7=" << std::abs(q5[0] - q7[0]) / q7[0] << '\n';
}
Coordinates MeanVolumeUnitBrick() {
  Coordinates x{};
  for (std::size_t a = 0; a < 8; ++a)
    for (std::size_t i = 0; i < 3; ++i) x[a][i] = .5 * (kSigns[a][i] + 1);
  return x;
}

TEST(Hex8Element,
     MeanVolumeNonaffinePolynomialEnergyForceHessianAndActualStrain) {
  const auto x = MeanVolumeUnitBrick();
  const auto element = Make(x, 12.0, .25);
  ASSERT_TRUE(element.HasValue());
  constexpr double mu = 4.8, bulk = 8.0, c = .125;
  // Exact monomial integrals of u_x=c*x*y over the unit cube. The other
  // translations vanish; this direction varies trace among all eight points.
  constexpr double a = 7 * mu / 18 + bulk / 8;
  constexpr double b = 5 * mu / 18 + bulk / 6;
  // Fixed two-point quadrature has Q[x^4]=7/36, not the continuum 1/5.
  constexpr double d = 11 * mu / 108 + bulk / 18;
  fesa::Vector u(24);
  for (std::size_t n = 0; n < 8; ++n) u[3 * n] = c * x[n][0] * x[n][1];
  for (bool nonlinear : {false, true}) {
    const auto response = nonlinear ? element.Value().EvaluateNonlinear(u)
                                    : element.Value().EvaluateLinear(u);
    ASSERT_TRUE(response.HasValue());
    const auto& r = response.Value();
    const double energy =
        a * c * c + (nonlinear ? b * c * c * c + d * c * c * c * c : 0);
    const double force =
        2 * a * c + (nonlinear ? 3 * b * c * c + 4 * d * c * c * c : 0);
    const double hessian = 2 * a + (nonlinear ? 6 * b * c + 12 * d * c * c : 0);
    EXPECT_NEAR(r.strain_energy, energy, 2e-14);
    EXPECT_NEAR(r.internal_force[6] + r.internal_force[18], force, 2e-13);
    EXPECT_NEAR(r.tangent(6, 6) + r.tangent(6, 18) + r.tangent(18, 6) +
                    r.tangent(18, 18),
                hessian, 2e-12);
    const double mean = c / 2 + (nonlinear ? c * c / 3 : 0);
    for (const auto& point : r.points) {
      const double px = point.reference_coordinates[0],
                   py = point.reference_coordinates[1];
      const double t = 1 + c * py;
      const double exx = c * py + (nonlinear ? .5 * c * c * py * py : 0);
      const double eyy = nonlinear ? .5 * c * c * px * px : 0;
      const double exy = .5 * c * px + (nonlinear ? .5 * c * c * px * py : 0);
      const double trace = exx + eyy;
      const double sxx = 2 * mu * (exx - trace / 3) + bulk * mean;
      const double syy = 2 * mu * (eyy - trace / 3) + bulk * mean;
      const double szz = -2 * mu * trace / 3 + bulk * mean;
      const double sxy = 2 * mu * exy;
      EXPECT_NEAR(point.strain[0], exx, 2e-15);
      EXPECT_NEAR(point.strain[1], eyy, 2e-15);
      EXPECT_NEAR(point.strain[3], exy, 2e-15);
      EXPECT_NEAR(
          point.cauchy_stress[0],
          nonlinear
              ? (t * t * sxx + 2 * t * c * px * sxy + c * c * px * px * syy) / t
              : sxx,
          2e-13);
      EXPECT_NEAR(point.cauchy_stress[1], nonlinear ? syy / t : syy, 2e-13);
      EXPECT_NEAR(point.cauchy_stress[2], nonlinear ? szz / t : szz, 2e-13);
      EXPECT_NEAR(point.cauchy_stress[3],
                  nonlinear ? (t * sxy + c * px * syy) / t : sxy, 2e-13);
    }
  }
  const auto k = element.Value().ComputeStiffness().Value().values;
  EXPECT_NEAR(k(0, 0), 41.0 / 18, 2e-14);
  EXPECT_NEAR(k(0, 1), 19.0 / 30, 2e-14);
}

TEST(Hex8Element, MeanVolumeWarpedStressUsesReferenceWeightedTrace) {
  auto element = Make(Warped());
  ASSERT_TRUE(element.HasValue());
  fesa::Vector u(24);
  u[1] = 1;  // Independent derivative coefficients separate the two means by
             // >1e-8.
  const auto response = element.Value().EvaluateLinear(u);
  ASSERT_TRUE(response.HasValue());
  const auto weights = element.Value().ReferenceWeights();
  double volume = 0, integral = 0, arithmetic = 0;
  for (std::size_t g = 0; g < 8; ++g) {
    const auto& e = response.Value().points[g].strain;
    const double trace = e[0] + e[1] + e[2];
    volume += weights[g];
    integral += weights[g] * trace;
    arithmetic += trace / 8;
  }
  const double mean = integral / volume;
  ASSERT_GT(std::abs(mean - arithmetic), 1e-8);
  for (const auto& point : response.Value().points) {
    const auto& s = point.cauchy_stress;
    EXPECT_NEAR((s[0] + s[1] + s[2]) / 3, 80 * mean, 2e-13);
  }
  EXPECT_NEAR(response.Value().strain_energy, WarpedQuadratureEnergies(u, 2)[0],
              2e-12);
}

TEST(Hex8Element,
     MeanVolumeDirectBulkSupportsBoundsWithoutLambdaOrThreeBulkGuard) {
  EXPECT_TRUE(
      Make(MeanVolumeUnitBrick(), 12, std::nextafter(-1.0, 0.0)).HasValue());
  EXPECT_TRUE(
      Make(MeanVolumeUnitBrick(), 1, std::nextafter(.5, 0.0)).HasValue());
  const auto huge = Make(MeanVolumeUnitBrick(), 1e308, .4);
  ASSERT_TRUE(huge.HasValue());
  const double direct_bulk = 1e308 / (3 * (1 - 2 * .4));
  ASSERT_TRUE(std::isfinite(direct_bulk));
  ASSERT_FALSE(std::isfinite(3 * direct_bulk));
  fesa::Vector u(24);
  constexpr double shear = 1e-150;
  const auto x = MeanVolumeUnitBrick();
  for (unsigned n = 0; n < 8; ++n) u[3 * n] = shear * x[n][1];
  const auto response = huge.Value().EvaluateLinear(u);
  ASSERT_TRUE(response.HasValue());
  const double expected = .5 * (1e308 / (2 * (1 + .4))) * shear * shear;
  EXPECT_NEAR(response.Value().strain_energy, expected, 1e-13 * expected);
  EXPECT_FALSE(Make(MeanVolumeUnitBrick(), 1e308, .499).HasValue());
  EXPECT_FALSE(Make(MeanVolumeUnitBrick(),
                    std::numeric_limits<double>::denorm_min(), -.5)
                   .HasValue());
}
TEST(Hex8Element,
     MeanVolumeWarpedIndependentEnergyHasCoupledNodeZeroDerivatives) {
  auto element = Make(Warped());
  ASSERT_TRUE(element.HasValue());
  fesa::Vector u(24), direction(24), zero(24);
  for (std::size_t i = 0; i < 24; ++i) u[i] = .03 * std::sin(.7 * i);
  direction[0] = 1;
  direction[8] = -.375;
  direction[22] = .25;
  constexpr double h = .01;
  for (bool nonlinear : {false, true}) {
    const auto response = nonlinear ? element.Value().EvaluateNonlinear(u)
                                    : element.Value().EvaluateLinear(u);
    const auto paired = nonlinear
                            ? element.Value().EvaluateNonlinearHighLow(u, zero)
                            : element.Value().EvaluateLinearHighLow(u, zero);
    ASSERT_TRUE(response.HasValue());
    ASSERT_TRUE(paired.HasValue());
    const double energy = WarpedQuadratureEnergies(u, 2, nonlinear)[0];
    EXPECT_NEAR(response.Value().strain_energy, energy, 2e-12);
    EXPECT_NEAR(paired.Value().high.strain_energy, energy, 2e-12);
    std::array<double, 2> first{}, second{};
    for (unsigned level = 0; level < 2; ++level) {
      const double step = h / (1U << level);
      auto plus = u, minus = u;
      plus.Axpy(step, direction);
      minus.Axpy(-step, direction);
      const double ep = WarpedQuadratureEnergies(plus, 2, nonlinear)[0];
      const double em = WarpedQuadratureEnergies(minus, 2, nonlinear)[0];
      first[level] = (ep - em) / (2 * step);
      second[level] = (ep - 2 * energy + em) / (step * step);
    }
    // This energy is quartic in nodal translations; Richardson removes the
    // exact centered-difference h² term, independent of either kernel path.
    const double df = (4 * first[1] - first[0]) / 3;
    const double ddf = (4 * second[1] - second[0]) / 3;
    EXPECT_NEAR(direction.Dot(response.Value().internal_force), df, 2e-10);
    EXPECT_NEAR(direction.Dot(paired.Value().high.internal_force), df, 2e-10);
    EXPECT_NEAR(direction.Dot(response.Value().tangent.Multiply(direction)),
                ddf, 2e-8);
    EXPECT_NEAR(direction.Dot(paired.Value().high.tangent.Multiply(direction)),
                ddf, 2e-8);
  }
}

TEST(Hex8Element, MeanVolumeWarpedLowOnlyStateRetainsWeightedCubicEnergy) {
  using namespace fesa::math_internal;
  auto element = Make(Warped());
  ASSERT_TRUE(element.HasValue());
  const double c = std::ldexp(1.0, -60);
  fesa::Vector high(24), low(24), zero(24), basis(24), minus(24);
  for (unsigned a = 0; a < 8; ++a)
    for (unsigned i = 0; i < 3; ++i) high[3 * a + i] = 1;
  low[1] = c;
  basis[1] = 1;
  minus[1] = -1;
  const auto linear = element.Value().EvaluateLinearHighLow(high, low);
  const auto nonlinear = element.Value().EvaluateNonlinearHighLow(high, low);
  const auto high_only = element.Value().EvaluateNonlinearHighLow(high, zero);
  ASSERT_TRUE(linear.HasValue());
  ASSERT_TRUE(nonlinear.HasValue());
  ASSERT_TRUE(high_only.HasValue());
  EXPECT_EQ(high_only.Value().high.strain_energy, 0);
  const double quadratic = WarpedQuadratureEnergies(basis, 2)[0];
  const double cubic = .5 * (WarpedQuadratureEnergies(basis, 2, true)[0] -
                             WarpedQuadratureEnergies(minus, 2, true)[0]);
  EXPECT_NEAR(linear.Value().high.strain_energy, quadratic * c * c,
              1e-12 * c * c);
  const auto actual_cubic = Subtract(
      {nonlinear.Value().high.strain_energy,
       nonlinear.Value().strain_energy_low},
      {linear.Value().high.strain_energy, linear.Value().strain_energy_low});
  EXPECT_NEAR(actual_cubic.high, cubic * c * c * c, 1e-11 * c * c * c);
  EXPECT_NE(actual_cubic.high, 0);
}
}  // namespace