#include "fesa/Analysis/Analysis.hpp"
#include "fesa/Assembly/Assembly.hpp"
#include "fesa/Boundary/Boundary.hpp"
#include "fesa/Core/Core.hpp"
#include "fesa/Element/Element.hpp"
#include "fesa/IO/IO.hpp"
#include "fesa/Load/Load.hpp"
#include "fesa/Math/Math.hpp"
#include "fesa/Material/Material.hpp"
#include "fesa/Property/Property.hpp"
#include "fesa/Results/Results.hpp"
#include "fesa/Util/Util.hpp"
#include "fesa/fesa.hpp"

#include <array>
#include <cmath>
#include <functional>
#include <iostream>
#include <type_traits>
#include <stdexcept>
#include <string>
#include <string_view>
#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;
}

void checkVecNear(const fesa::Vec3& actual, const fesa::Vec3& expected, fesa::Real tolerance) {
  FESA_CHECK_NEAR(actual.x, expected.x, tolerance);
  FESA_CHECK_NEAR(actual.y, expected.y, tolerance);
  FESA_CHECK_NEAR(actual.z, expected.z, tolerance);
}

void checkRightHandedOrthonormal(const fesa::Vec3& e1, const fesa::Vec3& e2, const fesa::Vec3& e3) {
  FESA_CHECK_NEAR(fesa::norm(e1), 1.0, 1.0e-14);
  FESA_CHECK_NEAR(fesa::norm(e2), 1.0, 1.0e-14);
  FESA_CHECK_NEAR(fesa::norm(e3), 1.0, 1.0e-14);
  FESA_CHECK_NEAR(fesa::dot(e1, e2), 0.0, 1.0e-14);
  FESA_CHECK_NEAR(fesa::dot(e1, e3), 0.0, 1.0e-14);
  FESA_CHECK_NEAR(fesa::dot(e2, e3), 0.0, 1.0e-14);
  FESA_CHECK_NEAR(fesa::dot(fesa::cross(e1, e2), e3), 1.0, 1.0e-14);
}

std::array<fesa::Real, 24> zeroElementDofs() {
  std::array<fesa::Real, 24> values{};
  values.fill(0.0);
  return values;
}

std::vector<fesa::Real> zeroElementVector() {
  return std::vector<fesa::Real>(24, 0.0);
}

fesa::Real quadraticEnergy(const fesa::DenseMatrix& stiffness, const std::vector<fesa::Real>& values) {
  const auto internal = stiffness.multiply(values);
  fesa::Real energy = 0.0;
  for (std::size_t i = 0; i < values.size(); ++i) {
    energy += values[i] * internal[i];
  }
  return energy;
}

std::array<fesa::Vec3, 4> flatUnitSquareCoordinates() {
  return {fesa::Vec3{0.0, 0.0, 0.0},
          fesa::Vec3{1.0, 0.0, 0.0},
          fesa::Vec3{1.0, 1.0, 0.0},
          fesa::Vec3{0.0, 1.0, 0.0}};
}

using MITC4DofField = std::function<fesa::Vec3(const fesa::Vec3&)>;

fesa::MITC4ElementDofVector elementDofsFromFields(const std::array<fesa::Vec3, 4>& coords,
                                                  MITC4DofField translation,
                                                  MITC4DofField rotation) {
  fesa::MITC4ElementDofVector values{};
  for (std::size_t node = 0; node < coords.size(); ++node) {
    const auto displacement = translation(coords[node]);
    const auto rotational_dof = rotation(coords[node]);
    values[node * 6 + 0] = displacement.x;
    values[node * 6 + 1] = displacement.y;
    values[node * 6 + 2] = displacement.z;
    values[node * 6 + 3] = rotational_dof.x;
    values[node * 6 + 4] = rotational_dof.y;
    values[node * 6 + 5] = rotational_dof.z;
  }
  return values;
}

std::vector<fesa::Real> toVector(const fesa::MITC4ElementDofVector& values) {
  return {values.begin(), values.end()};
}

std::string readTextFile(const std::string& path) {
  std::ifstream input(path);
  std::ostringstream buffer;
  buffer << input.rdbuf();
  FESA_CHECK(input.good() || !buffer.str().empty());
  return buffer.str();
}

void checkComparisonPass(const fesa::ComparisonResult& comparison) {
  if (!comparison.pass) {
    throw std::runtime_error("reference comparison failed: max_abs_error=" +
                             std::to_string(comparison.max_abs_error) +
                             ", max_rel_error=" + std::to_string(comparison.max_rel_error) +
                             ", diagnostics=" + std::to_string(comparison.diagnostics.size()));
  }
}

fesa::MITC4StrainVector localStrainAt(const fesa::MITC4Geometry& geometry,
                                      const fesa::MITC4ElementDofVector& values,
                                      fesa::Real xi,
                                      fesa::Real eta,
                                      fesa::Real zeta) {
  const auto rows = fesa::mitc4TiedCovariantStrainRows(geometry, xi, eta, zeta);
  FESA_CHECK(rows.ok());
  const auto convected_strain = fesa::evaluateMITC4StrainRows(rows, values);
  FESA_CHECK(convected_strain.ok());
  const auto basis = fesa::computeMITC4IntegrationBasis(geometry, xi, eta, zeta);
  FESA_CHECK(basis.ok());
  const auto transform = fesa::mitc4CovariantToLocalStrainTransform(basis);
  FESA_CHECK(transform.ok());
  return fesa::multiplyMITC4MaterialMatrix(transform.matrix, convected_strain.values);
}

fesa::Real singleElementCantileverTipUz(fesa::Real thickness) {
  fesa::Domain domain;
  domain.nodes[1] = {1, {0.0, 0.0, 0.0}};
  domain.nodes[2] = {2, {4.0, 0.0, 0.0}};
  domain.nodes[3] = {3, {4.0, 1.0, 0.0}};
  domain.nodes[4] = {4, {0.0, 1.0, 0.0}};
  domain.elements[1] = {1, fesa::ElementType::MITC4, {1, 2, 3, 4}, "SHELLS"};
  domain.node_sets["clamped"] = {"CLAMPED", {1, 4}};
  domain.element_sets["shells"] = {"SHELLS", {1}};
  domain.materials["shell_steel"] = {"SHELL_STEEL", 210000.0, 0.30};
  domain.shell_sections.push_back({"SHELLS", "SHELL_STEEL", thickness});
  domain.boundary_conditions.push_back({"CLAMPED", 1, 6, 0.0});
  domain.loads.push_back({"2", 3, -0.5});
  domain.loads.push_back({"3", 3, -0.5});

  fesa::LinearStaticAnalysis analysis;
  const auto result = analysis.run(domain);
  FESA_CHECK(result.ok());
  fesa::DofManager dofs(domain);
  const auto node2_uz = result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(2, fesa::Dof::UZ))];
  const auto node3_uz = result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(3, fesa::Dof::UZ))];
  return 0.5 * (node2_uz + node3_uz);
}

class RecordingSolver final : public fesa::LinearSolver {
 public:
  explicit RecordingSolver(std::vector<fesa::Real> solution) : solution_(std::move(solution)) {}

  fesa::SolveResult solve(fesa::DenseMatrix a, std::vector<fesa::Real> b) const override {
    called = true;
    captured_a = std::move(a);
    captured_b = std::move(b);
    return {solution_, {}};
  }

  mutable bool called = false;
  mutable fesa::DenseMatrix captured_a;
  mutable std::vector<fesa::Real> captured_b;

 private:
  std::vector<fesa::Real> solution_;
};

class FailingSolver final : public fesa::LinearSolver {
 public:
  fesa::SolveResult solve(fesa::DenseMatrix, std::vector<fesa::Real>) const override {
    return {{}, {fesa::makeDiagnostic(fesa::Severity::Error, "FESA-SINGULAR-SOLVER",
                                      "Injected reduced system singularity", "solver")}};
  }
};

fesa::Domain singleElementValidationDomain() {
  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, 4}, "EALL"};
  domain.element_sets["eall"] = {"EALL", {1}};
  domain.node_sets["all_nodes"] = {"ALL_NODES", {1, 2, 3, 4}};
  domain.materials["mat"] = {"MAT", 1000.0, 0.3};
  domain.shell_sections.push_back({"EALL", "MAT", 0.1});
  domain.loads.push_back({"2", 3, -1.0});
  return domain;
}

fesa::Domain noncontiguousDofDomain() {
  fesa::Domain domain;
  domain.nodes[30] = {30, {1, 1, 0}};
  domain.nodes[10] = {10, {0, 0, 0}};
  domain.nodes[40] = {40, {0, 1, 0}};
  domain.nodes[20] = {20, {1, 0, 0}};
  domain.elements[5] = {5, fesa::ElementType::MITC4, {10, 20, 30, 40}, "EALL"};
  domain.element_sets["eall"] = {"EALL", {5}};
  domain.node_sets["clamped"] = {"CLAMPED", {20}};
  domain.materials["mat"] = {"MAT", 1000.0, 0.3};
  domain.shell_sections.push_back({"EALL", "MAT", 0.1});
  domain.boundary_conditions.push_back({"CLAMPED", 1, 6, 0.0});
  domain.boundary_conditions.push_back({"30", 1, 2, 0.0});
  return domain;
}

}  // 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(module_scaffold_headers_are_include_compatible_with_umbrella) {
  const auto modules = fesa::architectureModules();
  FESA_CHECK(modules.size() == 12);
  FESA_CHECK(modules[0] == std::string_view("Analysis"));
  FESA_CHECK(modules[1] == std::string_view("Assembly"));
  FESA_CHECK(modules[2] == std::string_view("Boundary"));
  FESA_CHECK(modules[3] == std::string_view("Core"));
  FESA_CHECK(modules[4] == std::string_view("Element"));
  FESA_CHECK(modules[5] == std::string_view("IO"));
  FESA_CHECK(modules[6] == std::string_view("Load"));
  FESA_CHECK(modules[7] == std::string_view("Math"));
  FESA_CHECK(modules[8] == std::string_view("Material"));
  FESA_CHECK(modules[9] == std::string_view("Property"));
  FESA_CHECK(modules[10] == std::string_view("Results"));
  FESA_CHECK(modules[11] == std::string_view("Util"));
  FESA_CHECK(fesa::umbrellaFacadeHeader() == std::string_view("fesa/fesa.hpp"));
  FESA_CHECK(fesa::dofFromAbaqus(1).value() == fesa::Dof::UX);
}

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(analysis_model_activates_single_linear_static_step) {
  auto domain = parsedPhase1Domain();
  const auto node_count = domain.nodes.size();
  const auto element_count = domain.elements.size();
  const auto model = fesa::buildLinearStaticAnalysisModel(domain);

  FESA_CHECK(model.ok());
  FESA_CHECK(model.step.name == "Step-1");
  FESA_CHECK(model.step.analysis_type == "linear_static");
  FESA_CHECK(model.active_element_ids == std::vector<fesa::GlobalId>({1}));
  FESA_CHECK(model.active_boundary_condition_indices == std::vector<std::size_t>({0, 1, 2}));
  FESA_CHECK(model.active_load_indices == std::vector<std::size_t>({0, 1}));
  FESA_CHECK(model.active_shell_section_indices == std::vector<std::size_t>({0}));
  FESA_CHECK(model.active_material_keys == std::vector<std::string>({"steel"}));
  FESA_CHECK(domain.nodes.size() == node_count);
  FESA_CHECK(domain.elements.size() == element_count);
}

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(analysis_model_rejects_multiple_steps_for_phase1_execution) {
  const std::string text = phase1Input() + R"inp(
*Step, name=Step-2
*Static
*End Step
)inp";
  fesa::AbaqusInputParser parser;
  auto parsed = parser.parseString(text);
  FESA_CHECK(parsed.ok());
  FESA_CHECK(parsed.domain.steps.size() == 2);
  const auto model = fesa::buildLinearStaticAnalysisModel(parsed.domain);
  FESA_CHECK(!model.ok());
  FESA_CHECK(fesa::containsDiagnostic(model.diagnostics, "FESA-ANALYSIS-MULTIPLE-STEPS"));
}

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(domain_validation_reports_missing_sets_and_set_members) {
  auto domain = singleElementValidationDomain();
  domain.shell_sections.clear();
  domain.shell_sections.push_back({"MISSING_ELSET", "MISSING_MAT", 0.1});
  domain.node_sets["bad_nodes"] = {"BAD_NODES", {1, 99}};
  domain.element_sets["bad_elements"] = {"BAD_ELEMENTS", {1, 77}};
  domain.boundary_conditions.push_back({"MISSING_BOUNDARY_SET", 1, 6, 0.0});
  domain.boundary_conditions.push_back({"BAD_NODES", 1, 1, 0.0});
  domain.loads.push_back({"MISSING_LOAD_SET", 3, 1.0});
  domain.loads.push_back({"BAD_NODES", 3, 1.0});

  auto diagnostics = fesa::validateDomain(domain);
  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"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-NSET-MISSING-NODE"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-ELSET-MISSING-ELEMENT"));

  const auto* missing_set = findDiagnostic(diagnostics, "FESA-VALIDATION-MISSING-NSET");
  FESA_CHECK(missing_set != nullptr);
  FESA_CHECK(missing_set->message.find("MISSING_BOUNDARY_SET") != std::string::npos ||
             missing_set->message.find("MISSING_LOAD_SET") != std::string::npos);

  const auto* missing_node = findDiagnostic(diagnostics, "FESA-VALIDATION-NSET-MISSING-NODE");
  FESA_CHECK(missing_node != nullptr);
  FESA_CHECK(missing_node->message.find("BAD_NODES") != std::string::npos);
  FESA_CHECK(missing_node->message.find("99") != std::string::npos);
}

FESA_TEST(domain_validation_reports_nonpositive_thickness_and_invalid_dofs) {
  auto domain = singleElementValidationDomain();
  domain.shell_sections.front().thickness = 0.0;
  domain.boundary_conditions.push_back({"ALL_NODES", 0, 1, 0.0});
  domain.loads.push_back({"2", 7, 1.0});

  auto diagnostics = fesa::validateDomain(domain);
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-NONPOSITIVE-THICKNESS"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-BOUNDARY-DOF"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-VALIDATION-CLOAD-DOF"));

  const auto* thickness = findDiagnostic(diagnostics, "FESA-VALIDATION-NONPOSITIVE-THICKNESS");
  FESA_CHECK(thickness != nullptr);
  FESA_CHECK(thickness->message.find("EALL") != std::string::npos);

  const auto* load_dof = findDiagnostic(diagnostics, "FESA-VALIDATION-CLOAD-DOF");
  FESA_CHECK(load_dof != nullptr);
  FESA_CHECK(load_dof->message.find("DOF 7") != std::string::npos);
}

FESA_TEST(domain_validation_reports_no_active_elements_and_missing_load) {
  fesa::Domain domain;
  domain.nodes[1] = {1, {0, 0, 0}};
  auto diagnostics = fesa::validateDomain(domain);
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-SINGULAR-NO-ACTIVE-ELEMENTS"));
  FESA_CHECK(fesa::containsDiagnostic(diagnostics, "FESA-SINGULAR-NO-NONZERO-LOAD"));
  for (const auto& diagnostic : diagnostics) {
    FESA_CHECK(diagnostic.code.find("MESH") == std::string::npos);
  }
}

FESA_TEST(domain_validation_reports_no_free_dofs) {
  auto domain = singleElementValidationDomain();
  domain.boundary_conditions.push_back({"ALL_NODES", 1, 6, 0.0});

  auto diagnostics = fesa::validateDomain(domain);
  const auto* no_free = findDiagnostic(diagnostics, "FESA-SINGULAR-NO-FREE-DOFS");
  FESA_CHECK(no_free != nullptr);
  FESA_CHECK(no_free->source.keyword == "dof");
  FESA_CHECK(no_free->message.find("No free DOFs") != std::string::npos);
}

FESA_TEST(domain_validation_reports_free_untouched_dofs_for_isolated_nodes) {
  auto domain = singleElementValidationDomain();
  domain.nodes[99] = {99, {10, 0, 0}};
  domain.boundary_conditions.push_back({"ALL_NODES", 1, 6, 0.0});

  auto diagnostics = fesa::validateDomain(domain);
  const auto* untouched = findDiagnostic(diagnostics, "FESA-SINGULAR-DOF-UNTOUCHED");
  FESA_CHECK(untouched != nullptr);
  FESA_CHECK(untouched->source.keyword == "dof");
  FESA_CHECK(untouched->message.find("Node 99") != std::string::npos);
  FESA_CHECK(untouched->message.find("UX") != std::string::npos);
}

FESA_TEST(domain_validation_reports_weak_drilling_dof_smoke) {
  auto domain = singleElementValidationDomain();
  domain.boundary_conditions.push_back({"ALL_NODES", 1, 5, 0.0});

  auto diagnostics = fesa::validateDomain(domain);
  const auto* weak_drilling = findDiagnostic(diagnostics, "FESA-SINGULAR-WEAK-DRILLING-DOF");
  FESA_CHECK(weak_drilling != nullptr);
  FESA_CHECK(weak_drilling->severity == fesa::Severity::Warning);
  FESA_CHECK(weak_drilling->source.keyword == "dof");
  FESA_CHECK(weak_drilling->message.find("RZ") != std::string::npos);
}

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(dof_manager_preserves_global_order_for_noncontiguous_nodes) {
  auto domain = noncontiguousDofDomain();
  fesa::DofManager dofs(domain);

  FESA_CHECK(dofs.nodeIds() == std::vector<fesa::GlobalId>({10, 20, 30, 40}));
  FESA_CHECK(dofs.fullDofCount() == 24);
  for (std::size_t node_offset = 0; node_offset < dofs.nodeIds().size(); ++node_offset) {
    const fesa::GlobalId node_id = dofs.nodeIds()[node_offset];
    for (fesa::Dof dof : fesa::allDofs()) {
      const fesa::LocalIndex expected = static_cast<fesa::LocalIndex>(6 * node_offset + fesa::dofIndex(dof));
      FESA_CHECK(dofs.fullIndex(node_id, dof) == expected);
      const auto address = dofs.fullDof(expected);
      FESA_CHECK(address.node_id == node_id);
      FESA_CHECK(address.dof == dof);
    }
  }
}

FESA_TEST(dof_manager_partitions_constraints_and_equations_deterministically) {
  auto domain = noncontiguousDofDomain();
  fesa::DofManager dofs(domain);

  FESA_CHECK(dofs.constrainedDofCount() == 8);
  FESA_CHECK(dofs.freeDofCount() == 16);
  FESA_CHECK(dofs.constrainedFullIndices() == std::vector<fesa::LocalIndex>({6, 7, 8, 9, 10, 11, 12, 13}));
  FESA_CHECK(dofs.freeFullIndices() ==
             std::vector<fesa::LocalIndex>({0, 1, 2, 3, 4, 5, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}));
  FESA_CHECK(dofs.equation(10, fesa::Dof::UX) == 0);
  FESA_CHECK(dofs.equation(20, fesa::Dof::RZ) == -1);
  FESA_CHECK(dofs.equation(30, fesa::Dof::UX) == -1);
  FESA_CHECK(dofs.equation(30, fesa::Dof::UZ) == 6);
  FESA_CHECK(dofs.equation(40, fesa::Dof::RZ) == 15);
}

FESA_TEST(dof_manager_reduces_and_reconstructs_full_vectors) {
  auto domain = noncontiguousDofDomain();
  fesa::DofManager dofs(domain);

  std::vector<fesa::Real> full(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
  for (std::size_t i = 0; i < full.size(); ++i) {
    full[i] = static_cast<fesa::Real>(100 + i);
  }

  auto reduced = dofs.reduceFullVector(full);
  FESA_CHECK(reduced.size() == static_cast<std::size_t>(dofs.freeDofCount()));
  FESA_CHECK_NEAR(reduced[0], full[0], 1.0e-15);
  FESA_CHECK_NEAR(reduced[6], full[14], 1.0e-15);
  FESA_CHECK_NEAR(reduced.back(), full[23], 1.0e-15);

  auto reconstructed = dofs.reconstructFullVector(reduced);
  for (fesa::LocalIndex full_index = 0; full_index < dofs.fullDofCount(); ++full_index) {
    const auto address = dofs.fullDof(full_index);
    if (dofs.isConstrained(address.node_id, address.dof)) {
      FESA_CHECK_NEAR(reconstructed[static_cast<std::size_t>(full_index)], 0.0, 1.0e-15);
    } else {
      FESA_CHECK_NEAR(reconstructed[static_cast<std::size_t>(full_index)], full[static_cast<std::size_t>(full_index)], 1.0e-15);
    }
  }
}

FESA_TEST(dof_manager_provides_element_sparse_connectivity_inputs) {
  auto domain = parsedPhase1Domain();
  fesa::DofManager dofs(domain);
  const auto& element = domain.elements.at(1);

  auto full_indices = dofs.elementFullDofIndices(element);
  auto equation_ids = dofs.elementEquationIds(element);

  for (fesa::LocalIndex i = 0; i < 24; ++i) {
    FESA_CHECK(full_indices[static_cast<std::size_t>(i)] == i);
  }
  for (fesa::LocalIndex i = 0; i < 6; ++i) {
    FESA_CHECK(equation_ids[static_cast<std::size_t>(i)] == -1);
    FESA_CHECK(equation_ids[static_cast<std::size_t>(18 + i)] == -1);
  }
  FESA_CHECK(equation_ids[static_cast<std::size_t>(6 + fesa::dofIndex(fesa::Dof::UZ))] == 0);
  FESA_CHECK(equation_ids[static_cast<std::size_t>(12 + fesa::dofIndex(fesa::Dof::UZ))] == 1);
  FESA_CHECK(equation_ids[static_cast<std::size_t>(6 + fesa::dofIndex(fesa::Dof::UX))] == -1);
}

FESA_TEST(full_vector_reaction_recovery_uses_full_system_quantities) {
  auto domain = noncontiguousDofDomain();
  fesa::DofManager dofs(domain);
  fesa::DenseMatrix k_full(dofs.fullDofCount(), dofs.fullDofCount());
  std::vector<fesa::Real> u_full(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
  std::vector<fesa::Real> f_full(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);

  const fesa::LocalIndex support_ux = dofs.fullIndex(20, fesa::Dof::UX);
  const fesa::LocalIndex free_uz = dofs.fullIndex(30, fesa::Dof::UZ);
  k_full(support_ux, support_ux) = 10.0;
  k_full(support_ux, free_uz) = 2.0;
  k_full(free_uz, support_ux) = 2.0;
  k_full(free_uz, free_uz) = 5.0;
  u_full[static_cast<std::size_t>(free_uz)] = 3.0;
  f_full[static_cast<std::size_t>(support_ux)] = 1.0;
  f_full[static_cast<std::size_t>(free_uz)] = 7.0;

  auto reaction = fesa::recoverFullReaction(k_full, u_full, f_full);
  FESA_CHECK_NEAR(reaction[static_cast<std::size_t>(support_ux)], 5.0, 1.0e-15);
  FESA_CHECK_NEAR(reaction[static_cast<std::size_t>(free_uz)], 8.0, 1.0e-15);
}

FESA_TEST(reduced_sparse_pattern_is_deterministic_for_phase1_connectivity) {
  auto domain = parsedPhase1Domain();
  fesa::DofManager dofs(domain);
  const auto pattern = fesa::buildReducedSparsePattern(domain, dofs);

  FESA_CHECK(pattern.equation_count == dofs.freeDofCount());
  FESA_CHECK(pattern.nonzeroCount() == 4);
  FESA_CHECK(pattern.contains(0, 0));
  FESA_CHECK(pattern.contains(0, 1));
  FESA_CHECK(pattern.contains(1, 0));
  FESA_CHECK(pattern.contains(1, 1));
  FESA_CHECK(pattern.entries[0].row == 0);
  FESA_CHECK(pattern.entries[0].col == 0);
  FESA_CHECK(pattern.entries[1].row == 0);
  FESA_CHECK(pattern.entries[1].col == 1);
  FESA_CHECK(pattern.entries[2].row == 1);
  FESA_CHECK(pattern.entries[2].col == 0);
  FESA_CHECK(pattern.entries[3].row == 1);
  FESA_CHECK(pattern.entries[3].col == 1);
}

FESA_TEST(assembly_projection_uses_dof_manager_free_indices) {
  auto domain = parsedPhase1Domain();
  fesa::DofManager dofs(domain);
  fesa::AssemblyResult assembly;
  assembly.k_full = fesa::DenseMatrix(dofs.fullDofCount(), dofs.fullDofCount());
  assembly.f_full = std::vector<fesa::Real>(static_cast<std::size_t>(dofs.fullDofCount()), 0.0);
  assembly.reduced_pattern = fesa::buildReducedSparsePattern(domain, dofs);

  const auto node2_uz = dofs.fullIndex(2, fesa::Dof::UZ);
  const auto node3_uz = dofs.fullIndex(3, fesa::Dof::UZ);
  const auto support_uz = dofs.fullIndex(1, fesa::Dof::UZ);
  assembly.k_full(node2_uz, node2_uz) = 3.0;
  assembly.k_full(node2_uz, node3_uz) = 1.0;
  assembly.k_full(node3_uz, node2_uz) = 1.0;
  assembly.k_full(node3_uz, node3_uz) = 2.0;
  assembly.k_full(support_uz, node2_uz) = 7.0;
  assembly.f_full[static_cast<std::size_t>(node2_uz)] = -1.0;
  assembly.f_full[static_cast<std::size_t>(node3_uz)] = -2.0;

  const auto reduced = fesa::projectToReducedSystem(assembly, dofs);
  FESA_CHECK(reduced.ok());
  FESA_CHECK(reduced.k.rows() == 2);
  FESA_CHECK(reduced.k.cols() == 2);
  FESA_CHECK(reduced.free_full_indices == dofs.freeFullIndices());
  FESA_CHECK_NEAR(reduced.k(0, 0), 3.0, 1.0e-15);
  FESA_CHECK_NEAR(reduced.k(0, 1), 1.0, 1.0e-15);
  FESA_CHECK_NEAR(reduced.k(1, 0), 1.0, 1.0e-15);
  FESA_CHECK_NEAR(reduced.k(1, 1), 2.0, 1.0e-15);
  FESA_CHECK_NEAR(reduced.f[0], -1.0, 1.0e-15);
  FESA_CHECK_NEAR(reduced.f[1], -2.0, 1.0e-15);

  fesa::GaussianEliminationSolver solver;
  const auto solved = solver.solve(reduced.k, reduced.f);
  FESA_CHECK(solved.ok());
  FESA_CHECK_NEAR(solved.x[0], 0.0, 1.0e-12);
  FESA_CHECK_NEAR(solved.x[1], -1.0, 1.0e-12);
}

FESA_TEST(assembly_preserves_full_space_stiffness_load_and_reduced_pattern) {
  auto domain = parsedPhase1Domain();
  fesa::DofManager dofs(domain);
  const auto assembly = fesa::assembleSystem(domain, dofs);
  FESA_CHECK(!fesa::hasError(assembly.diagnostics));
  FESA_CHECK(assembly.k_full.rows() == dofs.fullDofCount());
  FESA_CHECK(assembly.k_full.cols() == dofs.fullDofCount());
  FESA_CHECK(assembly.f_full.size() == static_cast<std::size_t>(dofs.fullDofCount()));
  FESA_CHECK(assembly.reduced_pattern.equation_count == dofs.freeDofCount());
  FESA_CHECK(assembly.reduced_pattern.nonzeroCount() == 4);

  const auto node2_uz = dofs.fullIndex(2, fesa::Dof::UZ);
  const auto node3_uz = dofs.fullIndex(3, fesa::Dof::UZ);
  FESA_CHECK_NEAR(assembly.f_full[static_cast<std::size_t>(node2_uz)], -1.0, 1.0e-15);
  FESA_CHECK_NEAR(assembly.f_full[static_cast<std::size_t>(node3_uz)], -1.0, 1.0e-15);
  for (fesa::LocalIndex i = 0; i < assembly.k_full.rows(); ++i) {
    for (fesa::LocalIndex j = 0; j < assembly.k_full.cols(); ++j) {
      FESA_CHECK_NEAR(assembly.k_full(i, j), assembly.k_full(j, i), 1.0e-8);
    }
  }

  const auto reduced = fesa::projectToReducedSystem(assembly, dofs);
  FESA_CHECK(reduced.ok());
  const auto solved = fesa::GaussianEliminationSolver{}.solve(reduced.k, reduced.f);
  FESA_CHECK(solved.ok());
  const auto residual = reduced.k.multiply(solved.x);
  for (std::size_t i = 0; i < residual.size(); ++i) {
    FESA_CHECK_NEAR(residual[i], reduced.f[i], 1.0e-8);
  }
}

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.schema_version == 1);
  FESA_CHECK(result.solver_name == "FESA");
  FESA_CHECK(result.dof_convention == "UX,UY,UZ,RX,RY,RZ");
  FESA_CHECK(result.sign_convention == "Abaqus-compatible");
  FESA_CHECK(result.precision == "double");
  FESA_CHECK(result.index_type == "int64");
  FESA_CHECK(result.element_types == std::vector<std::string>({"MITC4"}));
  FESA_CHECK(result.steps.size() == 1);
  FESA_CHECK(result.steps[0].name == "Step-1");
  FESA_CHECK(result.steps[0].frames.size() == 1);
  const auto& frame = result.steps[0].frames[0];
  FESA_CHECK(frame.frame_id == 0);
  FESA_CHECK(frame.increment == 1);
  FESA_CHECK(frame.iteration == 0);
  FESA_CHECK(frame.converged);
  FESA_CHECK_NEAR(frame.step_time, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(frame.total_time, 1.0, 1.0e-15);
  FESA_CHECK(frame.field_outputs.count("U") == 1);
  FESA_CHECK(frame.field_outputs.count("RF") == 1);
  const auto& u_field = frame.field_outputs.at("U");
  const auto& rf_field = frame.field_outputs.at("RF");
  FESA_CHECK(u_field.name == "U");
  FESA_CHECK(u_field.position == "NODAL");
  FESA_CHECK(u_field.entity_type == "node");
  FESA_CHECK(u_field.basis == "GLOBAL");
  FESA_CHECK(u_field.component_labels == fesa::displacementComponentLabels());
  FESA_CHECK(rf_field.name == "RF");
  FESA_CHECK(rf_field.position == "NODAL");
  FESA_CHECK(rf_field.entity_type == "node");
  FESA_CHECK(rf_field.basis == "GLOBAL");
  FESA_CHECK(rf_field.component_labels == fesa::reactionComponentLabels());
  FESA_CHECK(u_field.component_labels[2] == "UZ");
  FESA_CHECK(rf_field.component_labels[2] == "RFZ");
  FESA_CHECK(u_field.entity_ids.size() == domain.nodes.size());
  FESA_CHECK(rf_field.entity_ids.size() == domain.nodes.size());
}

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_csv_loader_reports_required_header_errors) {
  auto table = fesa::loadDisplacementCsvFromString("Node Label,U-U1,U-U2,U-U3,UR-UR1,UR-UR2\n"
                                                   "1,0,0,0,0,0\n",
                                                   "missing-header.csv");
  FESA_CHECK(fesa::containsDiagnostic(table.diagnostics, "FESA-CSV-MISSING-COLUMN"));
}

FESA_TEST(displacement_csv_loader_reports_duplicate_node_rows) {
  auto table = fesa::loadDisplacementCsvFromString("Node Label,U-U1,U-U2,U-U3,UR-UR1,UR-UR2,UR-UR3\n"
                                                   "1,0,0,0,0,0,0\n"
                                                   "1,0,0,0,0,0,0\n",
                                                   "duplicate-node.csv");
  FESA_CHECK(fesa::containsDiagnostic(table.diagnostics, "FESA-CSV-DUPLICATE-NODE"));
}

FESA_TEST(displacement_csv_loader_reports_missing_and_non_numeric_node_rows) {
  auto table = fesa::loadDisplacementCsvFromString("Node Label,U-U1,U-U2,U-U3,UR-UR1,UR-UR2,UR-UR3\n"
                                                   ",0,0,0,0,0,0\n"
                                                   "two,0,0,0,0,0,0\n"
                                                   "3,0,not-a-number,0,0,0,0\n",
                                                   "invalid-node-row.csv");
  FESA_CHECK(diagnosticCount(table.diagnostics, "FESA-CSV-NODE") == 2);
  FESA_CHECK(fesa::containsDiagnostic(table.diagnostics, "FESA-CSV-NUMERIC"));
}

FESA_TEST(displacement_comparator_matches_by_node_id_not_row_order) {
  fesa::FieldOutput actual;
  actual.name = "U";
  actual.position = "NODAL";
  actual.entity_type = "node";
  actual.basis = "GLOBAL";
  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(displacement_comparator_uses_absolute_and_relative_tolerances) {
  fesa::FieldOutput actual;
  actual.name = "U";
  actual.position = "NODAL";
  actual.entity_type = "node";
  actual.basis = "GLOBAL";
  actual.entity_ids = {10};
  actual.component_labels = fesa::displacementComponentLabels();
  actual.values = {{{5.0e-7, 100.0005, 0, 0, 0, 0}}};
  fesa::CsvDisplacementTable expected;
  expected.rows[10] = {10, {0.0, 100.0, 0, 0, 0, 0}};

  auto loose = fesa::compareDisplacements(actual, expected, {1.0e-6, 1.0e-5, 1.0});
  FESA_CHECK(loose.pass);
  FESA_CHECK_NEAR(loose.max_abs_error, 5.0e-4, 1.0e-12);

  auto strict = fesa::compareDisplacements(actual, expected, {1.0e-8, 1.0e-8, 1.0});
  FESA_CHECK(!strict.pass);
  FESA_CHECK(fesa::containsDiagnostic(strict.diagnostics, "FESA-COMPARE-TOLERANCE"));
}

FESA_TEST(displacement_comparator_rejects_wrong_component_labels_and_missing_nodes) {
  fesa::FieldOutput actual;
  actual.name = "U";
  actual.position = "NODAL";
  actual.entity_type = "node";
  actual.basis = "GLOBAL";
  actual.entity_ids = {1};
  actual.component_labels = fesa::reactionComponentLabels();
  actual.values = {{{0, 0, 0, 0, 0, 0}}};
  fesa::CsvDisplacementTable expected;
  expected.rows[2] = {2, {0, 0, 0, 0, 0, 0}};

  auto compared = fesa::compareDisplacements(actual, expected, {1.0e-12, 1.0e-12, 1.0});
  FESA_CHECK(!compared.pass);
  FESA_CHECK(fesa::containsDiagnostic(compared.diagnostics, "FESA-COMPARE-COMPONENT-LABELS"));
  FESA_CHECK(fesa::containsDiagnostic(compared.diagnostics, "FESA-COMPARE-MISSING-ACTUAL"));
}

FESA_TEST(displacement_comparator_reports_duplicate_actual_nodes) {
  fesa::FieldOutput actual;
  actual.name = "U";
  actual.position = "NODAL";
  actual.entity_type = "node";
  actual.basis = "GLOBAL";
  actual.entity_ids = {1, 1};
  actual.component_labels = fesa::displacementComponentLabels();
  actual.values = {{{0, 0, 0, 0, 0, 0}}, {{0, 0, 0, 0, 0, 0}}};
  fesa::CsvDisplacementTable expected;
  expected.rows[1] = {1, {0, 0, 0, 0, 0, 0}};

  auto compared = fesa::compareDisplacements(actual, expected, {1.0e-12, 1.0e-12, 1.0});
  FESA_CHECK(!compared.pass);
  FESA_CHECK(fesa::containsDiagnostic(compared.diagnostics, "FESA-COMPARE-DUPLICATE-ACTUAL"));
}

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

  const auto normalized = parser.parseFile(sourceRoot() + "/references/quad_02_phase1.inp");
  FESA_CHECK(normalized.ok());
  FESA_CHECK(normalized.domain.nodes.size() == 121);
  FESA_CHECK(normalized.domain.elements.size() == 100);
  FESA_CHECK(normalized.domain.steps.size() == 1);
  FESA_CHECK(normalized.domain.steps.front().name == "Step-1");

  const auto reference = fesa::loadDisplacementCsv(sourceRoot() + "/references/quad_02_displacements.csv");
  FESA_CHECK(!fesa::hasError(reference.diagnostics));
  FESA_CHECK(reference.rows.size() == normalized.domain.nodes.size());
  for (const auto& [node_id, node] : normalized.domain.nodes) {
    (void)node;
    FESA_CHECK(reference.rows.count(node_id) == 1);
  }
}

FESA_TEST(quad02_phase1_stored_displacement_reference_regression) {
  const auto input_text = readTextFile(sourceRoot() + "/references/quad_02_phase1.inp");
  const auto analysis = fesa::runLinearStaticInputString(input_text, "quad_02_phase1.inp");
  FESA_CHECK(analysis.ok());
  FESA_CHECK(analysis.state.converged);
  FESA_CHECK(analysis.result_file.steps.size() == 1);
  const auto& frame = analysis.result_file.steps[0].frames[0];
  FESA_CHECK(frame.field_outputs.count("U") == 1);

  const auto expected = fesa::loadDisplacementCsv(sourceRoot() + "/references/quad_02_displacements.csv");
  FESA_CHECK(!fesa::hasError(expected.diagnostics));
  const auto comparison = fesa::compareDisplacements(frame.field_outputs.at("U"), expected, {1.0e-12, 1.0e-5, 1.0});
  checkComparisonPass(comparison);
  FESA_CHECK(comparison.max_abs_error <= 1.0e-5);
  FESA_CHECK(comparison.max_rel_error <= 1.0e-5);
}

FESA_TEST(mitc4_shape_functions_node_order_and_tying_points) {
  auto center = fesa::shapeFunctions(0.0, 0.0);
  const fesa::Real sum = center.n[0] + center.n[1] + center.n[2] + center.n[3];
  FESA_CHECK_NEAR(sum, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(center.dr[0], -0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.dr[1], 0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.dr[2], 0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.dr[3], -0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.ds[0], -0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.ds[1], -0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.ds[2], 0.25, 1.0e-15);
  FESA_CHECK_NEAR(center.ds[3], 0.25, 1.0e-15);

  const auto node_points = fesa::mitc4NodeNaturalCoordinates();
  FESA_CHECK_NEAR(node_points[0].xi, -1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[0].eta, -1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[1].xi, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[1].eta, -1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[2].xi, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[2].eta, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[3].xi, -1.0, 1.0e-15);
  FESA_CHECK_NEAR(node_points[3].eta, 1.0, 1.0e-15);

  for (std::size_t i = 0; i < 4; ++i) {
    const auto corner = fesa::shapeFunctions(node_points[i].xi, node_points[i].eta);
    FESA_CHECK_NEAR(corner.n[i], 1.0, 1.0e-15);
  }

  const auto tying_points = fesa::mitc4TyingPoints();
  FESA_CHECK(tying_points[0].label == "A");
  FESA_CHECK_NEAR(tying_points[0].natural.xi, 0.0, 1.0e-15);
  FESA_CHECK_NEAR(tying_points[0].natural.eta, -1.0, 1.0e-15);
  FESA_CHECK((tying_points[0].edge_node_indices == std::array<fesa::LocalIndex, 2>{0, 1}));
  FESA_CHECK(tying_points[1].label == "B");
  FESA_CHECK_NEAR(tying_points[1].natural.xi, -1.0, 1.0e-15);
  FESA_CHECK_NEAR(tying_points[1].natural.eta, 0.0, 1.0e-15);
  FESA_CHECK((tying_points[1].edge_node_indices == std::array<fesa::LocalIndex, 2>{0, 3}));
  FESA_CHECK(tying_points[2].label == "C");
  FESA_CHECK_NEAR(tying_points[2].natural.xi, 0.0, 1.0e-15);
  FESA_CHECK_NEAR(tying_points[2].natural.eta, 1.0, 1.0e-15);
  FESA_CHECK((tying_points[2].edge_node_indices == std::array<fesa::LocalIndex, 2>{3, 2}));
  FESA_CHECK(tying_points[3].label == "D");
  FESA_CHECK_NEAR(tying_points[3].natural.xi, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(tying_points[3].natural.eta, 0.0, 1.0e-15);
  FESA_CHECK((tying_points[3].edge_node_indices == std::array<fesa::LocalIndex, 2>{1, 2}));
}

FESA_TEST(mitc4_geometry_builds_flat_directors_and_integration_basis) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  checkVecNear(geometry.center_normal, {0, 0, 1}, 1.0e-14);
  checkVecNear(geometry.g1_center, {0.5, 0, 0}, 1.0e-14);
  checkVecNear(geometry.g2_center, {0, 0.5, 0}, 1.0e-14);
  for (const auto& frame : geometry.nodal_frames) {
    checkVecNear(frame.v1, {1, 0, 0}, 1.0e-14);
    checkVecNear(frame.v2, {0, 1, 0}, 1.0e-14);
    checkVecNear(frame.vn, {0, 0, 1}, 1.0e-14);
    checkRightHandedOrthonormal(frame.v1, frame.v2, frame.vn);
  }

  auto basis = fesa::computeMITC4IntegrationBasis(geometry, 0.0, 0.0, 0.0);
  FESA_CHECK(basis.ok());
  checkVecNear(basis.g1, {0.5, 0, 0}, 1.0e-14);
  checkVecNear(basis.g2, {0, 0.5, 0}, 1.0e-14);
  checkVecNear(basis.g3, {0, 0, 0.1}, 1.0e-14);
  checkVecNear(basis.local.e1, {1, 0, 0}, 1.0e-14);
  checkVecNear(basis.local.e2, {0, 1, 0}, 1.0e-14);
  checkVecNear(basis.local.e3, {0, 0, 1}, 1.0e-14);
  checkRightHandedOrthonormal(basis.local.e1, basis.local.e2, basis.local.e3);
  FESA_CHECK_NEAR(basis.jacobian, 0.025, 1.0e-14);
}

FESA_TEST(mitc4_geometry_uses_deterministic_director_fallback_axis) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 0, -1}, {0, 0, -1}}};
  auto geometry = fesa::buildMITC4Geometry(coords, 0.1);
  FESA_CHECK(geometry.ok());
  checkVecNear(geometry.center_normal, {0, 1, 0}, 1.0e-14);
  for (const auto& frame : geometry.nodal_frames) {
    checkVecNear(frame.v1, {-1, 0, 0}, 1.0e-14);
    checkVecNear(frame.v2, {0, 0, 1}, 1.0e-14);
    checkVecNear(frame.vn, {0, 1, 0}, 1.0e-14);
    checkRightHandedOrthonormal(frame.v1, frame.v2, frame.vn);
  }
}

FESA_TEST(mitc4_geometry_reports_singular_geometry_and_thickness) {
  const std::array<fesa::Vec3, 4> flat = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  auto invalid_thickness = fesa::buildMITC4Geometry(flat, 0.0);
  FESA_CHECK(!invalid_thickness.ok());
  FESA_CHECK(fesa::containsDiagnostic(invalid_thickness.diagnostics, "FESA-MITC4-THICKNESS"));

  const std::array<fesa::Vec3, 4> collinear = {{{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {3, 0, 0}}};
  auto singular = fesa::buildMITC4Geometry(collinear, 0.1);
  FESA_CHECK(!singular.ok());
  FESA_CHECK(fesa::containsDiagnostic(singular.diagnostics, "FESA-MITC4-SINGULAR-NORMAL"));

  const std::array<fesa::Vec3, 4> collapsed_corner = {{{0, 0, 0}, {0, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  auto geometry = fesa::buildMITC4Geometry(collapsed_corner, 0.1);
  FESA_CHECK(geometry.ok());
  auto corner_basis = fesa::computeMITC4IntegrationBasis(geometry, -1.0, -1.0, 0.0);
  FESA_CHECK(!corner_basis.ok());
  FESA_CHECK(fesa::containsDiagnostic(corner_basis.diagnostics, "FESA-MITC4-SINGULAR-JACOBIAN"));
}

FESA_TEST(mitc4_rotation_transform_and_displacement_interpolation) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const auto flat_rotation = fesa::mitc4LocalRotations(geometry.nodal_frames[0], {1.0, 2.0, 3.0});
  FESA_CHECK_NEAR(flat_rotation.alpha, 1.0, 1.0e-14);
  FESA_CHECK_NEAR(flat_rotation.beta, 2.0, 1.0e-14);
  FESA_CHECK_NEAR(flat_rotation.gamma, 3.0, 1.0e-14);
  checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 0.0, 5.0}), {0, 0, 0}, 1.0e-14);
  checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 2.0, 5.0}), {2, 0, 0}, 1.0e-14);

  std::array<fesa::Real, 24> dofs = zeroElementDofs();
  for (std::size_t node = 0; node < 4; ++node) {
    dofs[6 * node + 0] = 1.0;
    dofs[6 * node + 1] = 0.5;
    dofs[6 * node + 4] = 2.0;
    dofs[6 * node + 5] = 99.0;
  }
  const auto mid = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 0.0);
  FESA_CHECK(mid.ok());
  checkVecNear(mid.displacement, {1.0, 0.5, 0.0}, 1.0e-14);
  const auto top = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 1.0);
  FESA_CHECK(top.ok());
  checkVecNear(top.displacement, {1.2, 0.5, 0.0}, 1.0e-14);

  const std::array<fesa::Vec3, 4> fallback_coords = {{{0, 0, 0}, {1, 0, 0}, {1, 0, -1}, {0, 0, -1}}};
  const auto fallback_geometry = fesa::buildMITC4Geometry(fallback_coords, 0.1);
  FESA_CHECK(fallback_geometry.ok());
  const auto fallback_rotation = fesa::mitc4LocalRotations(fallback_geometry.nodal_frames[0], {2.0, 3.0, 5.0});
  FESA_CHECK_NEAR(fallback_rotation.alpha, -2.0, 1.0e-14);
  FESA_CHECK_NEAR(fallback_rotation.beta, 5.0, 1.0e-14);
  FESA_CHECK_NEAR(fallback_rotation.gamma, 3.0, 1.0e-14);
}

FESA_TEST(mitc4_direct_covariant_strain_rows_match_finite_difference) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  constexpr fesa::Real xi = 0.2;
  constexpr fesa::Real eta = -0.3;
  constexpr fesa::Real zeta = 0.4;
  const auto rows = fesa::mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta);
  FESA_CHECK(rows.ok());
  FESA_CHECK(fesa::mitc4StrainComponentLabels()[0] == "eps11");
  FESA_CHECK(fesa::mitc4StrainComponentLabels()[3] == "gamma23");
  const std::array<std::size_t, 7> checked_dofs = {0, 1, 2, 7, 9, 10, 11};
  const fesa::Real h = 1.0e-6;
  for (std::size_t dof : checked_dofs) {
    auto plus = zeroElementDofs();
    auto minus = zeroElementDofs();
    plus[dof] = h;
    minus[dof] = -h;
    const auto plus_strain = fesa::mitc4DirectCovariantStrain(geometry, plus, xi, eta, zeta);
    const auto minus_strain = fesa::mitc4DirectCovariantStrain(geometry, minus, xi, eta, zeta);
    FESA_CHECK(plus_strain.ok());
    FESA_CHECK(minus_strain.ok());
    for (std::size_t component = 0; component < 6; ++component) {
      const fesa::Real finite_difference = (plus_strain.values[component] - minus_strain.values[component]) / (2.0 * h);
      FESA_CHECK_NEAR(rows.rows[component][dof], finite_difference, 1.0e-8);
    }
  }

  const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);
  const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
  for (std::size_t node = 0; node < 4; ++node) {
    const std::size_t drilling_dof = 6 * node + 5;
    for (std::size_t component = 0; component < 6; ++component) {
      FESA_CHECK_NEAR(rows.rows[component][drilling_dof], 0.0, 1.0e-14);
    }
  }
  FESA_CHECK(gamma13 == 4);
  FESA_CHECK(gamma23 == 3);
}

FESA_TEST(mitc4_mitc_tying_rows_use_fesa_ac_bd_signs) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
  const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);

  const auto direct_a = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, 0.0);
  const auto direct_b = fesa::mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, 0.0);
  const auto direct_c = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, 0.0);
  const auto direct_d = fesa::mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, 0.0);
  const auto mitc_a = fesa::mitc4TiedCovariantStrainRows(geometry, 0.0, -1.0, 0.0);
  const auto mitc_b = fesa::mitc4TiedCovariantStrainRows(geometry, -1.0, 0.0, 0.0);
  const auto mitc_c = fesa::mitc4TiedCovariantStrainRows(geometry, 0.0, 1.0, 0.0);
  const auto mitc_d = fesa::mitc4TiedCovariantStrainRows(geometry, 1.0, 0.0, 0.0);
  FESA_CHECK(direct_a.ok());
  FESA_CHECK(direct_b.ok());
  FESA_CHECK(direct_c.ok());
  FESA_CHECK(direct_d.ok());
  FESA_CHECK(mitc_a.ok());
  FESA_CHECK(mitc_b.ok());
  FESA_CHECK(mitc_c.ok());
  FESA_CHECK(mitc_d.ok());
  for (std::size_t dof = 0; dof < 24; ++dof) {
    FESA_CHECK_NEAR(mitc_a.rows[gamma13][dof], direct_a.rows[gamma13][dof], 1.0e-14);
    FESA_CHECK_NEAR(mitc_c.rows[gamma13][dof], direct_c.rows[gamma13][dof], 1.0e-14);
    FESA_CHECK_NEAR(mitc_b.rows[gamma23][dof], direct_b.rows[gamma23][dof], 1.0e-14);
    FESA_CHECK_NEAR(mitc_d.rows[gamma23][dof], direct_d.rows[gamma23][dof], 1.0e-14);
  }
}

FESA_TEST(mitc4_mitc_gauss_shear_rows_are_interpolated_from_tying_rows) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  constexpr fesa::Real xi = 0.5;
  constexpr fesa::Real eta = 0.25;
  constexpr fesa::Real zeta = 0.0;
  const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
  const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);
  const auto direct_at_point = fesa::mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta);
  const auto direct_a = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, zeta);
  const auto direct_b = fesa::mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, zeta);
  const auto direct_c = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, zeta);
  const auto direct_d = fesa::mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, zeta);
  const auto tied = fesa::mitc4TiedCovariantStrainRows(geometry, xi, eta, zeta);
  FESA_CHECK(direct_at_point.ok());
  FESA_CHECK(tied.ok());

  const std::size_t node2_ry = 6 + 4;
  const fesa::Real expected_gamma13 =
      0.5 * (1.0 - eta) * direct_a.rows[gamma13][node2_ry] + 0.5 * (1.0 + eta) * direct_c.rows[gamma13][node2_ry];
  FESA_CHECK_NEAR(tied.rows[gamma13][node2_ry], expected_gamma13, 1.0e-14);
  FESA_CHECK(std::fabs(tied.rows[gamma13][node2_ry] - direct_at_point.rows[gamma13][node2_ry]) > 1.0e-4);

  const std::size_t node4_rx = 3 * 6 + 3;
  const fesa::Real expected_gamma23 =
      0.5 * (1.0 - xi) * direct_b.rows[gamma23][node4_rx] + 0.5 * (1.0 + xi) * direct_d.rows[gamma23][node4_rx];
  FESA_CHECK_NEAR(tied.rows[gamma23][node4_rx], expected_gamma23, 1.0e-14);
  FESA_CHECK(std::fabs(tied.rows[gamma23][node4_rx] - direct_at_point.rows[gamma23][node4_rx]) > 1.0e-4);

  for (std::size_t component : {std::size_t{0}, std::size_t{1}, std::size_t{2}, std::size_t{5}}) {
    for (std::size_t dof = 0; dof < 24; ++dof) {
      FESA_CHECK_NEAR(tied.rows[component][dof], direct_at_point.rows[component][dof], 1.0e-14);
    }
  }
}

FESA_TEST(mitc4_plane_stress_material_matrix_uses_documented_order_and_shear_correction) {
  constexpr fesa::Real elastic_modulus = 210.0;
  constexpr fesa::Real poisson_ratio = 0.3;
  constexpr fesa::Real kappa = 5.0 / 6.0;
  const auto law = fesa::mitc4PlaneStressMaterialMatrix(elastic_modulus, poisson_ratio, kappa);
  FESA_CHECK(law.ok());
  const auto eps11 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps11);
  const auto eps22 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps22);
  const auto eps33 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps33);
  const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23);
  const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13);
  const auto gamma12 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma12);
  const fesa::Real plane_stress_scale = elastic_modulus / (1.0 - poisson_ratio * poisson_ratio);
  const fesa::Real shear_modulus = elastic_modulus / (2.0 * (1.0 + poisson_ratio));

  FESA_CHECK(fesa::mitc4StrainComponentLabels() ==
             (std::array<std::string, 6>{"eps11", "eps22", "eps33", "gamma23", "gamma13", "gamma12"}));
  FESA_CHECK_NEAR(law.matrix[eps11][eps11], plane_stress_scale, 1.0e-12);
  FESA_CHECK_NEAR(law.matrix[eps11][eps22], poisson_ratio * plane_stress_scale, 1.0e-12);
  FESA_CHECK_NEAR(law.matrix[eps22][eps11], poisson_ratio * plane_stress_scale, 1.0e-12);
  FESA_CHECK_NEAR(law.matrix[eps22][eps22], plane_stress_scale, 1.0e-12);
  for (std::size_t component = 0; component < 6; ++component) {
    FESA_CHECK_NEAR(law.matrix[eps33][component], 0.0, 1.0e-15);
    FESA_CHECK_NEAR(law.matrix[component][eps33], 0.0, 1.0e-15);
  }
  FESA_CHECK_NEAR(law.matrix[gamma23][gamma23], kappa * shear_modulus, 1.0e-12);
  FESA_CHECK_NEAR(law.matrix[gamma13][gamma13], kappa * shear_modulus, 1.0e-12);
  FESA_CHECK_NEAR(law.matrix[gamma12][gamma12], shear_modulus, 1.0e-12);
  FESA_CHECK_NEAR(law.matrix[gamma13][gamma23], 0.0, 1.0e-15);

  fesa::MITC4StrainVector strain{};
  strain[gamma23] = 2.0;
  strain[gamma13] = 3.0;
  strain[gamma12] = 4.0;
  const auto stress = fesa::multiplyMITC4MaterialMatrix(law.matrix, strain);
  FESA_CHECK_NEAR(stress[gamma23], 2.0 * kappa * shear_modulus, 1.0e-12);
  FESA_CHECK_NEAR(stress[gamma13], 3.0 * kappa * shear_modulus, 1.0e-12);
  FESA_CHECK_NEAR(stress[gamma12], 4.0 * shear_modulus, 1.0e-12);
}

FESA_TEST(mitc4_material_matrix_reports_invalid_elastic_inputs) {
  auto invalid_e = fesa::mitc4PlaneStressMaterialMatrix(-1.0, 0.3);
  FESA_CHECK(!invalid_e.ok());
  FESA_CHECK(fesa::containsDiagnostic(invalid_e.diagnostics, "FESA-MITC4-MATERIAL"));

  auto invalid_nu = fesa::mitc4PlaneStressMaterialMatrix(1000.0, 0.5);
  FESA_CHECK(!invalid_nu.ok());
  FESA_CHECK(fesa::containsDiagnostic(invalid_nu.diagnostics, "FESA-MITC4-POISSON"));

  auto invalid_kappa = fesa::mitc4PlaneStressMaterialMatrix(1000.0, 0.25, 0.0);
  FESA_CHECK(!invalid_kappa.ok());
  FESA_CHECK(fesa::containsDiagnostic(invalid_kappa.diagnostics, "FESA-MITC4-SHEAR-CORRECTION"));
}

FESA_TEST(mitc4_gauss_integration_uses_2x2x2_points_and_unit_weights) {
  const auto points = fesa::mitc4GaussQuadrature2x2x2();
  FESA_CHECK(points.size() == 8);
  const fesa::Real gauss = 1.0 / std::sqrt(3.0);
  fesa::Real total_weight = 0.0;
  int positive_xi = 0;
  int positive_eta = 0;
  int positive_zeta = 0;
  for (const auto& point : points) {
    FESA_CHECK_NEAR(std::fabs(point.xi), gauss, 1.0e-15);
    FESA_CHECK_NEAR(std::fabs(point.eta), gauss, 1.0e-15);
    FESA_CHECK_NEAR(std::fabs(point.zeta), gauss, 1.0e-15);
    FESA_CHECK_NEAR(point.weight, 1.0, 1.0e-15);
    total_weight += point.weight;
    positive_xi += point.xi > 0.0 ? 1 : 0;
    positive_eta += point.eta > 0.0 ? 1 : 0;
    positive_zeta += point.zeta > 0.0 ? 1 : 0;
  }
  FESA_CHECK_NEAR(total_weight, 8.0, 1.0e-15);
  FESA_CHECK(positive_xi == 4);
  FESA_CHECK(positive_eta == 4);
  FESA_CHECK(positive_zeta == 4);
}

FESA_TEST(mitc4_covariant_to_local_transform_is_identity_for_orthonormal_basis) {
  fesa::MITC4IntegrationBasis basis;
  basis.g1 = {1.0, 0.0, 0.0};
  basis.g2 = {0.0, 1.0, 0.0};
  basis.g3 = {0.0, 0.0, 1.0};
  basis.local = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
  basis.jacobian = 1.0;
  const auto transform = fesa::mitc4CovariantToLocalStrainTransform(basis);
  FESA_CHECK(transform.ok());
  for (std::size_t row = 0; row < 6; ++row) {
    for (std::size_t col = 0; col < 6; ++col) {
      FESA_CHECK_NEAR(transform.matrix[row][col], row == col ? 1.0 : 0.0, 1.0e-15);
    }
  }
}

FESA_TEST(mitc4_flat_element_material_transform_scales_convected_strains_to_local_frame) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const auto basis = fesa::computeMITC4IntegrationBasis(geometry, 0.0, 0.0, 0.0);
  FESA_CHECK(basis.ok());
  const auto transform = fesa::mitc4CovariantToLocalStrainTransform(basis);
  FESA_CHECK(transform.ok());

  const std::array<fesa::Real, 6> expected_diagonal = {4.0, 4.0, 100.0, 20.0, 20.0, 4.0};
  for (std::size_t row = 0; row < 6; ++row) {
    for (std::size_t col = 0; col < 6; ++col) {
      FESA_CHECK_NEAR(transform.matrix[row][col], row == col ? expected_diagonal[row] : 0.0, 1.0e-13);
    }
  }

  const auto law = fesa::mitc4PlaneStressMaterialMatrix(1000.0, 0.25);
  FESA_CHECK(law.ok());
  const auto convected = fesa::mitc4TransformMaterialMatrix(law.matrix, transform.matrix);
  fesa::MITC4StrainVector local_strain{};
  local_strain[fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps11)] = 0.02;
  local_strain[fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps22)] = -0.01;
  local_strain[fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma12)] = 0.03;

  fesa::MITC4StrainVector convected_strain{};
  for (std::size_t component = 0; component < 6; ++component) {
    convected_strain[component] = local_strain[component] / expected_diagonal[component];
  }
  const auto local_stress = fesa::multiplyMITC4MaterialMatrix(law.matrix, local_strain);
  const auto convected_stress = fesa::multiplyMITC4MaterialMatrix(convected, convected_strain);
  const fesa::Real local_energy = fesa::dotMITC4Vector(local_strain, local_stress);
  const fesa::Real convected_energy = fesa::dotMITC4Vector(convected_strain, convected_stress);
  FESA_CHECK_NEAR(convected_energy, local_energy, 1.0e-12);
}

FESA_TEST(mitc4_material_integration_samples_carry_tied_rows_basis_and_transformed_material) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const auto data = fesa::mitc4BuildMaterialIntegrationData(geometry, 1000.0, 0.3);
  FESA_CHECK(data.ok());
  FESA_CHECK(data.samples.size() == 8);
  for (const auto& sample : data.samples) {
    FESA_CHECK(sample.ok());
    FESA_CHECK_NEAR(sample.point.weight, 1.0, 1.0e-15);
    FESA_CHECK(sample.basis.ok());
    FESA_CHECK(sample.strain_rows.ok());
    for (std::size_t i = 0; i < 6; ++i) {
      for (std::size_t j = 0; j < 6; ++j) {
        FESA_CHECK_NEAR(sample.convected_material[i][j], sample.convected_material[j][i], 1.0e-10);
      }
    }
  }
}

FESA_TEST(mitc4_element_stiffness_result_is_symmetric_and_uses_2x2x2_samples) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
  FESA_CHECK(result.ok());
  FESA_CHECK(result.integration_point_count == 8);
  FESA_CHECK(result.local_without_drilling.rows() == 24);
  FESA_CHECK(result.local_with_drilling.rows() == 24);
  FESA_CHECK(result.global.rows() == 24);
  FESA_CHECK(result.global.cols() == 24);
  for (fesa::LocalIndex i = 0; i < 24; ++i) {
    for (fesa::LocalIndex j = 0; j < 24; ++j) {
      FESA_CHECK_NEAR(result.local_without_drilling(i, j), result.local_without_drilling(j, i), 1.0e-8);
      FESA_CHECK_NEAR(result.local_with_drilling(i, j), result.local_with_drilling(j, i), 1.0e-8);
      FESA_CHECK_NEAR(result.global(i, j), result.global(j, i), 1.0e-8);
    }
  }
}

FESA_TEST(mitc4_drilling_stiffness_uses_minimum_positive_physical_local_diagonal) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
  FESA_CHECK(result.ok());
  fesa::Real expected_reference = std::numeric_limits<fesa::Real>::infinity();
  for (fesa::LocalIndex node = 0; node < 4; ++node) {
    for (fesa::LocalIndex local_dof = 0; local_dof < 5; ++local_dof) {
      const fesa::Real diagonal = result.local_without_drilling(6 * node + local_dof, 6 * node + local_dof);
      if (std::isfinite(diagonal) && diagonal > 0.0) {
        expected_reference = std::min(expected_reference, diagonal);
      }
    }
  }
  FESA_CHECK(std::isfinite(expected_reference));
  FESA_CHECK_NEAR(result.drilling_reference_diagonal, expected_reference, 1.0e-10);
  FESA_CHECK_NEAR(result.drilling_stiffness_scale, 1.0e-3, 1.0e-15);
  FESA_CHECK_NEAR(result.drilling_stiffness, 1.0e-3 * expected_reference, 1.0e-10);
  for (fesa::LocalIndex node = 0; node < 4; ++node) {
    const fesa::LocalIndex gamma = 6 * node + 5;
    FESA_CHECK_NEAR(result.local_with_drilling(gamma, gamma) - result.local_without_drilling(gamma, gamma),
                    result.drilling_stiffness, 1.0e-12);
  }

  fesa::ElementStiffnessOptions options;
  options.drilling_stiffness_scale = 2.0e-3;
  const auto scaled = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2, options);
  FESA_CHECK(scaled.ok());
  FESA_CHECK_NEAR(scaled.drilling_reference_diagonal, result.drilling_reference_diagonal, 1.0e-10);
  FESA_CHECK_NEAR(scaled.drilling_stiffness, 2.0 * result.drilling_stiffness, 1.0e-10);
}

FESA_TEST(mitc4_drilling_reports_invalid_scale_and_missing_reference_diagonal) {
  fesa::DenseMatrix zero_physical(24, 24);
  auto missing_reference = fesa::mitc4ApplyDrillingStabilization(zero_physical, 1.0e-3);
  FESA_CHECK(!missing_reference.ok());
  FESA_CHECK(fesa::containsDiagnostic(missing_reference.diagnostics, "FESA-MITC4-DRILLING-REFERENCE"));

  auto invalid_scale = fesa::mitc4ApplyDrillingStabilization(zero_physical, -1.0);
  FESA_CHECK(!invalid_scale.ok());
  FESA_CHECK(fesa::containsDiagnostic(invalid_scale.diagnostics, "FESA-MITC4-DRILLING-SCALE"));
}

FESA_TEST(mitc4_rigid_translation_energy_is_zero_and_drilling_energy_is_documented) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
  FESA_CHECK(result.ok());

  auto translate_x = zeroElementVector();
  auto translate_y = zeroElementVector();
  auto translate_z = zeroElementVector();
  for (int node = 0; node < 4; ++node) {
    translate_x[static_cast<std::size_t>(6 * node + 0)] = 1.0;
    translate_y[static_cast<std::size_t>(6 * node + 1)] = 1.0;
    translate_z[static_cast<std::size_t>(6 * node + 2)] = 1.0;
  }
  FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_x)) < 1.0e-8);
  FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_y)) < 1.0e-8);
  FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_z)) < 1.0e-8);

  auto drilling = zeroElementVector();
  for (int node = 0; node < 4; ++node) {
    drilling[static_cast<std::size_t>(6 * node + 5)] = 1.0;
  }
  FESA_CHECK_NEAR(quadraticEnergy(result.global, drilling), 4.0 * result.drilling_stiffness, 1.0e-8);
}

FESA_TEST(mitc4_constant_membrane_patch_strain_is_uniform) {
  const auto coords = flatUnitSquareCoordinates();
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const fesa::Real eps_x = 0.010;
  const fesa::Real eps_y = -0.004;
  const fesa::Real gamma_xy = 0.006;
  const auto dofs = elementDofsFromFields(
      coords,
      [&](const fesa::Vec3& p) {
        return fesa::Vec3{eps_x * p.x + 0.5 * gamma_xy * p.y,
                          eps_y * p.y + 0.5 * gamma_xy * p.x,
                          0.0};
      },
      [](const fesa::Vec3&) { return fesa::Vec3{}; });

  for (const auto& sample : fesa::mitc4GaussQuadrature2x2x2()) {
    const auto strain = localStrainAt(geometry, dofs, sample.xi, sample.eta, sample.zeta);
    FESA_CHECK_NEAR(strain[0], eps_x, 1.0e-12);
    FESA_CHECK_NEAR(strain[1], eps_y, 1.0e-12);
    FESA_CHECK_NEAR(strain[2], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[3], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[4], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[5], gamma_xy, 1.0e-12);
  }

  const auto stiffness = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, 0.2);
  FESA_CHECK(stiffness.ok());
  const auto energy = quadraticEnergy(stiffness.local_without_drilling, toVector(dofs));
  FESA_CHECK(energy > 0.0);
}

FESA_TEST(mitc4_pure_bending_patch_is_membrane_free_and_thickness_antisymmetric) {
  const auto coords = flatUnitSquareCoordinates();
  const fesa::Real thickness = 0.2;
  const auto geometry = fesa::buildMITC4Geometry(coords, thickness);
  FESA_CHECK(geometry.ok());
  const fesa::Real curvature_x = 0.025;
  const auto dofs = elementDofsFromFields(
      coords,
      [](const fesa::Vec3&) { return fesa::Vec3{}; },
      [&](const fesa::Vec3& p) { return fesa::Vec3{0.0, curvature_x * p.x, 0.0}; });

  const fesa::Real xi = -0.30;
  const fesa::Real eta = 0.45;
  const fesa::Real zeta = 1.0 / std::sqrt(3.0);
  const auto top = localStrainAt(geometry, dofs, xi, eta, zeta);
  const auto mid = localStrainAt(geometry, dofs, xi, eta, 0.0);
  const auto bottom = localStrainAt(geometry, dofs, xi, eta, -zeta);

  FESA_CHECK(std::fabs(top[0]) > 1.0e-6);
  FESA_CHECK_NEAR(top[0] + bottom[0], 0.0, 1.0e-12);
  FESA_CHECK_NEAR(mid[0], 0.0, 1.0e-12);
  FESA_CHECK_NEAR(top[1] + bottom[1], 0.0, 1.0e-12);
  FESA_CHECK_NEAR(mid[1], 0.0, 1.0e-12);

  const auto stiffness = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, thickness);
  FESA_CHECK(stiffness.ok());
  const auto energy = quadraticEnergy(stiffness.local_without_drilling, toVector(dofs));
  FESA_CHECK(energy > 0.0);
}

FESA_TEST(mitc4_pure_transverse_shear_patch_is_constant) {
  const auto coords = flatUnitSquareCoordinates();
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const fesa::Real gamma_xz = 0.012;
  const auto dofs = elementDofsFromFields(
      coords,
      [&](const fesa::Vec3& p) { return fesa::Vec3{0.0, 0.0, gamma_xz * p.x}; },
      [](const fesa::Vec3&) { return fesa::Vec3{}; });

  for (const auto& sample : fesa::mitc4GaussQuadrature2x2x2()) {
    const auto strain = localStrainAt(geometry, dofs, sample.xi, sample.eta, sample.zeta);
    FESA_CHECK_NEAR(strain[0], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[1], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[2], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[3], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[4], gamma_xz, 1.0e-12);
    FESA_CHECK_NEAR(strain[5], 0.0, 1.0e-12);
  }
}

FESA_TEST(mitc4_pure_twist_patch_reproduces_bilinear_transverse_shear) {
  const auto coords = flatUnitSquareCoordinates();
  const auto geometry = fesa::buildMITC4Geometry(coords, 0.2);
  FESA_CHECK(geometry.ok());
  const fesa::Real twist = 0.018;
  const auto dofs = elementDofsFromFields(
      coords,
      [&](const fesa::Vec3& p) { return fesa::Vec3{0.0, 0.0, twist * p.x * p.y}; },
      [](const fesa::Vec3&) { return fesa::Vec3{}; });

  for (const auto& sample : fesa::mitc4GaussQuadrature2x2x2()) {
    const auto x = 0.5 * (sample.xi + 1.0);
    const auto y = 0.5 * (sample.eta + 1.0);
    const auto strain = localStrainAt(geometry, dofs, sample.xi, sample.eta, sample.zeta);
    FESA_CHECK_NEAR(strain[0], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[1], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[2], 0.0, 1.0e-12);
    FESA_CHECK_NEAR(strain[3], twist * x, 1.0e-12);
    FESA_CHECK_NEAR(strain[4], twist * y, 1.0e-12);
    FESA_CHECK_NEAR(strain[5], 0.0, 1.0e-12);
  }
}

FESA_TEST(mitc4_membrane_patch_energy_is_not_hidden_by_drilling_scale) {
  const auto coords = flatUnitSquareCoordinates();
  const auto dofs = elementDofsFromFields(
      coords,
      [](const fesa::Vec3& p) { return fesa::Vec3{0.003 * p.x, -0.002 * p.y, 0.0}; },
      [](const fesa::Vec3&) { return fesa::Vec3{}; });

  fesa::ElementStiffnessOptions no_drilling;
  no_drilling.drilling_stiffness_scale = 0.0;
  fesa::ElementStiffnessOptions strong_drilling;
  strong_drilling.drilling_stiffness_scale = 1.0e-1;
  const auto physical = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, 0.2, no_drilling);
  const auto stabilized = fesa::mitc4ElementStiffness(coords, 210000.0, 0.30, 0.2, strong_drilling);
  FESA_CHECK(physical.ok());
  FESA_CHECK(stabilized.ok());

  const auto physical_energy = quadraticEnergy(physical.local_with_drilling, toVector(dofs));
  const auto stabilized_energy = quadraticEnergy(stabilized.local_with_drilling, toVector(dofs));
  FESA_CHECK(physical_energy > 0.0);
  FESA_CHECK_NEAR(stabilized_energy, physical_energy, physical_energy * 1.0e-12);
}

FESA_TEST(mitc4_single_element_cantilever_is_thickness_sensitive) {
  const auto thick_tip = singleElementCantileverTipUz(0.20);
  const auto thin_tip = singleElementCantileverTipUz(0.10);
  FESA_CHECK(thick_tip < 0.0);
  FESA_CHECK(thin_tip < thick_tip);
  const auto response_ratio = std::fabs(thin_tip / thick_tip);
  FESA_CHECK(response_ratio > 2.0);
}

FESA_TEST(mitc4_internal_force_is_stiffness_times_displacement_for_linear_phase1) {
  const std::array<fesa::Vec3, 4> coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
  const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2);
  FESA_CHECK(result.ok());
  std::vector<fesa::Real> displacement(24, 0.0);
  for (std::size_t i = 0; i < displacement.size(); ++i) {
    displacement[i] = 0.001 * static_cast<fesa::Real>(i + 1);
  }
  const auto expected = result.global.multiply(displacement);
  const auto actual = fesa::mitc4ElementInternalForce(result, displacement);
  FESA_CHECK(actual.size() == expected.size());
  for (std::size_t i = 0; i < expected.size(); ++i) {
    FESA_CHECK_NEAR(actual[i], expected[i], 1.0e-12);
  }
}

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.state.f_internal_full.size() == result.state.u_full.size());
  FESA_CHECK(result.state.f_external_full.size() == result.state.u_full.size());
  FESA_CHECK(result.state.reaction_full.size() == result.state.u_full.size());
  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);
}

FESA_TEST(linear_static_input_workflow_produces_schema_ready_u_and_rf_fields) {
  const auto result = fesa::runLinearStaticInputString(phase1Input(), "phase1-workflow.inp");
  FESA_CHECK(result.ok());
  FESA_CHECK(result.model.ok());
  FESA_CHECK(result.model.step.name == "Step-1");
  FESA_CHECK(result.model.active_element_ids == std::vector<fesa::GlobalId>({1}));
  FESA_CHECK(result.state.converged);
  FESA_CHECK(result.state.u_full.size() == 24);
  FESA_CHECK(result.state.f_external_full.size() == 24);
  FESA_CHECK(result.state.f_internal_full.size() == 24);
  FESA_CHECK(result.state.reaction_full.size() == 24);

  const auto& result_file = result.result_file;
  FESA_CHECK(result_file.schema_name == "FESA_RESULTS");
  FESA_CHECK(result_file.schema_version == 1);
  FESA_CHECK(result_file.solver_name == "FESA");
  FESA_CHECK(result_file.dof_convention == "UX,UY,UZ,RX,RY,RZ");
  FESA_CHECK(result_file.sign_convention == "Abaqus-compatible");
  FESA_CHECK(result_file.precision == "double");
  FESA_CHECK(result_file.index_type == "int64");
  FESA_CHECK(result_file.node_ids == std::vector<fesa::GlobalId>({1, 2, 3, 4}));
  FESA_CHECK(result_file.element_ids == std::vector<fesa::GlobalId>({1}));
  FESA_CHECK(result_file.element_types == std::vector<std::string>({"MITC4"}));
  FESA_CHECK(result_file.connectivity.front() == (std::array<fesa::GlobalId, 4>{1, 2, 3, 4}));
  FESA_CHECK(result_file.steps.size() == 1);
  FESA_CHECK(result_file.steps[0].name == "Step-1");
  FESA_CHECK(result_file.steps[0].frames.size() == 1);

  const auto& frame = result_file.steps[0].frames[0];
  FESA_CHECK(frame.frame_id == 0);
  FESA_CHECK(frame.increment == 1);
  FESA_CHECK(frame.iteration == 0);
  FESA_CHECK(frame.converged);
  FESA_CHECK_NEAR(frame.step_time, 1.0, 1.0e-15);
  FESA_CHECK_NEAR(frame.total_time, 1.0, 1.0e-15);
  FESA_CHECK(frame.field_outputs.count("U") == 1);
  FESA_CHECK(frame.field_outputs.count("RF") == 1);
  const auto& u = frame.field_outputs.at("U");
  const auto& rf = frame.field_outputs.at("RF");
  FESA_CHECK(u.component_labels == fesa::displacementComponentLabels());
  FESA_CHECK(rf.component_labels == fesa::reactionComponentLabels());
  FESA_CHECK(u.entity_ids == result_file.node_ids);
  FESA_CHECK(rf.entity_ids == result_file.node_ids);
  FESA_CHECK(u.values.size() == result_file.node_ids.size());
  FESA_CHECK(rf.values.size() == result_file.node_ids.size());

  fesa::Real total_rf_z = 0.0;
  for (const auto& values : rf.values) {
    total_rf_z += values[2];
  }
  FESA_CHECK_NEAR(total_rf_z, 2.0, 1.0e-8);
  FESA_CHECK_NEAR(u.values[0][2], 0.0, 1.0e-15);
  FESA_CHECK_NEAR(u.values[3][2], 0.0, 1.0e-15);
}

FESA_TEST(linear_static_input_workflow_routes_parse_errors) {
  const std::string text = R"inp(
*Part, name=P1
*Node
1, 0, 0, 0
*Element, type=S4R
1, 1, 2, 3, 4
*Step, nlgeom=YES
*Static
*End Step
)inp";
  const auto result = fesa::runLinearStaticInputString(text, "unsupported-workflow.inp");
  FESA_CHECK(!result.ok());
  FESA_CHECK(!result.state.converged);
  FESA_CHECK(fesa::containsDiagnostic(result.diagnostics, "FESA-PARSE-UNSUPPORTED-KEYWORD"));
  FESA_CHECK(fesa::containsDiagnostic(result.diagnostics, "FESA-PARSE-UNSUPPORTED-ELEMENT"));
  FESA_CHECK(fesa::containsDiagnostic(result.diagnostics, "FESA-PARSE-UNSUPPORTED-NLGEOM"));
}

FESA_TEST(linear_static_input_workflow_routes_validation_errors_without_bypassing_validator) {
  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
1, 1, 2, 3, 4
*Nset, nset=FIXED
1, 4
*Boundary
FIXED, 1, 6
*Cload
2, 3, -1
*Step, name=Step-1
*Static
*End Step
)inp";
  const auto result = fesa::runLinearStaticInputString(text, "missing-property-workflow.inp");
  FESA_CHECK(!result.ok());
  FESA_CHECK(!result.state.converged);
  FESA_CHECK(fesa::containsDiagnostic(result.diagnostics, "FESA-VALIDATION-MISSING-PROPERTY"));
}

FESA_TEST(linear_static_analysis_uses_solver_adapter_and_reconstructs_full_vectors) {
  auto domain = parsedPhase1Domain();
  RecordingSolver solver({0.25, -0.50});
  fesa::LinearStaticAnalysis analysis(&solver);
  const auto result = analysis.run(domain);
  FESA_CHECK(result.ok());
  FESA_CHECK(solver.called);
  FESA_CHECK(solver.captured_a.rows() == 2);
  FESA_CHECK(solver.captured_a.cols() == 2);
  FESA_CHECK(solver.captured_b.size() == 2);
  FESA_CHECK_NEAR(solver.captured_b[0], -1.0, 1.0e-15);
  FESA_CHECK_NEAR(solver.captured_b[1], -1.0, 1.0e-15);

  fesa::DofManager dofs(domain);
  FESA_CHECK_NEAR(result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(2, fesa::Dof::UZ))], 0.25, 1.0e-15);
  FESA_CHECK_NEAR(result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(3, fesa::Dof::UZ))], -0.50, 1.0e-15);
  FESA_CHECK_NEAR(result.state.u_full[static_cast<std::size_t>(dofs.fullIndex(1, fesa::Dof::UZ))], 0.0, 1.0e-15);
  for (std::size_t i = 0; i < result.state.reaction_full.size(); ++i) {
    FESA_CHECK_NEAR(result.state.reaction_full[i],
                    result.state.f_internal_full[i] - result.state.f_external_full[i], 1.0e-10);
  }
}

FESA_TEST(linear_static_analysis_propagates_solver_singular_diagnostic) {
  auto domain = parsedPhase1Domain();
  FailingSolver solver;
  fesa::LinearStaticAnalysis analysis(&solver);
  const auto result = analysis.run(domain);
  FESA_CHECK(!result.ok());
  FESA_CHECK(!result.state.converged);
  FESA_CHECK(fesa::containsDiagnostic(result.diagnostics, "FESA-SINGULAR-SOLVER"));
}

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;
}