#include "fesa/fesa.hpp"

#include <cmath>
#include <functional>
#include <iostream>
#include <type_traits>
#include <stdexcept>
#include <string>
#include <vector>

static_assert(std::is_same_v<fesa::Real, double>, "Real must remain double");
static_assert(std::is_same_v<fesa::GlobalId, std::int64_t>, "GlobalId must remain int64");
static_assert(std::is_same_v<fesa::LocalIndex, std::int64_t>, "LocalIndex must remain int64");
static_assert(std::is_same_v<fesa::EquationId, std::int64_t>, "EquationId must remain int64");
static_assert(std::is_same_v<fesa::SparseIndex, std::int64_t>, "SparseIndex must remain int64");

namespace {

using TestFn = std::function<void()>;

struct TestCase {
  std::string name;
  TestFn fn;
};

std::vector<TestCase>& registry() {
  static std::vector<TestCase> tests;
  return tests;
}

struct RegisterTest {
  RegisterTest(std::string name, TestFn fn) {
    registry().push_back({std::move(name), std::move(fn)});
  }
};

#define FESA_TEST(name)                                      \
  void name();                                               \
  RegisterTest register_##name(#name, name);                 \
  void name()

#define FESA_CHECK(expr)                                                                 \
  do {                                                                                   \
    if (!(expr)) {                                                                       \
      throw std::runtime_error(std::string("check failed: ") + #expr);                  \
    }                                                                                    \
  } while (false)

#define FESA_CHECK_NEAR(actual, expected, tol)                                            \
  do {                                                                                   \
    const auto actual_value = (actual);                                                   \
    const auto expected_value = (expected);                                               \
    if (std::fabs(actual_value - expected_value) > (tol)) {                              \
      throw std::runtime_error(std::string("near check failed: ") + #actual);           \
    }                                                                                    \
  } while (false)

std::string sourceRoot() {
#ifdef FESA_SOURCE_DIR
  return FESA_SOURCE_DIR;
#else
  return ".";
#endif
}

std::string phase1Input() {
  return R"inp(
*Node
1, 0, 0, 0
2, 1, 0, 0
3, 1, 1, 0
4, 0, 1, 0
*Element, type=S4, elset=EALL
1, 1, 2, 3, 4
*Nset, nset=LEFT
1, 4
*Nset, nset=RIGHT
2, 3
*Elset, elset=EALL
1
*Material, name=STEEL
*Elastic
1000.0, 0.3
*Shell Section, elset=EALL, material=STEEL
0.1
*Boundary
LEFT, 1, 6, 0
RIGHT, 1, 2, 0
RIGHT, 4, 6, 0
*Cload
2, 3, -1
3, 3, -1
*Step, name=Step-1
*Static
*End Step
)inp";
}

fesa::Domain parsedPhase1Domain() {
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(phase1Input());
  FESA_CHECK(parsed.ok());
  auto diagnostics = fesa::validateDomain(parsed.domain);
  FESA_CHECK(!fesa::hasError(diagnostics));
  return parsed.domain;
}

const fesa::Diagnostic* findDiagnostic(const std::vector<fesa::Diagnostic>& diagnostics, const std::string& code) {
  for (const auto& diagnostic : diagnostics) {
    if (diagnostic.code == code) {
      return &diagnostic;
    }
  }
  return nullptr;
}

std::size_t diagnosticCount(const std::vector<fesa::Diagnostic>& diagnostics, const std::string& code) {
  std::size_t count = 0;
  for (const auto& diagnostic : diagnostics) {
    if (diagnostic.code == code) {
      ++count;
    }
  }
  return count;
}

}  // namespace

FESA_TEST(core_types_and_dof_mapping_are_stable) {
  FESA_CHECK(sizeof(fesa::Real) == 8);
  FESA_CHECK(sizeof(fesa::GlobalId) == 8);
  FESA_CHECK(sizeof(fesa::LocalIndex) == 8);
  FESA_CHECK(sizeof(fesa::EquationId) == 8);
  FESA_CHECK(sizeof(fesa::SparseIndex) == 8);
  FESA_CHECK(std::numeric_limits<fesa::GlobalId>::is_signed);
  FESA_CHECK(std::numeric_limits<fesa::LocalIndex>::is_signed);
  FESA_CHECK(std::numeric_limits<fesa::EquationId>::is_signed);
  FESA_CHECK(std::numeric_limits<fesa::SparseIndex>::is_signed);

  const auto dofs = fesa::allDofs();
  FESA_CHECK(dofs.size() == 6);
  for (std::size_t i = 0; i < dofs.size(); ++i) {
    const int abaqus_number = static_cast<int>(i + 1);
    FESA_CHECK(fesa::abaqusDofNumber(dofs[i]) == abaqus_number);
    FESA_CHECK(fesa::dofFromAbaqus(abaqus_number).value() == dofs[i]);
    FESA_CHECK(std::string(fesa::dofLabel(dofs[i])) == fesa::displacementComponentLabels()[i]);
  }
  FESA_CHECK(!fesa::dofFromAbaqus(0).has_value());
  FESA_CHECK(!fesa::dofFromAbaqus(7).has_value());
}

FESA_TEST(parser_accepts_phase1_subset) {
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(phase1Input());
  FESA_CHECK(parsed.ok());
  FESA_CHECK(parsed.domain.nodes.size() == 4);
  FESA_CHECK(parsed.domain.elements.size() == 1);
  FESA_CHECK(parsed.domain.node_sets.at("left").node_ids.size() == 2);
  FESA_CHECK(parsed.domain.element_sets.at("eall").element_ids.size() == 1);
  FESA_CHECK(parsed.domain.materials.at("steel").elastic_modulus == 1000.0);
  FESA_CHECK(parsed.domain.shell_sections.front().thickness == 0.1);
  FESA_CHECK(parsed.domain.boundary_conditions.size() == 3);
  FESA_CHECK(parsed.domain.loads.size() == 2);
}

FESA_TEST(parser_accepts_repeated_and_generated_sets) {
  const std::string text = R"inp(
*Node
1, 0, 0, 0
2, 1, 0, 0
3, 1, 1, 0
4, 0, 1, 0
*Element, type=S4, elset=EALL
1, 1, 2, 3, 4
*Nset, nset=FIXED
1, 2, 2
3
*Nset, nset=FIXED, generate
3, 4, 1
*Nset, nset=LOADS, generate
2, 4, 2
*Elset, elset=EALL
1, 1
*Elset, elset=CHECK, generate
1, 5, 2
*Material, name=MAT
*Elastic
2.0D5, 0.25
*Shell Section, elset=EALL, material=MAT
0.2
*Boundary
FIXED, 1, 6
*Cload
LOADS, 3, -2.5
*Step, name=Step-A, nlgeom=NO
*Static
*End Step
)inp";
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(text);
  FESA_CHECK(parsed.ok());
  FESA_CHECK(parsed.domain.node_sets.at("fixed").node_ids == std::vector<fesa::GlobalId>({1, 2, 3, 4}));
  FESA_CHECK(parsed.domain.node_sets.at("loads").node_ids == std::vector<fesa::GlobalId>({2, 4}));
  FESA_CHECK(parsed.domain.element_sets.at("eall").element_ids == std::vector<fesa::GlobalId>({1}));
  FESA_CHECK(parsed.domain.element_sets.at("check").element_ids == std::vector<fesa::GlobalId>({1, 3, 5}));
  FESA_CHECK(parsed.domain.materials.at("mat").elastic_modulus == 2.0e5);
  FESA_CHECK(parsed.domain.steps.front().name == "Step-A");
}

FESA_TEST(parser_accepts_keyword_line_continuation) {
  const std::string text = R"inp(
*Node
1, 0, 0, 0
2, 1, 0, 0
3, 1, 1, 0
4, 0, 1, 0
*Element,
 type=S4, elset=EALL
1, 1, 2, 3, 4
*Nset,
 nset=FIXED, generate
1, 4, 3
*Elset,
 elset=EALL
1
*Material,
 name=MAT
*Elastic
2.0e5, 0.25
*Shell Section,
 elset=EALL, material=MAT
0.2
*Boundary
FIXED, 1, 6
*Cload
2, 3, -1.0
*Step,
 name=Step-1
*Static
*End Step
)inp";
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(text);
  FESA_CHECK(parsed.ok());
  FESA_CHECK(parsed.domain.elements.at(1).source_elset == "EALL");
  FESA_CHECK(parsed.domain.node_sets.at("fixed").node_ids == std::vector<fesa::GlobalId>({1, 4}));
  FESA_CHECK(parsed.domain.materials.count("mat") == 1);
  FESA_CHECK(parsed.domain.shell_sections.front().material == "MAT");
}

FESA_TEST(parser_rejects_unsupported_features) {
  const std::string text = R"inp(
*Part, name=P1
*Assembly, name=A1
*Instance, name=I1, part=P1
*Include, input=other.inp
*Node
1, 0, 0, 0
*Element, type=S4R
1, 1, 2, 3, 4
*Density
7850
*Step, nlgeom=YES
*End Step
)inp";
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(text);
  FESA_CHECK(!parsed.ok());
  FESA_CHECK(diagnosticCount(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-KEYWORD") >= 4);
  FESA_CHECK(fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-KEYWORD"));
  FESA_CHECK(fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-ELEMENT"));
  FESA_CHECK(fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-NLGEOM"));
}

FESA_TEST(parser_rejects_unsupported_keyword_parameters_and_modes) {
  const std::string text = R"inp(
*Node, input=nodes.csv
1, 0, 0, 0
2, 1, 0, 0
3, 1, 1, 0
4, 0, 1, 0
*Element, type=S4, elset=EALL, orientation=OR1
1, 1, 2, 3, 4
*Nset, nset=FIXED, unsorted
1, 4
*Material, name=MAT, description=bad
*Elastic, type=ENGINEERING CONSTANTS
2.0e5, 0.25
*Shell Section, elset=EALL, material=MAT, offset=SPOS
0.2, 5
*Boundary, op=NEW
FIXED, 1, 6
*Cload, amplitude=A1
2, 3, -1.0
*Step, name=Step-1, inc=100
*Static, stabilize
*End Step
)inp";
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(text, "unsupported_modes.inp");
  FESA_CHECK(!parsed.ok());
  FESA_CHECK(diagnosticCount(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-PARAMETER") >= 8);
  FESA_CHECK(fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-SHELL-SECTION-UNSUPPORTED"));
}

FESA_TEST(parser_diagnostics_include_file_line_and_keyword) {
  const std::string text = R"inp(
*Node
1, bad, 0, 0
*Boundary
FIXED, 7, 7
)inp";
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(text, "malformed.inp");
  FESA_CHECK(!parsed.ok());

  const fesa::Diagnostic* node = findDiagnostic(parsed.diagnostics, "FESA-PARSE-NODE-NUMERIC");
  FESA_CHECK(node != nullptr);
  FESA_CHECK(node->source.file == "malformed.inp");
  FESA_CHECK(node->source.line == 3);
  FESA_CHECK(node->source.keyword == "node");

  const fesa::Diagnostic* boundary = findDiagnostic(parsed.diagnostics, "FESA-PARSE-BOUNDARY-DOF");
  FESA_CHECK(boundary != nullptr);
  FESA_CHECK(boundary->source.file == "malformed.inp");
  FESA_CHECK(boundary->source.line == 5);
  FESA_CHECK(boundary->source.keyword == "boundary");
}

FESA_TEST(quad01_reference_input_remains_unsupported) {
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseFile(sourceRoot() + "/references/quad_01.inp");
  FESA_CHECK(!parsed.ok());
  FESA_CHECK(fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-KEYWORD") ||
             fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-ELEMENT"));
}

FESA_TEST(quad02_original_reference_input_remains_unsupported) {
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseFile(sourceRoot() + "/references/quad_02.inp");
  FESA_CHECK(!parsed.ok());
  FESA_CHECK(fesa::containsDiagnostic(parsed.diagnostics, "FESA-PARSE-UNSUPPORTED-KEYWORD"));
}

FESA_TEST(quad02_phase1_normalized_input_uses_supported_subset) {
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseFile(sourceRoot() + "/references/quad_02_phase1.inp");
  FESA_CHECK(parsed.ok());
  FESA_CHECK(parsed.domain.nodes.size() == 121);
  FESA_CHECK(parsed.domain.elements.size() == 100);
  FESA_CHECK(parsed.domain.node_sets.at("fixed_boundary").node_ids.size() == 40);
  FESA_CHECK(parsed.domain.node_sets.at("load_node").node_ids.size() == 1);
  FESA_CHECK(parsed.domain.element_sets.at("all_elements").element_ids.size() == 100);
  FESA_CHECK(parsed.domain.materials.at("material_1").elastic_modulus == 7.0e10);
  FESA_CHECK(parsed.domain.shell_sections.front().thickness == 1.0);
}

FESA_TEST(domain_validation_reports_missing_property_and_targets) {
  fesa::Domain domain;
  domain.nodes[1] = {1, {0, 0, 0}};
  domain.nodes[2] = {2, {1, 0, 0}};
  domain.nodes[3] = {3, {1, 1, 0}};
  domain.nodes[4] = {4, {0, 1, 0}};
  domain.elements[1] = {1, fesa::ElementType::MITC4, {1, 2, 3, 99}, ""};
  domain.shell_sections.push_back({"MISSING_ELSET", "MISSING_MAT", 0.1});
  domain.loads.push_back({"MISSING", 3, 1.0});
  auto diagnostics = fesa::validateDomain(domain);
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-ELEMENT-MISSING-NODE"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-MISSING-PROPERTY"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-MISSING-ELSET"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-MISSING-MATERIAL"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-MISSING-NSET"));
  for (const auto& diagnostic : diagnostics) {
    FESA_CHECK(!diagnostic.code.empty());
    FESA_CHECK(!diagnostic.message.empty());
    FESA_CHECK(!diagnostic.source.keyword.empty());
  }
}

FESA_TEST(dof_manager_owns_equation_numbering_and_reconstruction) {
  auto domain = parsedPhase1Domain();
  fesa::DofManager dofs(domain);
  FESA_CHECK(dofs.fullDofCount() == 24);
  FESA_CHECK(dofs.freeDofCount() == 2);
  FESA_CHECK(dofs.isConstrained(1, fesa::Dof::UX));
  FESA_CHECK(dofs.equation(2, fesa::Dof::UZ) == 0);
  FESA_CHECK(dofs.equation(3, fesa::Dof::UZ) == 1);
  auto full = dofs.reconstructFullVector({-0.1, -0.2});
  FESA_CHECK_NEAR(full[static_cast<std::size_t>(dofs.fullIndex(2, fesa::Dof::UZ))], -0.1, 1.0e-15);
  FESA_CHECK_NEAR(full[static_cast<std::size_t>(dofs.fullIndex(1, fesa::Dof::UX))], 0.0, 1.0e-15);
}

FESA_TEST(gaussian_solver_solves_and_diagnoses_singular_systems) {
  fesa::DenseMatrix a(2, 2);
  a(0, 0) = 2.0;
  a(0, 1) = 1.0;
  a(1, 0) = 1.0;
  a(1, 1) = 3.0;
  fesa::GaussianEliminationSolver solver;
  auto solved = solver.solve(a, {1.0, 2.0});
  FESA_CHECK(solved.ok());
  FESA_CHECK_NEAR(solved.x[0], 0.2, 1.0e-12);
  FESA_CHECK_NEAR(solved.x[1], 0.6, 1.0e-12);

  fesa::DenseMatrix singular(2, 2);
  singular(0, 0) = 1.0;
  singular(0, 1) = 2.0;
  singular(1, 0) = 2.0;
  singular(1, 1) = 4.0;
  auto failed = solver.solve(singular, {1.0, 2.0});
  FESA_CHECK(!failed.ok());
  FESA_CHECK(fesa::containsDiagnostic(failed.diagnostics, "FESA-SINGULAR-SOLVER"));
}

FESA_TEST(results_writer_uses_step_frame_fields_for_u_and_rf) {
  auto domain = parsedPhase1Domain();
  fesa::DofManager dofs(domain);
  std::vector<fesa::Real> u(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
  std::vector<fesa::Real> rf(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
  u[static_cast<std::size_t>(dofs.fullIndex(2, fesa::Dof::UZ))] = -0.1;
  rf[static_cast<std::size_t>(dofs.fullIndex(1, fesa::Dof::UZ))] = 1.0;
  fesa::InMemoryResultsWriter writer;
  writer.writeLinearStatic(domain, dofs, u, rf);
  const auto& result = writer.result();
  FESA_CHECK(result.schema_name == "FESA_RESULTS");
  FESA_CHECK(result.steps.size() == 1);
  FESA_CHECK(result.steps[0].frames[0].field_outputs.count("U") == 1);
  FESA_CHECK(result.steps[0].frames[0].field_outputs.count("RF") == 1);
  FESA_CHECK(result.steps[0].frames[0].field_outputs.at("U").component_labels[2] == "UZ");
  FESA_CHECK(result.steps[0].frames[0].field_outputs.at("RF").component_labels[2] == "RFZ");
}

FESA_TEST(displacement_csv_loader_accepts_quad01_format) {
  auto table = fesa::loadDisplacementCsv(sourceRoot() + "/references/quad_01_displacements.csv");
  FESA_CHECK(!fesa::hasError(table.diagnostics));
  FESA_CHECK(table.rows.size() == 121);
  FESA_CHECK(table.rows.count(1) == 1);
}

FESA_TEST(displacement_csv_loader_accepts_quad02_format) {
  auto table = fesa::loadDisplacementCsv(sourceRoot() + "/references/quad_02_displacements.csv");
  FESA_CHECK(!fesa::hasError(table.diagnostics));
  FESA_CHECK(table.rows.size() == 121);
  FESA_CHECK(table.rows.count(2) == 1);
  FESA_CHECK(table.rows.at(2).values[2] < 0.0);
}

FESA_TEST(displacement_comparator_matches_by_node_id_not_row_order) {
  fesa::FieldOutput actual;
  actual.name = "U";
  actual.entity_ids = {2, 1};
  actual.component_labels = fesa::displacementComponentLabels();
  actual.values = {{{2, 0, 0, 0, 0, 0}}, {{1, 0, 0, 0, 0, 0}}};
  fesa::CsvDisplacementTable expected;
  expected.rows[1] = {1, {1, 0, 0, 0, 0, 0}};
  expected.rows[2] = {2, {2, 0, 0, 0, 0, 0}};
  auto compared = fesa::compareDisplacements(actual, expected, {1.0e-12, 1.0e-12, 1.0});
  FESA_CHECK(compared.pass);
}

FESA_TEST(mitc4_shape_functions_and_stiffness_baseline) {
  auto shape = fesa::shapeFunctions(0.25, -0.5);
  const fesa::Real sum = shape.n[0] + shape.n[1] + shape.n[2] + shape.n[3];
  FESA_CHECK_NEAR(sum, 1.0, 1.0e-15);

  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  fesa::MITC4ElementKernel kernel;
  auto k = kernel.stiffness(coords, 1000.0, 0.3, 0.1);
  FESA_CHECK(k.rows() == 24);
  FESA_CHECK(k.cols() == 24);
  for (fesa::LocalIndex i = 0; i < 24; ++i) {
    for (fesa::LocalIndex j = 0; j < 24; ++j) {
      FESA_CHECK_NEAR(k(i, j), k(j, i), 1.0e-8);
    }
  }

  std::vector<fesa::Real> uniform_translation(24, 0.0);
  for (int node = 0; node < 4; ++node) {
    uniform_translation[static_cast<std::size_t>(6 * node + 0)] = 1.0;
  }
  auto internal = k.multiply(uniform_translation);
  fesa::Real norm = 0.0;
  for (auto value : internal) {
    norm += std::fabs(value);
  }
  FESA_CHECK(norm < 1.0e-8);
}

FESA_TEST(linear_static_analysis_solves_u_and_recovers_full_vector_rf) {
  auto domain = parsedPhase1Domain();
  fesa::LinearStaticAnalysis analysis;
  auto result = analysis.run(domain);
  FESA_CHECK(result.ok());
  FESA_CHECK(result.state.converged);
  FESA_CHECK(result.result_file.steps.size() == 1);
  const auto& frame = result.result_file.steps[0].frames[0];
  FESA_CHECK(frame.field_outputs.count("U") == 1);
  FESA_CHECK(frame.field_outputs.count("RF") == 1);
  fesa::Real total_rf_z = 0.0;
  fesa::DofManager dofs(domain);
  for (auto node_id : dofs.nodeIds()) {
    total_rf_z += result.state.reaction_full[static_cast<std::size_t>(dofs.fullIndex(node_id, fesa::Dof::UZ))];
  }
  FESA_CHECK_NEAR(total_rf_z, 2.0, 1.0e-8);
}

int main() {
  int failed = 0;
  for (const auto& test : registry()) {
    try {
      test.fn();
      std::cout << "[PASS] " << test.name << '\n';
    } catch (const std::exception& error) {
      ++failed;
      std::cerr << "[FAIL] " << test.name << ": " << error.what() << '\n';
    }
  }
  if (failed != 0) {
    std::cerr << failed << " test(s) failed\n";
    return 1;
  }
  std::cout << registry().size() << " test(s) passed\n";
  return 0;
}