#include "fesa/elements/hex8_element.h"

#include <gtest/gtest.h>

#include <cmath>

#include "fesa/materials/isotropic_linear_elastic_material.h"
#include "../../../src/fesa/math/high_low_arithmetic.h"
#include "../../../src/fesa/elements/hex8_geometry.h"

namespace {
using fesa::math_internal::HighLow;
using fesa::math_internal::Subtract;

fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) {
  fesa::ElementDofLayout layout;
  layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7};
  layout.components_per_node = {fesa::DofComponent::kUx,
                                fesa::DofComponent::kUy,
                                fesa::DofComponent::kUz};
  fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0},
                                             {1, 1, 0}, {0, 1, 0},
                                             {0, 0, 1}, {1, 0, 1},
                                             {1, 1, 1}, {0, 1, 1}}};
  for (auto& x : coordinates) x[0] += offset;
  fesa::IsotropicLinearElasticMaterial material("test", 12, nu, {});
  return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value();
}

TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) {
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1;
  for (auto a : {1, 2, 5, 6}) low[3 * a] = strain;
  const auto response = brick.EvaluateLinearHighLow(high, low);
  ASSERT_TRUE(response.HasValue());
  for (std::size_t a = 0; a < 8; ++a) {
    const bool positive = a == 1 || a == 2 || a == 5 || a == 6;
    EXPECT_NEAR(response.Value().high.internal_force[3 * a],
                (positive ? 3 : -3) * strain, 1e-13 * strain);
  }
  EXPECT_NEAR(response.Value().high.strain_energy, 6 * strain * strain,
              1e-13 * strain * strain);
  for (const auto& row : response.Value().high.points)
    EXPECT_NEAR(row.strain[0], strain, 1e-13 * strain);
}

TEST(Hex8HighLow, NonlinearLowGradientProducesSecondOrderForce) {
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1;
  for (auto a : {1, 2, 5, 6}) low[3 * a] = strain;
  const auto linear = brick.EvaluateLinearHighLow(high, low);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(linear.HasValue());
  ASSERT_TRUE(nonlinear.HasValue());
  const auto& l = linear.Value();
  const auto& n = nonlinear.Value();
  const auto increment = Subtract(
      HighLow{n.high.internal_force[3], n.internal_force_low[3]},
      HighLow{l.high.internal_force[3], l.internal_force_low[3]});
  EXPECT_NEAR(increment.high, 4.5 * strain * strain,
              1e-12 * strain * strain);
}

TEST(Hex8HighLow, VolumetricAndTransverseTermsRetainLowStrain) {
  auto brick = UnitBrick(0, 0.25);
  fesa::Vector high(24), low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) high[3*a] = 1;
  for (auto a : {1, 2, 5, 6}) low[3*a] = strain;
  const auto linear = brick.EvaluateLinearHighLow(high, low);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(linear.HasValue()); ASSERT_TRUE(nonlinear.HasValue());
  // Uniform axial strain: lambda=mu=4.8, each nodal face integral is 1/4.
  EXPECT_NEAR(linear.Value().high.internal_force[3], 3.6*strain, 1e-13*strain);
  EXPECT_NEAR(linear.Value().high.internal_force[7], 1.2*strain, 1e-13*strain);
  EXPECT_NEAR(linear.Value().high.strain_energy, 7.2*strain*strain,
              1e-13*strain*strain);
  const auto quadratic = Subtract(
      {nonlinear.Value().high.internal_force[3], nonlinear.Value().internal_force_low[3]},
      {linear.Value().high.internal_force[3], linear.Value().internal_force_low[3]});
  EXPECT_NEAR(quadratic.high, 5.4*strain*strain, 1e-12*strain*strain);
}

TEST(Hex8HighLow, RejectsInvalidLowShapeAndOverlappingComponents) {
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue());
  high[0] = 1;
  low[0] = 1;
  EXPECT_FALSE(brick.EvaluateNonlinearHighLow(high, low).HasValue());
}

TEST(Hex8HighLow, LinearAndNonlinearZeroStateHaveSameFullHessian) {
  auto brick = UnitBrick();
  fesa::Vector zero(24);
  const auto linear = brick.EvaluateLinearHighLow(zero, zero);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(zero, zero);
  ASSERT_TRUE(linear.HasValue());
  ASSERT_TRUE(nonlinear.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(linear.Value().high.tangent(i, j),
                nonlinear.Value().high.tangent(i, j));
      EXPECT_EQ(linear.Value().tangent_low(i, j),
                nonlinear.Value().tangent_low(i, j));
    }
}

TEST(Hex8HighLow, CachedLinearHessianIsIndependentOfTrialAndPriorFailure) {
  auto brick = UnitBrick(0, 0.25);
  fesa::Vector zero(24), high(24), low(24);
  high[3] = 0.1;
  low[3] = std::ldexp(1.0, -60);
  const auto reference = brick.EvaluateLinearHighLow(zero, zero);
  ASSERT_TRUE(reference.HasValue());
  EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue());
  const auto strained = brick.EvaluateLinearHighLow(high, low);
  ASSERT_TRUE(strained.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(reference.Value().high.tangent(i,j), strained.Value().high.tangent(i,j));
      EXPECT_EQ(reference.Value().tangent_low(i,j), strained.Value().tangent_low(i,j));
    }
  const auto nonlinear_zero = brick.EvaluateNonlinearHighLow(zero, zero);
  ASSERT_TRUE(nonlinear_zero.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(reference.Value().high.tangent(i,j), nonlinear_zero.Value().high.tangent(i,j));
      EXPECT_EQ(reference.Value().tangent_low(i,j), nonlinear_zero.Value().tangent_low(i,j));
    }
}

TEST(Hex8HighLow, CoupledBricksTransmitLowStrainAcrossSharedFace) {
  const auto left = UnitBrick(), right = UnitBrick(1);
  fesa::Vector high(24), left_low(24), right_low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) {
    high[3 * a] = 1;
    const bool positive = a == 1 || a == 2 || a == 5 || a == 6;
    left_low[3 * a] = positive ? strain : 0;
    right_low[3 * a] = positive ? 2 * strain : strain;
  }
  const auto l = left.EvaluateLinearHighLow(high, left_low);
  const auto r = right.EvaluateLinearHighLow(high, right_low);
  ASSERT_TRUE(l.HasValue());
  ASSERT_TRUE(r.HasValue());
  const std::array<std::size_t, 4> face_left{1, 2, 5, 6};
  const std::array<std::size_t, 4> face_right{0, 3, 4, 7};
  for (std::size_t a = 0; a < 4; ++a) {
    const auto i = 3 * face_left[a], j = 3 * face_right[a];
    const auto assembled = fesa::math_internal::Add(
        {l.Value().high.internal_force[i], l.Value().internal_force_low[i]},
        {r.Value().high.internal_force[j], r.Value().internal_force_low[j]});
    EXPECT_NEAR(assembled.high, 0, 1e-13 * strain);
  }
  EXPECT_NEAR(l.Value().high.strain_energy + r.Value().high.strain_energy,
              12 * strain * strain, 1e-13 * strain * strain);
}

TEST(Hex8HighLow, NodeZeroSignedDirectionMatchesForceAndEnergyDerivatives) {
  using namespace fesa::math_internal;
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  for (std::size_t a = 0; a < 8; ++a) {
    high[3 * a] = 1;
    low[3 * a] = (a == 1 || a == 2 || a == 5 || a == 6)
                     ? std::ldexp(1.0, -60) : 0;
  }
  const auto baseline = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(baseline.HasValue());
  const auto& base = baseline.Value();
  for (int exponent : {-65, -70, -75}) {
    const double step = std::ldexp(1.0, exponent);
    auto plus_high = high, minus_high = high;
    auto plus_low = low, minus_low = low;
    const auto plus = Add({high[0], low[0]}, HighLow{step});
    const auto minus = Subtract({high[0], low[0]}, HighLow{step});
    plus_high[0] = plus.high; plus_low[0] = plus.low;
    minus_high[0] = minus.high; minus_low[0] = minus.low;
    const auto positive = brick.EvaluateNonlinearHighLow(plus_high, plus_low);
    const auto negative = brick.EvaluateNonlinearHighLow(minus_high, minus_low);
    ASSERT_TRUE(positive.HasValue());
    ASSERT_TRUE(negative.HasValue());
    for (std::size_t i = 0; i < 24; ++i) {
      const auto derivative = Divide(Subtract(
          {positive.Value().high.internal_force[i], positive.Value().internal_force_low[i]},
          {negative.Value().high.internal_force[i], negative.Value().internal_force_low[i]}),
          HighLow{2 * step});
      const auto error = Subtract(derivative,
          {base.high.tangent(i, 0), base.tangent_low(i, 0)});
      EXPECT_LE(std::abs(error.high), 1e-25);
    }
    const auto gradient = Divide(Subtract(
        {positive.Value().high.strain_energy, positive.Value().strain_energy_low},
        {negative.Value().high.strain_energy, negative.Value().strain_energy_low}),
        HighLow{2 * step});
    const auto energy_error = Subtract(gradient,
        {base.high.internal_force[0], base.internal_force_low[0]});
    // Centered energy differences retain the StVK cubic term: O(step^2).
    // For this unit E=12 brick its directional coefficient is bounded by 32.
    EXPECT_LE(std::abs(energy_error.high), 32 * step * step + 1e-44);
  }
}

TEST(Hex8HighLow, GeometryCertificateIncludesLowOnlyCurrentCellWidth) {
  fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0},
                                             {1, 1, 0}, {0, 1, 0},
                                             {0, 0, 1}, {1, 0, 1},
                                             {1, 1, 1}, {0, 1, 1}}};
  fesa::Vector high(24), positive_low(24), negative_low(24);
  for (std::size_t a = 0; a < 8; ++a) {
    high[3 * a] = -coordinates[a][0];
    positive_low[3 * a] = std::ldexp(coordinates[a][0], -60);
    negative_low[3 * a] = -positive_low[3 * a];
  }
  using fesa::hex8_internal::Certificate;
  EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &positive_low),
            Certificate::kPositive);
  EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &negative_low),
            Certificate::kNonpositive);
  EXPECT_NE(fesa::hex8_internal::Certify(coordinates, &high),
            Certificate::kPositive);
}

TEST(Hex8HighLow, SignedZeroMaterialCoefficientHasIdenticalResponse) {
  auto positive = UnitBrick(0, 0.0), negative = UnitBrick(0, -0.0);
  fesa::Vector high(24), low(24);
  high[3] = 0.01;
  for (bool nonlinear : {false, true}) {
    const auto p = nonlinear ? positive.EvaluateNonlinearHighLow(high, low)
                             : positive.EvaluateLinearHighLow(high, low);
    const auto n = nonlinear ? negative.EvaluateNonlinearHighLow(high, low)
                             : negative.EvaluateLinearHighLow(high, low);
    ASSERT_TRUE(p.HasValue());
    ASSERT_TRUE(n.HasValue());
    EXPECT_EQ(p.Value().high.strain_energy, n.Value().high.strain_energy);
    EXPECT_EQ(p.Value().strain_energy_low, n.Value().strain_energy_low);
    for (std::size_t i = 0; i < 24; ++i) {
      EXPECT_EQ(p.Value().high.internal_force[i], n.Value().high.internal_force[i]);
      EXPECT_EQ(p.Value().internal_force_low[i], n.Value().internal_force_low[i]);
      for (std::size_t j = 0; j < 24; ++j) {
        EXPECT_EQ(p.Value().high.tangent(i,j), n.Value().high.tangent(i,j));
        EXPECT_EQ(p.Value().tangent_low(i,j), n.Value().tangent_low(i,j));
      }
    }
  }
}

TEST(Hex8HighLow, FiniteRotationRetainsTinyShearForceCovariance) {
  using namespace fesa::math_internal;
  auto brick = UnitBrick();
  fesa::Hex8Element::Coordinates x{{{0, 0, 0}, {1, 0, 0},
                                  {1, 1, 0}, {0, 1, 0},
                                  {0, 0, 1}, {1, 0, 1},
                                  {1, 1, 1}, {0, 1, 1}}};
  fesa::Vector zero(24), rotation(24), shear(24), rotated_shear(24);
  const double gamma = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) {
    // Q is an exact quarter turn; Q*(gamma*y,0,0)=(0,gamma*y,0).
    rotation[3*a] = HighLow{-x[a][1] - x[a][0]}.high;
    rotation[3*a+1] = x[a][0] - x[a][1];
    shear[3*a] = gamma*x[a][1];
    rotated_shear[3*a+1] = gamma*x[a][1];
  }
  fesa::Vector turned_high(24), turned_low(24);
  for (std::size_t j = 0; j < 24; ++j) {
    const auto full = Add(HighLow{rotation[j]}, HighLow{rotated_shear[j]});
    turned_high[j] = full.high; turned_low[j] = full.low;
  }
  const auto origin = brick.EvaluateNonlinearHighLow(zero, zero);
  const auto deformed = brick.EvaluateNonlinearHighLow(shear, zero);
  const auto turned = brick.EvaluateNonlinearHighLow(rotation, zero);
  const auto turned_deformed = brick.EvaluateNonlinearHighLow(turned_high, turned_low);
  ASSERT_TRUE(origin.HasValue()); ASSERT_TRUE(deformed.HasValue());
  ASSERT_TRUE(turned.HasValue()); ASSERT_TRUE(turned_deformed.HasValue());
  // The fixed F64 reference gradients leave coefficient-level rotation error.
  // Subtract the same rigid baseline; do not project that force to zero.
  for (std::size_t a = 0; a < 8; ++a) {
    std::array<HighLow, 3> baseline{}, rotated{};
    for (std::size_t i = 0; i < 3; ++i) {
      const auto j = 3*a+i;
      baseline[i] = Subtract(
          {deformed.Value().high.internal_force[j], deformed.Value().internal_force_low[j]},
          {origin.Value().high.internal_force[j], origin.Value().internal_force_low[j]});
      rotated[i] = Subtract(
          {turned_deformed.Value().high.internal_force[j], turned_deformed.Value().internal_force_low[j]},
          {turned.Value().high.internal_force[j], turned.Value().internal_force_low[j]});
    }
    EXPECT_LE(std::abs(Add(rotated[0], baseline[1]).high), 1e-13*gamma);
    EXPECT_LE(std::abs(Subtract(rotated[1], baseline[0]).high), 1e-13*gamma);
    EXPECT_LE(std::abs(Subtract(rotated[2], baseline[2]).high), 1e-13*gamma);
  }
}
}  // namespace