#include "fesa/analysis/linear_frequency_procedure.h"

#include <gtest/gtest.h>

#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
#include <utility>

#include "fesa/analysis/analysis_model.h"
#include "fesa/analysis/linear_static_procedure.h"
#include "fesa/analysis/prepared_analysis.h"
#include "fesa/assembly/parallel_for.h"
#include "fesa/assembly/sparse_assembler.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/model/domain.h"
#include "fesa/model/model_types.h"
#include "fesa/solvers/linear/mkl_pardiso_solver.h"

TEST(LinearFrequencyProcedureExecutor, RejectsDomainWithoutFrequencyStep) {
  fesa::ModelDefinition definition;
  definition.source_path = "frequency-missing.inp";
  definition.source_content_identity = "frequency-missing";
  auto domain = fesa::Domain::Create(std::move(definition));
  ASSERT_TRUE(domain.HasValue());

  const fesa::LinearFrequencyProcedureExecutor executor;
  const auto result = executor.Execute(domain.Value());

  ASSERT_FALSE(result.HasValue());
  EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kSolver);
  ASSERT_EQ(result.GetStatus().Diagnostics().size(), 1U);
  EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code,
            "invalid-frequency-procedure");
}

namespace {

std::string StatusText(const fesa::Status& status) {
  std::string text;
  for (const auto& d : status.Diagnostics())
    text += d.code + " " + d.message + "\n";
  return text;
}

fesa::Result<fesa::Domain> PanelDomain(bool rotated, bool free_shell,
                                       bool statics = false,
                                       std::size_t requested = 1U,
                                       double length_scale = 1.0,
                                       double density = 1.0,
                                       double modulus = 12.0) {
  const std::array<std::array<double, 3>, 4> points{{{{0.8, -2.0, -0.6}},
                                                     {{0.8, 2.0, -0.6}},
                                                     {{-0.8, 2.0, 0.6}},
                                                     {{-0.8, -2.0, 0.6}}}};
  std::ostringstream deck;
  deck << std::setprecision(17) << "*Part,name=Panel\n*Node\n";
  const double c = std::sqrt(0.5);
  for (std::size_t i = 0; i < points.size(); ++i) {
    const auto& p = points[i];
    deck << 11U * (i + 1U) << ','
         << length_scale * (rotated ? c * (p[0] - p[1]) : p[0]) << ','
         << length_scale * (rotated ? c * (p[0] + p[1]) : p[1]) << ','
         << length_scale * p[2] << '\n';
  }
  deck << R"inp(*Element,type=S4
71,11,22,33,44
*Elset,elset=Shell
71
*Shell Section,elset=Shell,material=Steel
)inp" << length_scale
       << R"inp(
*End Part
*Assembly,name=Assembly
*Instance,name=Panel-1,part=Panel
*End Instance
*Nset,nset=Fixed,instance=Panel-1
22,33,44
*Nset,nset=Tip,instance=Panel-1
11
*End Assembly
*Material,name=Steel
*Elastic
)inp" << modulus
       << ",0.\n";
  if (!statics) deck << "*Density\n" << density << '\n';
  if (!free_shell) deck << "*Boundary\nFixed,1,6\nTip,1,5\n";
  if (statics) {
    deck << "*Step,name=Static,nlgeom=NO\n*Static\n0.1,1.,0.01,1.\n*Cload\nTip,"
            "4,"
         << (-4.0 / 3.0) * (rotated ? c : 1.0) << "\nTip,5,"
         << (rotated ? (-4.0 / 3.0) * c : 0.0) << "\nTip,6,1.\n";
  } else
    deck << "*Step,name=Step-2,nlgeom=NO,perturbation\n*Frequency\n"
         << requested << '\n';
  deck << "*End Step\n";
  static std::atomic<unsigned long long> sequence{0U};
  const auto path =
      std::filesystem::temp_directory_path() /
      ("fesa-eigen-procedure-" +
       std::to_string(
           std::chrono::steady_clock::now().time_since_epoch().count()) +
       "-" + std::to_string(sequence.fetch_add(1U)) + ".inp");
  {
    std::ofstream stream{path};
    stream << deck.str();
  }
  auto prepared = fesa::PrepareAnalysisInput(path);
  std::error_code ignored;
  std::filesystem::remove(path, ignored);
  if (!prepared.HasValue())
    return fesa::Result<fesa::Domain>::Failure(prepared.GetStatus());
  return fesa::Result<fesa::Domain>::Success(
      std::move(prepared.Value().domain));
}

}  // namespace

TEST(LinearFrequencyProcedureExecutor,
     MatchesPartialBoundaryPencilAndFullResidualCovariance) {
  std::array<double, 24> original_residual{};
  for (const bool rotated : {false, true}) {
    auto domain = PanelDomain(rotated, false);
    ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus());
    const auto result =
        fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value());
    ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus());
    const auto& state = result.Value();
    ASSERT_EQ(state.modes.size(), 1U);
    EXPECT_EQ(state.mass_null_dimension, 0U);
    EXPECT_EQ(state.finite_dimension, 1U);
    EXPECT_NEAR(state.modes[0].eigenvalue, 10369.0 / 450.0, 1e-10);
    EXPECT_EQ(state.modes[0].classification, 1U);
    EXPECT_EQ(state.vectors.Rows(), 1U);
    EXPECT_EQ(state.vectors.Columns(), 24U);
    EXPECT_NEAR(state.vectors(0U, 5U), std::sqrt(50.0 / 3.0), 1e-12);
    for (std::size_t dof = 0; dof < 24U; ++dof)
      if (dof != 5U) EXPECT_EQ(state.vectors(0U, dof), 0.0);
    ASSERT_EQ(state.eigenpair_checks.size(), 1U);
    EXPECT_NEAR(state.eigenpair_checks[0].modal_mass, 1.0, 1e-12);
    EXPECT_LE(state.eigenpair_checks[0].eta_r, 1e-10);
    EXPECT_LE(state.eigenpair_checks[0].eta_energy, 1e-10);
    EXPECT_LE(state.eigenpair_checks[0].rho_b, 1e-10);
    EXPECT_EQ(state.backend.eigenvectors.Rows(), 0U);
    auto static_domain = PanelDomain(rotated, false, true);
    ASSERT_TRUE(static_domain.HasValue())
        << StatusText(static_domain.GetStatus());
    const fesa::SerialParallelFor serial;
    const auto model = fesa::AnalysisModel::Create(static_domain.Value());
    ASSERT_TRUE(model.HasValue());
    const auto dofs = fesa::DofManager::Create(model.Value());
    ASSERT_TRUE(dofs.HasValue());
    const auto stiffness = fesa::SparseAssembler::AssembleStiffness(
        model.Value(), dofs.Value(), serial);
    ASSERT_TRUE(stiffness.HasValue());
    fesa::Vector unit{24U};
    unit[5U] = 1.0;
    const double kff = stiffness.Value().Multiply(unit)[5U];
    EXPECT_NEAR(kff, 10369.0 / 7500.0, 1e-12);
    fesa::Vector displacement{24U};
    displacement[5U] = 1.0 / kff;
    EXPECT_NEAR(displacement[5U], 7500.0 / 10369.0, 1e-12);
    auto residual = stiffness.Value().Multiply(displacement);
    const double c = std::sqrt(0.5);
    residual[3U] -= (-4.0 / 3.0) * (rotated ? c : 1.0);
    residual[4U] -= rotated ? (-4.0 / 3.0) * c : 0.0;
    residual[5U] -= 1.0;
    EXPECT_NEAR(residual[5U], 0.0, 1e-12);
    for (std::size_t i = 0U; i < 24U; ++i) {
      if (!rotated)
        original_residual[i] = residual[i];
      else {
        const std::size_t base = (i / 3U) * 3U;
        const double expected =
            i % 3U == 0U
                ? c * (original_residual[base] - original_residual[base + 1U])
                : (i % 3U == 1U ? c * (original_residual[base] +
                                       original_residual[base + 1U])
                                : original_residual[i]);
        EXPECT_NEAR(residual[i], expected, 1e-12);
      }
    }
  }
}

// E2-C09 is a pencil oracle. The public static contract conflict is recorded
// separately; direct covariance does not certify spatial moment equilibrium.
TEST(LinearFrequencyProcedureExecutor,
     PartialBoundaryStaticRejectsUnbalancedDrillingMoment) {
  const auto domain = PanelDomain(false, false, true);
  ASSERT_TRUE(domain.HasValue());
  const fesa::SerialParallelFor serial;
  fesa::MklPardisoSolver solver;
  const auto result =
      fesa::LinearStaticProcedureExecutor{serial, solver}.Execute(
          domain.Value());
  ASSERT_FALSE(result.HasValue());
  EXPECT_NE(StatusText(result.GetStatus())
                .find("global-equilibrium-tolerance-failure"),
            std::string::npos);
}

TEST(LinearFrequencyProcedureExecutor,
     FreePanelRetainsSixRigidModesAndAllComputedEvidence) {
  auto domain = PanelDomain(false, true, false, 6U);
  ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus());
  const auto result =
      fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value());
  ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus());
  const auto& state = result.Value();
  EXPECT_EQ(state.mass_null_dimension, 4U);
  EXPECT_EQ(state.finite_dimension, 20U);
  ASSERT_EQ(state.modes.size(), 6U);
  EXPECT_GE(state.eigenpair_checks.size(), 6U);
  for (const auto& mode : state.modes) {
    EXPECT_EQ(mode.classification, 0U);
    EXPECT_EQ(mode.cluster_id, 1U);
  }
  EXPECT_LE(state.eta_gram, 1e-10);
  EXPECT_LE(state.max_eta_orth, 1e-10);
  EXPECT_LE(state.eta_mz, 1e-12);
  EXPECT_LE(state.eta_kz, 1e-12);
  for (const auto& check : state.eigenpair_checks) {
    EXPECT_LE(check.eta_r, 1e-10);
    EXPECT_LE(check.eta_norm, 1e-10);
    EXPECT_LE(check.eta_energy, 1e-10);
  }
}

TEST(LinearFrequencyProcedureExecutor,
     ExcessRequestedModesFailsWithSourceContext) {
  auto domain = PanelDomain(false, false, false, 2U);
  ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus());
  const auto result =
      fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value());
  ASSERT_FALSE(result.HasValue());
  EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kModel);
  ASSERT_FALSE(result.GetStatus().Diagnostics().empty());
  EXPECT_EQ(result.GetStatus().Diagnostics()[0].code,
            "insufficient-finite-modes");
  EXPECT_GT(result.GetStatus().Diagnostics()[0].location.line, 0U);
}

TEST(LinearFrequencyProcedureExecutor, PreservesDensityAndPhysicalUnitScaling) {
  const auto baseline_domain = PanelDomain(false, false);
  const auto dense_domain = PanelDomain(false, false, false, 1U, 1.0, 4.0);
  constexpr double length = 1024.0;
  const auto units_domain =
      PanelDomain(false, false, false, 1U, length,
                  1.0 / (length * length * length), 12.0 / length);
  ASSERT_TRUE(baseline_domain.HasValue());
  ASSERT_TRUE(dense_domain.HasValue());
  ASSERT_TRUE(units_domain.HasValue());
  const fesa::LinearFrequencyProcedureExecutor executor;
  const auto baseline = executor.Execute(baseline_domain.Value());
  const auto dense = executor.Execute(dense_domain.Value());
  const auto units = executor.Execute(units_domain.Value());
  ASSERT_TRUE(baseline.HasValue()) << StatusText(baseline.GetStatus());
  ASSERT_TRUE(dense.HasValue()) << StatusText(dense.GetStatus());
  ASSERT_TRUE(units.HasValue()) << StatusText(units.GetStatus());
  EXPECT_NEAR(
      dense.Value().modes[0].eigenvalue / baseline.Value().modes[0].eigenvalue,
      0.25, 1e-12);
  EXPECT_NEAR(dense.Value().vectors(0U, 5U) / baseline.Value().vectors(0U, 5U),
              0.5, 1e-12);
  EXPECT_NEAR(
      units.Value().modes[0].eigenvalue / baseline.Value().modes[0].eigenvalue,
      1.0, 1e-12);
  EXPECT_NEAR(
      units.Value().vectors(0U, 5U) * length / baseline.Value().vectors(0U, 5U),
      1.0, 1e-12);
  EXPECT_LE(units.Value().eigenpair_checks[0].eta_r, 1e-10);
  EXPECT_LE(units.Value().eigenpair_checks[0].eta_energy, 1e-10);
}