#include "fesa/analysis/solid_static_procedure.h"

#include <gtest/gtest.h>

#include <cmath>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>

#include "fesa/assembly/parallel_for.h"
#include "fesa/io/abaqus/domain_mapper.h"
#include "fesa/io/abaqus/input_reader.h"
#include "fesa/math/sparse_matrix.h"
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"

namespace {
const double kX[8][3] = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
                         {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
std::string Deck(const std::string& boundary, const std::string& load = "",
                 bool nonlinear = false, double h = 0.25, double period = 1.0) {
  std::ostringstream out;
  out << "*Part,name=P\n*Node\n";
  for (int a = 0; a < 8; ++a)
    out << a + 1 << ',' << kX[a][0] << ',' << kX[a][1] << ',' << kX[a][2]
        << '\n';
  out << "*Element,type=C3D8,elset=E\n99,1,2,3,4,5,6,7,8\n"
         "*Solid Section,elset=E,material=M\n,\n*End Part\n"
         "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
         "*End Assembly\n*Material,name=M\n*Elastic\n120,.25\n"
         "*Step,name=Solid,inc=100";
  if (nonlinear) out << ",nlgeom=YES";
  out << "\n*Static,direct\n" << h << ',' << period << '\n';
  if (!boundary.empty()) out << "*Boundary\n" << boundary;
  if (!load.empty()) out << "*Cload\n" << load;
  out << "*End Step\n";
  return out.str();
}
fesa::Result<fesa::Domain> Map(const std::string& text) {
  const auto path = std::filesystem::temp_directory_path() /
                    "fesa-solid-static-procedure.inp";
  {
    std::ofstream out{path};
    out << text;
  }
  auto syntax = fesa::AbaqusInputReader{}.Read(path);
  std::filesystem::remove(path);
  if (!syntax.HasValue())
    return fesa::Result<fesa::Domain>::Failure(syntax.GetStatus());
  return fesa::AbaqusDomainMapper{}.Map(syntax.Value());
}
std::string Message(const fesa::Status& status) {
  return status.Diagnostics().empty() ? "" : status.Diagnostics()[0].message;
}
fesa::Result<fesa::SolidStaticExecutionResult> ExecuteSolid(
    const fesa::Domain& domain) {
  fesa::SerialParallelFor serial;
  fesa::MklPardisoSolver linear;
  fesa::MklPardisoTangentSolver tangent;
  return fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
      .Execute(domain);
}
// Symmetry supports permit the exact affine uniaxial solution
// u=(e*x,-nu*e*y,-nu*e*z).
std::string Supports(double left = 0.0) {
  std::ostringstream out;
  for (int a = 0; a < 8; ++a) {
    if (kX[a][0] == 0) out << a + 1 << ",1,1," << left << '\n';
    if (kX[a][1] == 0) out << a + 1 << ",2,2,0\n";
    if (kX[a][2] == 0) out << a + 1 << ",3,3,0\n";
  }
  return out.str();
}
std::string Tension(double force) {
  std::ostringstream out;
  for (int a = 0; a < 8; ++a)
    if (kX[a][0] == 1) out << a + 1 << ",1," << force / 4.0 << '\n';
  return out.str();
}
std::string AffineBoundary(double ex, double shear = 0.0, double angle = 0.0) {
  std::ostringstream out;
  out.precision(17);
  for (int a = 0; a < 8; ++a) {
    const double x = (1 + ex) * kX[a][0] + shear * kX[a][1];
    const double y = kX[a][1];
    const double u[3] = {std::cos(angle) * x - std::sin(angle) * y - kX[a][0],
                         std::sin(angle) * x + std::cos(angle) * y - kX[a][1],
                         0};
    for (int i = 0; i < 3; ++i)
      out << a + 1 << ',' << i + 1 << ',' << i + 1 << ',' << u[i] << '\n';
  }
  return out.str();
}
class FailedTangent final : public fesa::TangentSolver {
 public:
  enum class Mode { kBackend, kIteration, kNonfinite, kGeometry };
  explicit FailedTangent(Mode mode, bool after_accepted = false)
      : mode_(mode), after_accepted_(after_accepted) {}
  fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
    // First increment is a genuine accepted solve, then inject a backend fault.
    ++calls_;
    if (mode_ == Mode::kBackend && calls_ > 3U)
      return fesa::Status::Failure(fesa::FailureCategory::kSolver,
                                   {{fesa::Severity::kError,
                                     "injected-backend",
                                     {},
                                     "STATIC",
                                     {},
                                     "Backend failed."}});
    return solver_.FactorizeFresh(matrix);
  }
  fesa::Status SolveCorrection(const fesa::Vector& rhs,
                               fesa::Vector& du) const override {
    if (mode_ == Mode::kBackend || (after_accepted_ && calls_ <= 2U))
      return solver_.SolveCorrection(rhs, du);
    for (std::size_t i = 0; i < du.Size(); ++i)
      du[i] = mode_ == Mode::kNonfinite
                  ? std::numeric_limits<double>::quiet_NaN()
              : mode_ == Mode::kGeometry ? -100.0
                                         : 0.0;
    return fesa::Status::Ok();
  }

 private:
  Mode mode_;
  bool after_accepted_;
  unsigned calls_{0U};
  fesa::MklPardisoTangentSolver solver_;
};

class ObservedIndefiniteTangent final : public fesa::TangentSolver {
 public:
  fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
    rayleigh_sum = 0.0;
    for (double value : matrix.Values()) rayleigh_sum += value;
    return solver_.FactorizeFresh(matrix);
  }
  fesa::Status SolveCorrection(const fesa::Vector& rhs,
                               fesa::Vector& du) const override {
    return solver_.SolveCorrection(rhs, du);
  }
  double rayleigh_sum{0.0};

 private:
  fesa::MklPardisoTangentSolver solver_;
};
}  // namespace

TEST(SolidStaticProcedure, LinearForceNonzeroSupportAndCompleteRecovery) {
  auto domain = Map(Deck(Supports(0.02), Tension(1.2)));
  ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 1U);
  const auto& frame = result.Value().state.AcceptedFrames()[0];
  EXPECT_EQ(frame.frame_index, 0U);
  EXPECT_DOUBLE_EQ(frame.load_factor, 1.0);
  ASSERT_EQ(frame.displacement.Size(), 24U);
  ASSERT_EQ(frame.solid_results.size(), 8U);
  for (std::size_t a = 0; a < 8; ++a) {
    EXPECT_NEAR(frame.displacement[3 * a], 0.02 + 0.01 * kX[a][0], 1e-13);
    EXPECT_NEAR(frame.displacement[3 * a + 1], -0.0025 * kX[a][1], 1e-13);
    EXPECT_NEAR(frame.displacement[3 * a + 2], -0.0025 * kX[a][2], 1e-13);
    if (kX[a][0] == 1)
      EXPECT_DOUBLE_EQ(frame.reaction[3 * a], 0.0);
    else
      EXPECT_NEAR(frame.reaction[3 * a], -0.3, 1e-13);
  }
  EXPECT_NEAR(frame.physical_strain_energy, 0.006, 1e-13);
  for (const auto& point : frame.solid_results) {
    EXPECT_NEAR(point.cauchy_stress[0], 1.2, 1e-12);
    EXPECT_NEAR(point.cauchy_stress[1], 0.0, 1e-12);
    EXPECT_NEAR(point.strain[0], 0.01, 1e-13);
  }
  for (double value : frame.global_equilibrium) EXPECT_NEAR(value, 0.0, 1e-12);
}

TEST(SolidStaticProcedure,
     AllConstrainedShearHasTensorStrainAndLinearMomentEvidence) {
  auto domain = Map(Deck(AffineBoundary(0.02, 0.03)));
  ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  const auto& frame = result.Value().state.CommittedState();
  for (const auto& point : frame.solid_results) {
    EXPECT_NEAR(point.strain[3], 0.015, 1e-13);
    EXPECT_NEAR(point.cauchy_stress[3], 1.44, 1e-12);
  }
  EXPECT_NEAR(frame.physical_strain_energy, 0.0504, 1e-13);
  EXPECT_GT(std::abs(frame.balance.current_internal_moment[2]), 1e-5);
  EXPECT_NEAR(frame.balance.current_internal_moment[2],
              frame.balance.linearization_moment[2], 1e-12);
  EXPECT_NEAR(frame.balance.reference_internal_moment[2], 0.0, 1e-12);
}

TEST(SolidStaticProcedure,
     NonlinearDirectScheduleScalesForceAndAbsoluteTargets) {
  auto domain = Map(Deck(Supports(0.02), Tension(1.2), true, 0.5, 2.0));
  ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 4U);
  for (std::size_t n = 0; n < 4; ++n) {
    const auto& frame = result.Value().state.AcceptedFrames()[n];
    EXPECT_EQ(frame.frame_index, n);
    EXPECT_EQ(frame.increment_index, n + 1U);
    EXPECT_DOUBLE_EQ(frame.load_factor, 0.25 * (n + 1));
    EXPECT_DOUBLE_EQ(frame.displacement[0], 0.02 * frame.load_factor);
    EXPECT_DOUBLE_EQ(frame.external_force[3], 0.3 * frame.load_factor);
    EXPECT_GT(frame.iteration_count, 1U);
    EXPECT_EQ(frame.is_factor_one, n == 3);
    for (double value : frame.balance.current_internal_moment)
      EXPECT_NEAR(value, 0.0, 1e-12);
  }
  ASSERT_FALSE(result.Value().state.IterationRecords().empty());
  EXPECT_TRUE(result.Value().state.IterationRecords().back().accepted);
}

TEST(SolidStaticProcedure, PurePrescribedNonlinearAndSmallLoadLinearLimit) {
  auto pure = Map(Deck(
      Supports() + "2,1,1,.01\n3,1,1,.01\n6,1,1,.01\n7,1,1,.01\n", "", true));
  ASSERT_TRUE(pure.HasValue());
  auto result = ExecuteSolid(pure.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  for (const auto& record : result.Value().state.IterationRecords())
    if (record.convergence_evaluated) {
      EXPECT_EQ(record.force_residual.branch,
                fesa::NewtonConvergenceBranch::kZeroScale);
      EXPECT_DOUBLE_EQ(record.force_residual.denominator_scale, 0.0);
    }
  auto linear = Map(Deck(Supports(), Tension(1.2e-5)));
  auto nonlinear = Map(Deck(Supports(), Tension(1.2e-5), true));
  ASSERT_TRUE(linear.HasValue());
  ASSERT_TRUE(nonlinear.HasValue());
  auto l = ExecuteSolid(linear.Value());
  auto n = ExecuteSolid(nonlinear.Value());
  ASSERT_TRUE(l.HasValue());
  ASSERT_TRUE(n.HasValue()) << Message(n.GetStatus());
  for (std::size_t i = 0; i < 24; ++i)
    EXPECT_NEAR(l.Value().state.CommittedState().displacement[i],
                n.Value().state.CommittedState().displacement[i], 2e-13);
}

TEST(SolidStaticProcedure, StrainedThenRotatedTLCauchyAndEnergyOracle) {
  const double angle = 1.2;
  auto domain = Map(Deck(AffineBoundary(0.03, 0, angle), "", true));
  ASSERT_TRUE(domain.HasValue());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  const auto& frame = result.Value().state.CommittedState();
  const double green = (1.03 * 1.03 - 1) / 2;
  const double s0 = 144 * green * 1.03;
  const double s1 = 48 * green / 1.03;
  for (const auto& point : frame.solid_results) {
    EXPECT_NEAR(point.strain[0], green, 1e-13);
    EXPECT_NEAR(
        point.cauchy_stress[0],
        std::pow(std::cos(angle), 2) * s0 + std::pow(std::sin(angle), 2) * s1,
        1e-11);
    EXPECT_NEAR(point.cauchy_stress[3],
                std::sin(angle) * std::cos(angle) * (s0 - s1), 1e-11);
  }
  EXPECT_NEAR(frame.physical_strain_energy, 72 * green * green, 1e-12);
  for (double value : frame.balance.current_internal_moment)
    EXPECT_NEAR(value, 0.0, 1e-12);
}

TEST(SolidStaticProcedure, FailedBackendKeepsActualAcceptedHistory) {
  auto domain = Map(Deck(Supports(), Tension(1.2), true));
  ASSERT_TRUE(domain.HasValue());
  fesa::SerialParallelFor serial;
  fesa::MklPardisoSolver linear;
  FailedTangent tangent{FailedTangent::Mode::kBackend};
  fesa::SolidStaticState state;
  auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                    .ExecuteInto(domain.Value(), state);
  ASSERT_FALSE(status.IsOk());
  ASSERT_EQ(state.AcceptedFrames().size(), 1U);
  EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.25);
  EXPECT_EQ(state.CommittedState().displacement.Size(), 24U);
  EXPECT_EQ(state.CommittedState().solid_results.size(), 8U);
  EXPECT_EQ(state.IterationRecords().back().diagnostic_code,
            "injected-backend");
  auto bad = state.CommittedState();
  bad.frame_index = 1;
  bad.increment_index = 2;
  bad.load_factor = 0.5;
  bad.displacement[0] = std::numeric_limits<double>::quiet_NaN();
  EXPECT_FALSE(state.CommitAcceptedFrame(bad).IsOk());
  bad = state.CommittedState();
  bad.frame_index = 1;
  bad.increment_index = 2;
  bad.load_factor = .5;
  bad.solid_results.pop_back();
  EXPECT_FALSE(state.CommitAcceptedFrame(bad).IsOk());
  EXPECT_EQ(state.AcceptedFrames().size(), 1U);
  EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.25);
}

TEST(SolidStaticProcedure,
     GeometryIterationAndNonfiniteFailuresPreserveBaseline) {
  auto domain = Map(Deck(Supports(), Tension(1.2), true));
  ASSERT_TRUE(domain.HasValue());
  for (auto mode :
       {FailedTangent::Mode::kGeometry, FailedTangent::Mode::kIteration,
        FailedTangent::Mode::kNonfinite}) {
    fesa::SerialParallelFor serial;
    fesa::MklPardisoSolver linear;
    FailedTangent tangent{mode};
    fesa::SolidStaticState state;
    auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                      .ExecuteInto(domain.Value(), state);
    EXPECT_FALSE(status.IsOk());
    EXPECT_TRUE(state.AcceptedFrames().empty());
    EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.0);
    ASSERT_EQ(state.CommittedState().displacement.Size(), 24U);
    for (std::size_t i = 0; i < 24; ++i)
      EXPECT_DOUBLE_EQ(state.CommittedState().displacement[i], 0.0);
    ASSERT_FALSE(state.IterationRecords().empty());
    EXPECT_FALSE(state.IterationRecords().back().accepted);
    EXPECT_FALSE(state.IterationRecords().back().diagnostic_code.empty());
    if (mode == FailedTangent::Mode::kIteration)
      EXPECT_EQ(state.IterationRecords().back().iteration, 100U);
  }
}

TEST(SolidStaticProcedure, UnconstrainedSolidFailsWithoutRegularization) {
  auto domain = Map(Deck("", Tension(1.2)));
  ASSERT_TRUE(domain.HasValue());
  auto result = ExecuteSolid(domain.Value());
  EXPECT_FALSE(result.HasValue());
}

TEST(SolidStaticProcedure,
     EightElementAffinePatchHasSolvedInteriorAndStableRows) {
  std::ostringstream out;
  out << "*Part,name=P\n*Node\n";
  const auto id = [](int x, int y, int z) { return 1 + x + 3 * y + 9 * z; };
  for (int z = 0; z < 3; ++z)
    for (int y = 0; y < 3; ++y)
      for (int x = 0; x < 3; ++x)
        out << id(x, y, z) << ',' << x << ',' << y << ',' << z << '\n';
  out << "*Element,type=C3D8,elset=E\n";
  int e = 0;
  for (int z = 0; z < 2; ++z)
    for (int y = 0; y < 2; ++y)
      for (int x = 0; x < 2; ++x)
        out << ++e << ',' << id(x, y, z) << ',' << id(x + 1, y, z) << ','
            << id(x + 1, y + 1, z) << ',' << id(x, y + 1, z) << ','
            << id(x, y, z + 1) << ',' << id(x + 1, y, z + 1) << ','
            << id(x + 1, y + 1, z + 1) << ',' << id(x, y + 1, z + 1) << '\n';
  out << "*Solid Section,elset=E,material=M\n,\n*End Part\n"
         "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
         "*End Assembly\n*Material,name=M\n*Elastic\n120,.25\n"
         "*Step,name=Patch\n*Static\n1,1,1,1\n*Boundary\n";
  for (int z = 0; z < 3; ++z)
    for (int y = 0; y < 3; ++y)
      for (int x = 0; x < 3; ++x) {
        if (x == 1 && y == 1 && z == 1) continue;
        out << id(x, y, z) << ",1,1," << 0.02 * x + 0.03 * y << '\n';
        out << id(x, y, z) << ",2,3,0\n";
      }
  out << "*End Step\n";
  auto domain = Map(out.str());
  ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  const auto& frame = result.Value().state.CommittedState();
  ASSERT_EQ(frame.solid_results.size(), 64U);
  EXPECT_NEAR(frame.displacement[3 * 13], 0.05, 1e-13);
  EXPECT_NEAR(frame.displacement[3 * 13 + 1], 0.0, 1e-13);
  EXPECT_NEAR(frame.displacement[3 * 13 + 2], 0.0, 1e-13);
  EXPECT_NEAR(frame.physical_strain_energy, 8 * 0.0504, 1e-12);
  for (std::size_t g = 0; g < 64; ++g) {
    EXPECT_EQ(frame.solid_results[g].element, g / 8);
    EXPECT_EQ(frame.solid_results[g].integration_point,
              static_cast<int>(g % 8 + 1));
    EXPECT_NEAR(frame.solid_results[g].cauchy_stress[0], 2.88, 1e-12);
    EXPECT_NEAR(frame.solid_results[g].cauchy_stress[3], 1.44, 1e-12);
  }
  fesa::TbbParallelFor parallel;
  fesa::MklPardisoSolver linear;
  fesa::MklPardisoTangentSolver tangent;
  auto parallel_result =
      fesa::SolidStaticProcedureExecutor(parallel, linear, tangent)
          .Execute(domain.Value());
  ASSERT_TRUE(parallel_result.HasValue());
  const auto& parallel_frame = parallel_result.Value().state.CommittedState();
  EXPECT_DOUBLE_EQ(parallel_frame.physical_strain_energy,
                   frame.physical_strain_energy);
  for (std::size_t i = 0; i < frame.displacement.Size(); ++i) {
    EXPECT_DOUBLE_EQ(parallel_frame.displacement[i], frame.displacement[i]);
    EXPECT_DOUBLE_EQ(parallel_frame.internal_force[i], frame.internal_force[i]);
    EXPECT_DOUBLE_EQ(parallel_frame.residual[i], frame.residual[i]);
  }
  for (std::size_t g = 0; g < 64; ++g) {
    EXPECT_EQ(parallel_frame.solid_results[g].cauchy_stress,
              frame.solid_results[g].cauchy_stress);
    EXPECT_EQ(parallel_frame.solid_results[g].strain,
              frame.solid_results[g].strain);
  }
}

TEST(SolidStaticProcedure,
     LaterInvertedPrescribedTrialPreservesThreeAcceptedFrames) {
  auto domain = Map(Deck(AffineBoundary(-1.2), "", true));
  ASSERT_TRUE(domain.HasValue());
  fesa::SerialParallelFor serial;
  fesa::MklPardisoSolver linear;
  fesa::MklPardisoTangentSolver tangent;
  fesa::SolidStaticState state;
  auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                    .ExecuteInto(domain.Value(), state);
  ASSERT_FALSE(status.IsOk());
  ASSERT_EQ(state.AcceptedFrames().size(), 3U);
  EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.75);
  EXPECT_NEAR(state.CommittedState().displacement[3], -0.9, 1e-15);
  ASSERT_EQ(state.CommittedState().solid_results.size(), 8U);
  EXPECT_EQ(state.IterationRecords().back().increment_index, 4U);
  EXPECT_FALSE(state.IterationRecords().back().accepted);
}

TEST(SolidStaticProcedure, PartiallySupportedRigidMechanismAlsoFailsUnloaded) {
  for (const auto& load : {std::string{}, Tension(1.2)}) {
    auto domain = Map(Deck("1,1,3,0\n", load));
    ASSERT_TRUE(domain.HasValue());
    auto result = ExecuteSolid(domain.Value());
    EXPECT_FALSE(result.HasValue());
    if (!result.HasValue())
      EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kSolver);
  }
}

TEST(SolidStaticProcedure, SerialAndTbbFramesAndFailureEvidenceAreIdentical) {
  auto domain = Map(Deck(Supports(), Tension(1.2), true));
  ASSERT_TRUE(domain.HasValue());
  fesa::SerialParallelFor serial;
  fesa::TbbParallelFor parallel;
  fesa::MklPardisoSolver linear_a, linear_b;
  fesa::MklPardisoTangentSolver tangent_a, tangent_b;
  auto a = fesa::SolidStaticProcedureExecutor(serial, linear_a, tangent_a)
               .Execute(domain.Value());
  auto b = fesa::SolidStaticProcedureExecutor(parallel, linear_b, tangent_b)
               .Execute(domain.Value());
  ASSERT_TRUE(a.HasValue());
  ASSERT_TRUE(b.HasValue());
  ASSERT_EQ(a.Value().state.AcceptedFrames().size(),
            b.Value().state.AcceptedFrames().size());
  for (std::size_t n = 0; n < a.Value().state.AcceptedFrames().size(); ++n) {
    const auto& af = a.Value().state.AcceptedFrames()[n];
    const auto& bf = b.Value().state.AcceptedFrames()[n];
    EXPECT_DOUBLE_EQ(af.physical_strain_energy, bf.physical_strain_energy);
    EXPECT_EQ(af.global_equilibrium, bf.global_equilibrium);
    for (std::size_t i = 0; i < 24; ++i) {
      EXPECT_DOUBLE_EQ(af.displacement[i], bf.displacement[i]);
      EXPECT_DOUBLE_EQ(af.residual[i], bf.residual[i]);
      EXPECT_DOUBLE_EQ(af.reaction[i], bf.reaction[i]);
    }
    for (std::size_t g = 0; g < 8; ++g) {
      EXPECT_EQ(af.solid_results[g].cauchy_stress,
                bf.solid_results[g].cauchy_stress);
      EXPECT_EQ(af.solid_results[g].strain, bf.solid_results[g].strain);
    }
  }
  auto invalid = Map(Deck(AffineBoundary(-1.2), "", true));
  ASSERT_TRUE(invalid.HasValue());
  fesa::SolidStaticState sa, sb;
  auto as = fesa::SolidStaticProcedureExecutor(serial, linear_a, tangent_a)
                .ExecuteInto(invalid.Value(), sa);
  auto bs = fesa::SolidStaticProcedureExecutor(parallel, linear_b, tangent_b)
                .ExecuteInto(invalid.Value(), sb);
  ASSERT_FALSE(as.IsOk());
  ASSERT_FALSE(bs.IsOk());
  EXPECT_EQ(as.Diagnostics()[0].code, bs.Diagnostics()[0].code);
  EXPECT_EQ(sa.AcceptedFrames().size(), sb.AcceptedFrames().size());
  EXPECT_EQ(sa.IterationRecords().back().stage,
            sb.IterationRecords().back().stage);
  EXPECT_DOUBLE_EQ(sa.CommittedState().physical_strain_energy,
                   sb.CommittedState().physical_strain_energy);
}

TEST(SolidStaticProcedure,
     InvalidRecoveryAndRepeatedExecutionKeepAcceptedFrame) {
  auto domain = Map(Deck(Supports(), Tension(1.2), true));
  ASSERT_TRUE(domain.HasValue());
  fesa::SerialParallelFor serial;
  fesa::MklPardisoSolver linear;
  FailedTangent tangent{FailedTangent::Mode::kBackend};
  fesa::SolidStaticState state;
  ASSERT_FALSE(fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                   .ExecuteInto(domain.Value(), state)
                   .IsOk());
  ASSERT_EQ(state.AcceptedFrames().size(), 1U);
  const auto valid = state.CommittedState();
  for (int mutation = 0; mutation < 5; ++mutation) {
    auto bad = valid;
    // Supply the next valid identity so rejection must examine recovery.
    bad.frame_index = 1U;
    bad.increment_index = 2U;
    bad.load_factor = .5;
    if (mutation == 0)
      bad.solid_results[0].cauchy_stress[2] =
          std::numeric_limits<double>::infinity();
    if (mutation == 1) bad.solid_results[0].integration_point = 2;
    if (mutation == 2) bad.solid_results[0].reference_coordinates[0] += 1;
    if (mutation == 3) bad.reaction[0] += 1;
    if (mutation == 4) bad.nodal_source_ids[0].source_label = 999;
    EXPECT_FALSE(state.CommitAcceptedFrame(bad).IsOk());
    ASSERT_EQ(state.AcceptedFrames().size(), 1U);
    EXPECT_DOUBLE_EQ(state.CommittedState().physical_strain_energy,
                     valid.physical_strain_energy);
    EXPECT_EQ(state.CommittedState().solid_results[0].cauchy_stress,
              valid.solid_results[0].cauchy_stress);
  }
  EXPECT_FALSE(fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                   .ExecuteInto(domain.Value(), state)
                   .IsOk());
  ASSERT_EQ(state.AcceptedFrames().size(), 1U);
}

TEST(SolidStaticProcedure, StressedIndefiniteFreeTangentUsesCapableBackend) {
  // F=diag(.2,1,1): the homogeneous transverse right-face mode has
  // second variation mu + S11 = 48 - 69.12 < 0. It is not an SPD tangent.
  std::ostringstream boundary;
  std::ostringstream load;
  for (int a = 0; a < 8; ++a) {
    boundary << a + 1 << ",1,1," << -.8 * kX[a][0] << '\n';
    boundary << a + 1 << ",3,3,0\n";
    if (kX[a][0] == 0) boundary << a + 1 << ",2,2,0\n";
    // P22 = S22 = lambda*E11 = -23.04; each Y-face node gets P22/4.
    // Hold the transverse affine state with its independently known traction.
    else
      load << a + 1 << ",2," << (kX[a][1] == 0 ? 5.76 : -5.76) << '\n';
  }
  auto domain = Map(Deck(boundary.str(), load.str(), true, 1.0));
  ASSERT_TRUE(domain.HasValue());
  fesa::SerialParallelFor serial;
  fesa::MklPardisoSolver linear;
  ObservedIndefiniteTangent tangent;
  auto result = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                    .Execute(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  EXPECT_NEAR(tangent.rayleigh_sum, -21.12, 1e-12);
  ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 1U);
  const auto& frame = result.Value().state.CommittedState();
  for (std::size_t a = 0; a < 8; ++a)
    EXPECT_NEAR(frame.displacement[3 * a + 1], 0.0, 1e-13);
  EXPECT_NEAR(frame.solid_results[0].strain[0], -.48, 1e-13);
}

TEST(SolidStaticProcedure,
     FullyConstrainedZeroDisplacementRetainsSupportLoadResidual) {
  auto domain = Map(Deck(AffineBoundary(0.0), "1,1,7\n"));
  ASSERT_TRUE(domain.HasValue());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue());
  const auto& frame = result.Value().state.CommittedState();
  EXPECT_DOUBLE_EQ(frame.internal_force[0], 0.0);
  EXPECT_DOUBLE_EQ(frame.residual[0], -7.0);
  EXPECT_DOUBLE_EQ(frame.reaction[0], -7.0);
  for (double value : frame.global_equilibrium) EXPECT_DOUBLE_EQ(value, 0.0);
  EXPECT_DOUBLE_EQ(frame.physical_strain_energy, 0.0);
  ASSERT_EQ(frame.solid_results.size(), 8U);
}

TEST(SolidStaticProcedure, PureLargeRigidRotationEndsAtZeroStrainEnergy) {
  auto domain = Map(Deck(AffineBoundary(0.0, 0.0, 1.2), "", true));
  ASSERT_TRUE(domain.HasValue());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  const auto& frame = result.Value().state.CommittedState();
  EXPECT_NEAR(frame.physical_strain_energy, 0.0, 1e-25);
  for (const auto& point : frame.solid_results) {
    for (double value : point.strain) EXPECT_NEAR(value, 0.0, 1e-14);
    for (double value : point.cauchy_stress) EXPECT_NEAR(value, 0.0, 1e-12);
  }
}

TEST(SolidStaticProcedure, OrdinaryStaticKeepsExistingFinalShortFactorPath) {
  std::string deck = Deck(AffineBoundary(.03), "", true);
  const std::string from = "*Static,direct\n0.25,1";
  const auto position = deck.find(from);
  ASSERT_NE(position, std::string::npos);
  deck.replace(position, from.size(), "*Static\n.3,1,.1,.3");
  auto domain = Map(deck);
  ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
  auto result = ExecuteSolid(domain.Value());
  ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
  ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 4U);
  const double factors[4] = {.3, .6, .9, 1.0};
  for (std::size_t i = 0; i < 4; ++i)
    EXPECT_NEAR(result.Value().state.AcceptedFrames()[i].load_factor,
                factors[i], 1e-15);
  EXPECT_DOUBLE_EQ(result.Value().state.CommittedState().load_factor, 1.0);
}

namespace {
// Actual slender 3D cantilever. The beam limit is a useful trend oracle,
// not an assertion that the coarse 3D discretization equals Euler theory.
std::string Cantilever(int nx, double poisson) {
  const int ny = 2, nz = 2;
  const auto id = [nx](int x, int y, int z) {
    return 1 + x + (nx + 1) * (y + 3 * z);
  };
  std::ostringstream out;
  out << "*Part,name=P\n*Node\n";
  for (int z = 0; z <= nz; ++z)
    for (int y = 0; y <= ny; ++y)
      for (int x = 0; x <= nx; ++x)
        out << id(x, y, z) << ',' << 8.0 * x / nx << ',' << .5 * y << ','
            << .5 * z << '\n';
  out << "*Element,type=C3D8,elset=E\n";
  int e = 0;
  for (int z = 0; z < nz; ++z)
    for (int y = 0; y < ny; ++y)
      for (int x = 0; x < nx; ++x)
        out << ++e << ',' << id(x, y, z) << ',' << id(x + 1, y, z) << ','
            << id(x + 1, y + 1, z) << ',' << id(x, y + 1, z) << ','
            << id(x, y, z + 1) << ',' << id(x + 1, y, z + 1) << ','
            << id(x + 1, y + 1, z + 1) << ',' << id(x, y + 1, z + 1) << '\n';
  out << "*Solid Section,elset=E,material=M\n,\n*End Part\n"
         "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n"
         "*End Assembly\n*Material,name=M\n*Elastic\n120,"
      << poisson << '\n';
  out << "*Step,name=Bending\n*Static\n1,1,1,1\n*Boundary\n";
  for (int z = 0; z <= nz; ++z)
    for (int y = 0; y <= ny; ++y) out << id(0, y, z) << ",1,3,0\n";
  out << "*Cload\n";
  for (int z = 0; z <= nz; ++z)
    for (int y = 0; y <= ny; ++y) {
      const double weight = (y == 1 ? .5 : .25) * (z == 1 ? .5 : .25);
      out << id(nx, y, z) << ",2," << 1e-4 * weight << '\n';
    }
  out << "*End Step\n";
  return out.str();
}
}  // namespace
TEST(SolidStaticProcedure,
     BoundedSpatialRefinementReportsFullIntegrationLocking) {
  const double beam_limit = 4 * 1e-4 * 8 * 8 * 8 / 120;
  double previous[2] = {0, 0};
  double final[2] = {0, 0};
  for (int material = 0; material < 2; ++material) {
    const double poisson = material == 0 ? .25 : .49;
    for (int nx : {2, 4, 8}) {
      auto domain = Map(Cantilever(nx, poisson));
      ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
      auto result = ExecuteSolid(domain.Value());
      ASSERT_TRUE(result.HasValue()) << Message(result.GetStatus());
      const auto& frame = result.Value().state.CommittedState();
      // End-center index x=nx,y=1,z=1 in node-major source ordering.
      const auto node = static_cast<std::size_t>(nx + 4 * (nx + 1));
      const double tip = frame.displacement[3 * node + 1];
      std::cout << "solid-static-refinement nu=" << poisson << " nx=" << nx
                << " tip=" << tip << " beam-limit=" << beam_limit << '\n';
      EXPECT_GT(tip, previous[material]);
      EXPECT_LT(tip, 1.05 * beam_limit);
      previous[material] = tip;
      final[material] = tip;
    }
  }
  // Near-incompressibility remains stiffer with this full-integration mesh.
  EXPECT_LT(final[1], final[0]);
  EXPECT_GT(std::abs(final[1] - beam_limit), std::abs(final[0] - beam_limit));
}

TEST(SolidStaticProcedure,
     LaterIterationNonfiniteAndGeometryFailureKeepEveryAcceptedQuantity) {
  auto domain = Map(Deck(Supports(), Tension(1.2), true));
  ASSERT_TRUE(domain.HasValue());
  auto successful = ExecuteSolid(domain.Value());
  ASSERT_TRUE(successful.HasValue());
  const auto& expected = successful.Value().state.AcceptedFrames()[0];
  ASSERT_EQ(expected.iteration_count, 2U);
  for (auto mode :
       {FailedTangent::Mode::kGeometry, FailedTangent::Mode::kIteration,
        FailedTangent::Mode::kNonfinite}) {
    fesa::SerialParallelFor serial;
    fesa::MklPardisoSolver linear;
    FailedTangent tangent{mode, true};
    fesa::SolidStaticState state;
    auto status = fesa::SolidStaticProcedureExecutor(serial, linear, tangent)
                      .ExecuteInto(domain.Value(), state);
    ASSERT_FALSE(status.IsOk());
    ASSERT_EQ(state.AcceptedFrames().size(), 1U);
    const auto& actual = state.CommittedState();
    EXPECT_EQ(actual.global_equilibrium, expected.global_equilibrium);
    EXPECT_DOUBLE_EQ(actual.physical_strain_energy,
                     expected.physical_strain_energy);
    for (std::size_t i = 0; i < 24; ++i) {
      EXPECT_DOUBLE_EQ(actual.displacement[i], expected.displacement[i]);
      EXPECT_DOUBLE_EQ(actual.internal_force[i], expected.internal_force[i]);
      EXPECT_DOUBLE_EQ(actual.external_force[i], expected.external_force[i]);
      EXPECT_DOUBLE_EQ(actual.residual[i], expected.residual[i]);
      EXPECT_DOUBLE_EQ(actual.reaction[i], expected.reaction[i]);
    }
    for (std::size_t g = 0; g < 8; ++g) {
      EXPECT_EQ(actual.solid_results[g].strain,
                expected.solid_results[g].strain);
      EXPECT_EQ(actual.solid_results[g].cauchy_stress,
                expected.solid_results[g].cauchy_stress);
    }
    EXPECT_EQ(state.IterationRecords().back().increment_index, 2U);
    EXPECT_FALSE(state.IterationRecords().back().accepted);
  }
}

TEST(SolidStaticProcedure,
     ExecutesBothUnmodifiedStaticSourceDecksThroughDirectApi) {
  auto root = std::filesystem::path{__FILE__}.parent_path();
  for (int i = 0; i < 3; ++i) root = root.parent_path();
  for (const std::string name :
       {"hexa_linearstatic.inp", "hexa_nonlinearstatic.inp"}) {
    auto syntax =
        fesa::AbaqusInputReader{}.Read(root / "reference/hexa" / name);
    ASSERT_TRUE(syntax.HasValue()) << Message(syntax.GetStatus());
    auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value());
    ASSERT_TRUE(domain.HasValue()) << Message(domain.GetStatus());
    auto result = ExecuteSolid(domain.Value());
    ASSERT_TRUE(result.HasValue())
        << name << ' ' << Message(result.GetStatus());
    const auto& state = result.Value().state;
    const auto& frame = state.CommittedState();
    EXPECT_DOUBLE_EQ(frame.load_factor, 1.0);
    EXPECT_EQ(frame.increment_index, 10U);
    EXPECT_EQ(frame.displacement.Size(), 442U * 3U);
    EXPECT_EQ(frame.solid_results.size(), 300U * 8U);
    EXPECT_EQ(state.AcceptedFrames().size(),
              domain.Value().Steps()[0].IsNonlinearStatic() ? 10U : 1U);
    std::cout << "solid-static-source " << name
              << " frames=" << state.AcceptedFrames().size()
              << " iterations=" << state.IterationRecords().size()
              << " energy=" << frame.physical_strain_energy << '\n';
  }
}