#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>

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 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, 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());
}
TEST(Hex8Element, RejectsRoundedLossOfPositiveVolumetricModulus) {
  // Binary64 Lame arithmetic at this valid nu yields lambda+2mu/3 == 0
  // even though exact E/[3(1-2nu)] is approximately1/9, which is positive.
  const double nu = std::nextafter(-1.0, 0.0);
  auto element = Make(Brick(), 1.0, nu);
  ASSERT_FALSE(element.HasValue());
  EXPECT_NE(element.GetStatus().Diagnostics()[0].message.find("bulk"),
            std::string::npos);
}
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) {
  const auto rule = LegendreRule(count);
  std::array<double, 2> result{};
  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];
        double trace = 0.0, strain_squared = 0.0, velocity_squared = 0.0;
        for (std::size_t i = 0; i < 3; ++i) {
          trace += gradient[i][i];
          velocity_squared += velocity[i] * velocity[i];
          for (std::size_t j = 0; j < 3; ++j) {
            const double strain = .5 * (gradient[i][j] + gradient[j][i]);
            strain_squared += strain * strain;
          }
        }
        const double weight = xi[1] * eta[1] * zeta[1] * det;
        result[0] += weight * (24.0 * trace * trace + 48.0 * strain_squared);
        result[1] += weight * velocity_squared;  // .5*rho*v^2 with rho=2.
      }
  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';
}
}  // namespace