#include "reference_comparison.hpp"

#include "fesa/analysis/analysis_model.hpp"
#include "fesa/analysis/analysis_state.hpp"
#include "fesa/fem/dof_manager.hpp"
#include "fesa/io/abaqus/input_reader.hpp"
#include "fesa/io/hdf5/hdf5_results_writer.hpp"
#include "fesa/model/domain.h"

#include <gtest/gtest.h>

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <fstream>
#include <iterator>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

#ifndef FESA_TEST_SOURCE_DIR
#error FESA_TEST_SOURCE_DIR must identify the repository root.
#endif

#ifndef FESA_TEST_BINARY_DIR
#error FESA_TEST_BINARY_DIR must identify the CMake binary root.
#endif

namespace {

constexpr const char* kInputName = "cantilever beam.inp";
constexpr const char* kDisplacementName = "cantilever beam displacements.csv";
constexpr const char* kReactionName = "cantilever beam reactions.csv";
constexpr const char* kSectionName = "cantilever beam elemental forces.csv";
constexpr const char* kInstanceName = "PART-1_1-1";
constexpr std::size_t kNodeCount = 11U;
constexpr std::size_t kElementCount = 10U;
constexpr std::size_t kExpectedRowCount = 176U;
constexpr std::size_t kExpectedMetricCount = 16U;

using NodalValues = std::array<std::array<double, 6>, kNodeCount>;
using EndpointValues =
    std::array<std::array<std::array<double, 4>, 2>, kElementCount>;

struct ComparisonValues {
    NodalValues displacement{};
    NodalValues reaction{};
    EndpointValues sectionResultants{};
};

const std::filesystem::path& sourceRoot() {
    static const std::filesystem::path root{FESA_TEST_SOURCE_DIR};
    return root;
}

const std::filesystem::path& binaryRoot() {
    static const std::filesystem::path root{FESA_TEST_BINARY_DIR};
    return root;
}

std::string readBytes(const std::filesystem::path& path) {
    std::ifstream stream{path, std::ios::binary};
    if (!stream) {
        throw std::runtime_error{"Unable to read fixture: " + path.string()};
    }
    return {std::istreambuf_iterator<char>{stream},
            std::istreambuf_iterator<char>{}};
}

void writeBytes(const std::filesystem::path& path, const std::string& contents) {
    std::ofstream stream{path, std::ios::binary | std::ios::trunc};
    if (!stream) {
        throw std::runtime_error{"Unable to write build-local fixture: " +
                                 path.string()};
    }
    stream.write(contents.data(), static_cast<std::streamsize>(contents.size()));
    if (!stream) {
        throw std::runtime_error{"Unable to finish build-local fixture write."};
    }
}

std::vector<std::string> readLines(const std::filesystem::path& path) {
    std::ifstream stream{path};
    if (!stream) {
        throw std::runtime_error{"Unable to read fixture lines."};
    }
    std::vector<std::string> lines;
    for (std::string line; std::getline(stream, line);) {
        if (!line.empty() && line.back() == '\r') {
            line.pop_back();
        }
        lines.push_back(std::move(line));
    }
    return lines;
}

void writeLines(
    const std::filesystem::path& path,
    const std::vector<std::string>& lines) {
    std::ofstream stream{path, std::ios::trunc};
    if (!stream) {
        throw std::runtime_error{"Unable to write build-local fixture lines."};
    }
    for (const auto& line : lines) {
        stream << line << '\n';
    }
    if (!stream) {
        throw std::runtime_error{"Unable to finish build-local line write."};
    }
}

void replaceFirst(
    std::string& contents,
    const std::string& from,
    const std::string& to) {
    const std::size_t position = contents.find(from);
    if (position == std::string::npos) {
        throw std::runtime_error{"Fixture token was not found: " + from};
    }
    contents.replace(position, from.size(), to);
}

ComparisonValues referenceValues() {
    ComparisonValues values{};
    const std::array<double, kNodeCount> uz = {
        -1.0e-30,
        -2.761905780e-4,
        -1.066667140e-3,
        -2.314286540e-3,
        -3.961906300e-3,
        -5.952383390e-3,
        -8.228574880e-3,
        -1.073333810e-2,
        -1.340952890e-2,
        -1.620000600e-2,
        -1.904762720e-2};
    const std::array<double, kNodeCount> ury = {
        1.0e-29,
        5.428573350e-4,
        1.028571860e-3,
        1.457143460e-3,
        1.828572130e-3,
        2.142857990e-3,
        2.400001050e-3,
        2.600000940e-3,
        2.742858140e-3,
        2.828572640e-3,
        2.857143990e-3};
    std::array<std::array<double, 4>, kNodeCount> stations{};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        values.displacement[node][2U] = uz[node];
        values.displacement[node][4U] = ury[node];
        stations[node][2U] = node < kElementCount
                                 ? 1.0e7 - 1.0e6 * static_cast<double>(node)
                                 : -1.56e-2;
    }
    values.reaction[0U][2U] = 1.0e6;
    values.reaction[0U][4U] = -1.0e7;
    for (std::size_t element = 0U; element < kElementCount; ++element) {
        values.sectionResultants[element][0U] = stations[element];
        values.sectionResultants[element][1U] = stations[element + 1U];
    }
    return values;
}

fesa::ModelDefinition makeDefinition(
    const std::filesystem::path& input,
    std::string sourceContentIdentity) {
    fesa::ModelDefinition definition{};
    definition.source_path = input;
    definition.source_content_identity = std::move(sourceContentIdentity);
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        const auto label = static_cast<std::int64_t>(node + 1U);
        definition.nodes.push_back({
            {kInstanceName, label, std::to_string(label)},
            {static_cast<double>(node), 0.0, 0.0},
            {input, node + 1U}});
    }
    definition.materials.push_back(
        {"Material-1", 2.1e11, 0.3, {input, 20U}});
    definition.sections.push_back({
        "Section-1",
        1.0,
        0.0833333,
        0.0,
        0.0833333,
        0.140833,
        {0.0, 1.0, 0.0},
        {},
        {input, 30U}});
    for (std::size_t element = 0U; element < kElementCount; ++element) {
        const auto label = static_cast<std::int64_t>(element + 1U);
        definition.elements.push_back({
            {kInstanceName, label, std::to_string(label)},
            {static_cast<fesa::EntityIndex>(element),
             static_cast<fesa::EntityIndex>(element + 1U)},
            0U,
            0U,
            {input, 40U + element}});
    }
    definition.steps.push_back(
        {"Step-1", {}, {}, 1.0, 1.0, 1.0e-5, 1.0, {input, 60U}});
    return definition;
}

void writeResultsFixture(
    const std::filesystem::path& output,
    const std::filesystem::path& input,
    const ComparisonValues& values) {
    auto parsedInput = fesa::AbaqusInputReader{}.read(input);
    if (!parsedInput.HasValue()) {
        throw std::runtime_error{"Reference fixture input identity read failed."};
    }
    auto domainResult = fesa::Domain::Create(
        makeDefinition(input, parsedInput.Value().sourceContentIdentity));
    if (!domainResult.HasValue()) {
        throw std::runtime_error{"Reference fixture Domain construction failed."};
    }
    fesa::Domain domain = std::move(domainResult.Value());
    auto modelResult = fesa::AnalysisModel::create(domain);
    if (!modelResult.HasValue()) {
        throw std::runtime_error{"Reference fixture AnalysisModel construction failed."};
    }
    fesa::AnalysisModel model = std::move(modelResult.Value());
    auto dofsResult = fesa::DofManager::create(model);
    if (!dofsResult.HasValue()) {
        throw std::runtime_error{"Reference fixture DofManager construction failed."};
    }
    fesa::DofManager dofs = std::move(dofsResult.Value());
    fesa::AnalysisState state =
        fesa::AnalysisState::create(dofs, {"Step-1", 0U});

    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        for (std::size_t component = 0U; component < 6U; ++component) {
            const std::size_t index = node * 6U + component;
            state.displacement()[index] = values.displacement[node][component];
            state.reaction()[index] = values.reaction[node][component];
            state.residual()[index] = values.reaction[node][component];
        }
    }
    for (std::size_t element = 0U; element < kElementCount; ++element) {
        for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
            const std::size_t node = element + endpoint;
            state.endpointResults().push_back({
                static_cast<fesa::EntityIndex>(element),
                static_cast<int>(endpoint),
                domain.Nodes()[node].source_id,
                {},
                values.sectionResultants[element][endpoint]});
        }
        state.gaussResults().push_back(
            {static_cast<fesa::EntityIndex>(element), 1, {}, {}});
        state.gaussResults().push_back(
            {static_cast<fesa::EntityIndex>(element), 2, {}, {}});
        state.stressResults().push_back({
            static_cast<fesa::EntityIndex>(element),
            1,
            0U,
            0.0,
            0.0,
            0.0,
            "fesa-default"});
        state.stressResults().push_back({
            static_cast<fesa::EntityIndex>(element),
            2,
            0U,
            0.0,
            0.0,
            0.0,
            "fesa-default"});
    }

    fesa::Hdf5ResultsWriter writer;
    const fesa::Status status = writer.write(output, domain, state, {});
    if (!status.IsOk()) {
        throw std::runtime_error{"Reference fixture HDF5 write failed."};
    }
}

class ContractFixture {
public:
    ContractFixture(std::string label, const ComparisonValues& values) {
        static std::atomic<std::uint64_t> sequence{0U};
        root_ = binaryRoot() / "reference" / "contract-fixtures" /
                (std::move(label) + "-" +
                 std::to_string(sequence.fetch_add(1U)));
        legacy_ = root_ / "cantilever beam";
        std::error_code error;
        std::filesystem::remove_all(root_, error);
        error.clear();
        if (!std::filesystem::create_directories(legacy_, error) || error) {
            throw std::runtime_error{"Unable to create contract fixture directory."};
        }
        const auto approved = sourceRoot() / "reference" / "cantilever beam";
        for (const char* name :
             {kInputName, kDisplacementName, kReactionName, kSectionName}) {
            std::filesystem::copy_file(
                approved / name,
                legacy_ / name,
                std::filesystem::copy_options::overwrite_existing);
        }
        results_ = root_ / "results.h5";
        writeResultsFixture(results_, legacy_ / kInputName, values);
    }

    ContractFixture(const ContractFixture&) = delete;
    ContractFixture& operator=(const ContractFixture&) = delete;

    ~ContractFixture() {
        std::error_code ignored;
        std::filesystem::remove_all(root_, ignored);
    }

    const std::filesystem::path& root() const noexcept { return root_; }
    const std::filesystem::path& legacy() const noexcept { return legacy_; }
    const std::filesystem::path& results() const noexcept { return results_; }

private:
    std::filesystem::path root_;
    std::filesystem::path legacy_;
    std::filesystem::path results_;
};

void expectFailureCode(
    const fesa::Result<fesa::test::ComparisonReport>& result,
    const std::string& expectedCode) {
    ASSERT_FALSE(result.HasValue());
    ASSERT_FALSE(result.GetStatus().IsOk());
    ASSERT_FALSE(result.GetStatus().Diagnostics().empty());
    EXPECT_EQ(result.GetStatus().Diagnostics().front().code, expectedCode);
}

const fesa::test::RowDecision* findRow(
    const fesa::test::ComparisonReport& report,
    const fesa::test::ComparisonQuantity quantity,
    const std::int64_t sourceNodeLabel,
    const std::string& component) {
    const auto found = std::find_if(
        report.rows.begin(),
        report.rows.end(),
        [&](const fesa::test::RowDecision& row) {
            return row.reference.quantity == quantity &&
                   row.reference.sourceNodeLabel == sourceNodeLabel &&
                   row.reference.component == component;
        });
    return found == report.rows.end() ? nullptr : &*found;
}

const fesa::test::ComponentMetrics* findMetric(
    const fesa::test::ComparisonReport& report,
    const fesa::test::ComparisonQuantity quantity,
    const std::string& component) {
    const auto found = std::find_if(
        report.metrics.begin(),
        report.metrics.end(),
        [&](const fesa::test::ComponentMetrics& metric) {
            return metric.quantity == quantity && metric.component == component;
        });
    return found == report.metrics.end() ? nullptr : &*found;
}

void expectExactRowInventory(const fesa::test::ComparisonReport& report) {
    ASSERT_EQ(report.rows.size(), kExpectedRowCount);
    std::size_t rowIndex = 0U;
    const auto expectRow = [&](const fesa::test::ComparisonQuantity quantity,
                               const std::size_t node,
                               const std::string& component,
                               const std::string& unit,
                               const std::string& coordinateSystem,
                               const std::string& datasetPath) {
        ASSERT_LT(rowIndex, report.rows.size());
        const auto& row = report.rows[rowIndex++];
        for (const auto* side : {&row.fesa, &row.reference}) {
            EXPECT_EQ(side->modelId, "cantilever-beam-b33");
            EXPECT_EQ(side->stepName, "Step-1");
            EXPECT_EQ(side->frameIndex, 0U);
            EXPECT_EQ(side->instanceName, kInstanceName);
            EXPECT_EQ(
                side->sourceNodeLabel,
                static_cast<std::int64_t>(node + 1U));
            EXPECT_EQ(side->quantity, quantity);
            EXPECT_EQ(side->component, component);
            EXPECT_EQ(side->unitDimension, unit);
            EXPECT_EQ(side->coordinateSystem, coordinateSystem);
            EXPECT_EQ(side->hdf5DatasetPath, datasetPath);
        }
    };

    const std::array<std::string, 6> displacementComponents = {
        "UX", "UY", "UZ", "URX", "URY", "URZ"};
    const std::array<std::string, 6> displacementUnits = {
        "length", "length", "length", "radian", "radian", "radian"};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        for (std::size_t component = 0U;
             component < displacementComponents.size();
             ++component) {
            expectRow(
                fesa::test::ComparisonQuantity::displacement,
                node,
                displacementComponents[component],
                displacementUnits[component],
                "global-cartesian",
                "/steps/Step-1/frames/0/nodal/displacement");
        }
    }

    const std::array<std::string, 6> reactionComponents = {
        "RF1", "RF2", "RF3", "RM1", "RM2", "RM3"};
    const std::array<std::string, 6> reactionUnits = {
        "force",
        "force",
        "force",
        "force*length",
        "force*length",
        "force*length"};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        for (std::size_t component = 0U;
             component < reactionComponents.size();
             ++component) {
            expectRow(
                fesa::test::ComparisonQuantity::reaction,
                node,
                reactionComponents[component],
                reactionUnits[component],
                "global-cartesian",
                "/steps/Step-1/frames/0/nodal/reaction");
        }
    }

    const std::array<std::string, 4> sectionComponents = {
        "N", "T", "My", "Mz"};
    const std::array<std::string, 4> sectionUnits = {
        "force", "force*length", "force*length", "force*length"};
    for (std::size_t node = 0U; node < kNodeCount; ++node) {
        for (std::size_t component = 0U;
             component < sectionComponents.size();
             ++component) {
            expectRow(
                fesa::test::ComparisonQuantity::sectionResultant,
                node,
                sectionComponents[component],
                sectionUnits[component],
                "beam-local",
                "/steps/Step-1/frames/0/element/section_resultant");
        }
    }
    EXPECT_EQ(rowIndex, report.rows.size());
}

} // namespace

TEST(ReferenceComparisonContract,
     PrecheckRejectsMissingSchemaDuplicateAndNonfiniteRows) {
    auto mismatchedValues = referenceValues();
    mismatchedValues.displacement[0U][0U] = 1.0;

    {
        ContractFixture fixture{"missing-file", mismatchedValues};
        ASSERT_TRUE(std::filesystem::remove(
            fixture.legacy() / kDisplacementName));
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "needs-reference-artifacts");
    }
    {
        ContractFixture fixture{"b31", mismatchedValues};
        auto input = readBytes(fixture.legacy() / kInputName);
        replaceFirst(input, "type=B33", "type=B31");
        writeBytes(fixture.legacy() / kInputName, input);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "needs-reference-artifacts");
    }
    {
        ContractFixture fixture{"header", mismatchedValues};
        auto csv = readBytes(fixture.legacy() / kDisplacementName);
        replaceFirst(csv, "U-U1", "U1");
        writeBytes(fixture.legacy() / kDisplacementName, csv);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "schema-mismatch");
    }
    {
        ContractFixture fixture{"missing-row", mismatchedValues};
        auto lines = readLines(fixture.legacy() / kDisplacementName);
        ASSERT_EQ(lines.size(), kNodeCount + 1U);
        lines.pop_back();
        writeLines(fixture.legacy() / kDisplacementName, lines);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "schema-mismatch");
    }
    {
        ContractFixture fixture{"extra-row", mismatchedValues};
        auto lines = readLines(fixture.legacy() / kDisplacementName);
        ASSERT_EQ(lines.size(), kNodeCount + 1U);
        std::string extra = lines.back();
        replaceFirst(extra, ",PART-1_1-1,11,", ",PART-1_1-1,12,");
        lines.push_back(std::move(extra));
        writeLines(fixture.legacy() / kDisplacementName, lines);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "schema-mismatch");
    }
    {
        ContractFixture fixture{"duplicate-row", mismatchedValues};
        auto lines = readLines(fixture.legacy() / kReactionName);
        ASSERT_EQ(lines.size(), kNodeCount + 1U);
        lines.push_back(lines[1U]);
        writeLines(fixture.legacy() / kReactionName, lines);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "schema-mismatch");
    }
    {
        ContractFixture fixture{"nonfinite-row", mismatchedValues};
        auto csv = readBytes(fixture.legacy() / kReactionName);
        replaceFirst(csv, "0.000000000E+00", "NaN");
        writeBytes(fixture.legacy() / kReactionName, csv);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "schema-mismatch");
    }
    {
        ContractFixture fixture{"identity", mismatchedValues};
        auto csv = readBytes(fixture.legacy() / kSectionName);
        replaceFirst(csv, "PART-1_1-1", "WRONG-INSTANCE");
        writeBytes(fixture.legacy() / kSectionName, csv);
        expectFailureCode(
            fesa::test::ReferenceComparison::compare(
                fixture.results(), fixture.legacy()),
            "schema-mismatch");
    }
}

TEST(ReferenceComparisonContract,
     AppliesAbaqusOnlyComponentScaleWithoutClampOrDrop) {
    auto values = referenceValues();
    values.displacement[1U][0U] = 0.999e-9;
    values.displacement[2U][0U] = 1.001e-9;
    values.sectionResultants[0U][0U][2U] = 1.0e7 + 9.0;
    values.sectionResultants[9U][1U][2U] = 0.0;
    ContractFixture fixture{"tolerance", values};

    auto result = fesa::test::ReferenceComparison::compare(
        fixture.results(), fixture.legacy());
    ASSERT_TRUE(result.HasValue());
    const auto& report = result.Value();
    EXPECT_FALSE(report.passed);
    EXPECT_EQ(report.rows.size(), kExpectedRowCount);
    EXPECT_EQ(report.metrics.size(), kExpectedMetricCount);

    const auto* zero = findRow(
        report, fesa::test::ComparisonQuantity::displacement, 1, "UX");
    const auto* nearZero = findRow(
        report, fesa::test::ComparisonQuantity::displacement, 2, "UX");
    const auto* deliberateFailure = findRow(
        report, fesa::test::ComparisonQuantity::displacement, 3, "UX");
    ASSERT_NE(zero, nullptr);
    ASSERT_NE(nearZero, nullptr);
    ASSERT_NE(deliberateFailure, nullptr);
    EXPECT_DOUBLE_EQ(zero->reference.value, 0.0);
    EXPECT_DOUBLE_EQ(zero->fesa.value, 0.0);
    EXPECT_DOUBLE_EQ(zero->tolerance, 1.0e-9);
    EXPECT_TRUE(zero->passed);
    EXPECT_TRUE(nearZero->passed);
    EXPECT_FALSE(deliberateFailure->passed);

    const auto* myMetric = findMetric(
        report, fesa::test::ComparisonQuantity::sectionResultant, "My");
    ASSERT_NE(myMetric, nullptr);
    EXPECT_DOUBLE_EQ(myMetric->referenceScale, 1.0e7);
    const auto* scaled = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 1, "My");
    const auto* residue = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 11, "My");
    ASSERT_NE(scaled, nullptr);
    ASSERT_NE(residue, nullptr);
    EXPECT_DOUBLE_EQ(scaled->tolerance, 10.001);
    EXPECT_DOUBLE_EQ(scaled->absoluteError, 9.0);
    EXPECT_TRUE(scaled->passed);
    EXPECT_DOUBLE_EQ(residue->reference.value, -1.56e-2);
    EXPECT_DOUBLE_EQ(residue->fesa.value, 0.0);
    EXPECT_DOUBLE_EQ(residue->absoluteError, 1.56e-2);
    EXPECT_DOUBLE_EQ(residue->tolerance, 10.001);
    EXPECT_TRUE(residue->passed);
}

TEST(ReferenceComparisonContract,
     ReportsEveryRowAndAggregateMetricDeterministically) {
    auto values = referenceValues();
    values.displacement[0U][0U] = 0.5e-9;
    values.displacement[1U][0U] = -1.0e-9;
    ContractFixture fixture{"metrics", values};

    auto result = fesa::test::ReferenceComparison::compare(
        fixture.results(), fixture.legacy());
    ASSERT_TRUE(result.HasValue());
    const auto& report = result.Value();
    ASSERT_TRUE(report.passed);
    expectExactRowInventory(report);
    ASSERT_EQ(report.metrics.size(), kExpectedMetricCount);
    EXPECT_TRUE(std::all_of(
        report.rows.begin(),
        report.rows.end(),
        [](const fesa::test::RowDecision& row) { return row.passed; }));

    const auto* metric = findMetric(
        report, fesa::test::ComparisonQuantity::displacement, "UX");
    const auto* worst = findRow(
        report, fesa::test::ComparisonQuantity::displacement, 2, "UX");
    ASSERT_NE(metric, nullptr);
    ASSERT_NE(worst, nullptr);
    EXPECT_DOUBLE_EQ(metric->referenceScale, 0.0);
    EXPECT_DOUBLE_EQ(metric->maximumAbsoluteError, 1.0e-9);
    EXPECT_DOUBLE_EQ(metric->maximumNormalizedError, 1.0);
    EXPECT_NEAR(
        metric->rmsError,
        std::sqrt(1.25 / static_cast<double>(kNodeCount)) * 1.0e-9,
        1.0e-21);
    EXPECT_NEAR(metric->normError, std::sqrt(1.25) * 1.0e-9, 1.0e-21);
    EXPECT_EQ(
        metric->worstRow,
        static_cast<std::size_t>(worst - report.rows.data()));

    EXPECT_FALSE(report.stressComparisonApplicable);
    EXPECT_NE(report.stressComparisonReason.find("N/A"), std::string::npos);
    EXPECT_NE(report.stressComparisonReason.find("HDF5"), std::string::npos);
    EXPECT_DOUBLE_EQ(report.physicsEvidence.freeResidualNorm, 0.0);
    EXPECT_EQ(
        report.physicsEvidence.appliedForce,
        (std::array<double, 3>{0.0, 0.0, -1.0e6}));
    EXPECT_EQ(
        report.physicsEvidence.reactionForce,
        (std::array<double, 3>{0.0, 0.0, 1.0e6}));
    EXPECT_EQ(
        report.physicsEvidence.appliedMomentAboutOrigin,
        (std::array<double, 3>{0.0, 1.0e7, 0.0}));
    EXPECT_EQ(
        report.physicsEvidence.reactionMomentAboutOrigin,
        (std::array<double, 3>{0.0, -1.0e7, 0.0}));
    EXPECT_TRUE(report.physicsEvidence.endpointConsistencyPassed);

    const auto jsonA = fixture.root() / "comparison-a.json";
    const auto jsonB = fixture.root() / "comparison-b.json";
    ASSERT_TRUE(
        fesa::test::ReferenceComparison::writeDeterministicJson(report, jsonA)
            .IsOk());
    ASSERT_TRUE(
        fesa::test::ReferenceComparison::writeDeterministicJson(report, jsonB)
            .IsOk());
    const std::string first = readBytes(jsonA);
    EXPECT_EQ(first, readBytes(jsonB));
    for (const char* required : {
             "\"rows\"",
             "\"metrics\"",
             "\"stress_comparison_applicable\":false",
             "\"stress_comparison_reason\"",
             "\"physics_evidence\"",
             "\"free_residual_norm\"",
             "\"applied_force\"",
             "\"reaction_force\"",
             "\"applied_moment_about_origin\"",
             "\"reaction_moment_about_origin\"",
             "\"endpoint_consistency_passed\""}) {
        EXPECT_NE(first.find(required), std::string::npos) << required;
    }
}

TEST(ReferenceComparisonContract,
     NormalizesEligibleStationsWithoutAveraging) {
    auto values = referenceValues();
    values.sectionResultants[0U][0U] = {2.0e-4, 3.0e-4, 1.0e7, 4.0e-4};
    values.sectionResultants[0U][1U][2U] = 9.0e6 - 5.0;
    values.sectionResultants[1U][0U][2U] = 9.0e6 + 5.0;
    ContractFixture fixture{"stations", values};

    auto result = fesa::test::ReferenceComparison::compare(
        fixture.results(), fixture.legacy());
    ASSERT_TRUE(result.HasValue());
    const auto& report = result.Value();
    ASSERT_TRUE(report.passed);
    EXPECT_TRUE(report.physicsEvidence.endpointConsistencyPassed);

    const auto* n = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 1, "N");
    const auto* t = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 1, "T");
    const auto* my = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 1, "My");
    const auto* mz = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 1, "Mz");
    ASSERT_NE(n, nullptr);
    ASSERT_NE(t, nullptr);
    ASSERT_NE(my, nullptr);
    ASSERT_NE(mz, nullptr);
    EXPECT_DOUBLE_EQ(n->fesa.value, 2.0e-4);
    EXPECT_DOUBLE_EQ(t->fesa.value, 3.0e-4);
    EXPECT_DOUBLE_EQ(my->fesa.value, 1.0e7);
    EXPECT_DOUBLE_EQ(mz->fesa.value, 4.0e-4);

    const auto* interior = findRow(
        report, fesa::test::ComparisonQuantity::sectionResultant, 2, "My");
    ASSERT_NE(interior, nullptr);
    EXPECT_DOUBLE_EQ(interior->fesa.value, 9.0e6 - 5.0);
    EXPECT_DOUBLE_EQ(interior->reference.value, 9.0e6);
    EXPECT_DOUBLE_EQ(interior->absoluteError, 5.0);

    auto mismatchValues = referenceValues();
    mismatchValues.sectionResultants[0U][1U][2U] = 9.0e6 - 6.0;
    mismatchValues.sectionResultants[1U][0U][2U] = 9.0e6 + 6.0;
    ContractFixture mismatch{"station-mismatch", mismatchValues};
    expectFailureCode(
        fesa::test::ReferenceComparison::compare(
            mismatch.results(), mismatch.legacy()),
        "tolerance-failure");
}