#include "reference_comparison.hpp"

#include "fesa/analysis/analysis_model.hpp"
#include "fesa/assembly/load_assembler.hpp"
#include "fesa/fem/dof_manager.hpp"
#include "fesa/io/abaqus/domain_mapper.hpp"
#include "fesa/io/abaqus/input_reader.hpp"
#include "fesa/results/result_recovery.hpp"

#include <hdf5.h>

#include <algorithm>
#include <array>
#include <cerrno>
#include <cctype>
#include <charconv>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iterator>
#include <limits>
#include <locale>
#include <sstream>
#include <stdexcept>
#include <string>
#include <system_error>
#include <utility>
#include <vector>

namespace fesa::test {
namespace {

constexpr const char* kModelId = "cantilever-beam-b33";
constexpr const char* kStepName = "Step-1";
constexpr std::size_t kFrameIndex = 0U;
constexpr const char* kFrameText = "Increment 1: Step Time = 1.000";
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* kDisplacementPath =
    "/steps/Step-1/frames/0/nodal/displacement";
constexpr const char* kReactionPath =
    "/steps/Step-1/frames/0/nodal/reaction";
constexpr const char* kSectionPath =
    "/steps/Step-1/frames/0/element/section_resultant";
constexpr const char* kStressPath =
    "/steps/Step-1/frames/0/element/stress_s11";
constexpr double kKinematicFloor = 1.0e-9;
constexpr double kForceMomentFloor = 1.0e-3;
constexpr double kRelativeCoefficient = 1.0e-6;

class ComparisonFailure final : public std::runtime_error {
public:
    ComparisonFailure(std::string code, std::string message)
        : std::runtime_error{std::move(message)}, code_{std::move(code)} {}

    const std::string& code() const noexcept { return code_; }

private:
    std::string code_;
};

[[noreturn]] void fail(const std::string& code, const std::string& message) {
    throw ComparisonFailure{code, message};
}

Status comparisonFailureStatus(
    const std::string& code, const std::string& message) {
    return Status::failure(
        FailureCategory::model,
        {{Severity::error, code, {}, "", kModelId, message}});
}

std::string trim(const std::string& value) {
    const auto isSpace = [](const unsigned char character) {
        return std::isspace(character) != 0;
    };
    const auto begin = std::find_if_not(
        value.begin(), value.end(), [&](const char character) {
            return isSpace(static_cast<unsigned char>(character));
        });
    const auto end = std::find_if_not(
        value.rbegin(), value.rend(), [&](const char character) {
            return isSpace(static_cast<unsigned char>(character));
        }).base();
    return begin < end ? std::string{begin, end} : std::string{};
}

std::string collapseWhitespace(const std::string& value) {
    std::string result;
    bool pendingSpace = false;
    for (const char character : trim(value)) {
        if (std::isspace(static_cast<unsigned char>(character)) != 0) {
            pendingSpace = !result.empty();
        } else {
            if (pendingSpace) {
                result.push_back(' ');
            }
            result.push_back(character);
            pendingSpace = false;
        }
    }
    return result;
}

std::string asciiLower(std::string value) {
    std::transform(
        value.begin(), value.end(), value.begin(), [](const char character) {
            if (character >= 'A' && character <= 'Z') {
                return static_cast<char>(character - 'A' + 'a');
            }
            return character;
        });
    return value;
}

std::vector<std::string> splitCsvLine(const std::string& line) {
    std::vector<std::string> fields;
    std::size_t start = 0U;
    while (true) {
        const std::size_t comma = line.find(',', start);
        fields.push_back(trim(line.substr(start, comma - start)));
        if (comma == std::string::npos) {
            break;
        }
        start = comma + 1U;
    }
    return fields;
}

std::int64_t parsePositiveLabel(const std::string& field) {
    std::int64_t value = 0;
    const char* const begin = field.data();
    const char* const end = begin + field.size();
    const auto parsed = std::from_chars(begin, end, value);
    if (parsed.ec != std::errc{} || parsed.ptr != end || value <= 0) {
        fail("schema-mismatch", "A CSV or HDF5 source-node label is invalid.");
    }
    return value;
}

double parseFiniteDouble(const std::string& field) {
    if (field.empty()) {
        fail("schema-mismatch", "A reference numeric field is empty.");
    }
    errno = 0;
    char* end = nullptr;
    const double value = std::strtod(field.c_str(), &end);
    if (errno == ERANGE || end == field.c_str() || end == nullptr ||
        *end != '\0' || !std::isfinite(value)) {
        fail("schema-mismatch", "A reference numeric field is nonfinite or invalid.");
    }
    return value;
}

struct WideReferenceRow {
    std::string instanceName;
    std::int64_t sourceNodeLabel;
    std::vector<double> values;
};

struct ReferenceTable {
    std::vector<WideReferenceRow> rows;
};

ReferenceTable readReferenceCsv(
    const std::filesystem::path& path,
    const std::vector<std::string>& expectedHeader) {
    std::ifstream stream{path};
    if (!stream) {
        fail("needs-reference-artifacts", "An approved reference CSV is missing.");
    }
    std::string line;
    if (!std::getline(stream, line)) {
        fail("schema-mismatch", "An approved reference CSV is empty.");
    }
    if (!line.empty() && line.back() == '\r') {
        line.pop_back();
    }
    if (splitCsvLine(line) != expectedHeader) {
        fail("schema-mismatch", "An approved reference CSV header is not exact.");
    }

    ReferenceTable table;
    while (std::getline(stream, line)) {
        if (!line.empty() && line.back() == '\r') {
            line.pop_back();
        }
        if (line.empty()) {
            fail("schema-mismatch", "Blank reference CSV rows are not allowed.");
        }
        const auto fields = splitCsvLine(line);
        if (fields.size() != expectedHeader.size() ||
            collapseWhitespace(fields[0U]) != kFrameText ||
            fields[1U].empty()) {
            fail("schema-mismatch", "A reference CSV row has invalid schema or frame identity.");
        }
        WideReferenceRow row{};
        row.instanceName = fields[1U];
        row.sourceNodeLabel = parsePositiveLabel(fields[2U]);
        row.values.reserve(fields.size() - 3U);
        for (std::size_t field = 3U; field < fields.size(); ++field) {
            row.values.push_back(parseFiniteDouble(fields[field]));
        }
        const auto duplicate = std::find_if(
            table.rows.begin(),
            table.rows.end(),
            [&](const WideReferenceRow& existing) {
                return asciiLower(existing.instanceName) ==
                           asciiLower(row.instanceName) &&
                       existing.sourceNodeLabel == row.sourceNodeLabel;
            });
        if (duplicate != table.rows.end()) {
            fail("schema-mismatch", "A reference CSV row identity is duplicated.");
        }
        table.rows.push_back(std::move(row));
    }
    if (table.rows.empty()) {
        fail("schema-mismatch", "An approved reference CSV has no data rows.");
    }
    return table;
}

void requireExactArtifactInventory(
    const std::filesystem::path& legacyDirectory) {
    std::error_code error;
    if (!std::filesystem::is_directory(legacyDirectory, error) || error) {
        fail("needs-reference-artifacts", "The approved legacy directory is missing.");
    }
    std::vector<std::string> names;
    for (std::filesystem::directory_iterator iterator{legacyDirectory, error}, end;
         iterator != end && !error;
         iterator.increment(error)) {
        if (!iterator->is_regular_file(error) || error) {
            fail("needs-reference-artifacts", "The legacy bundle contains a non-file entry.");
        }
        names.push_back(iterator->path().filename().string());
    }
    if (error) {
        fail("needs-reference-artifacts", "The legacy bundle cannot be inspected.");
    }
    std::sort(names.begin(), names.end());
    std::vector<std::string> expected = {
        kDisplacementName, kInputName, kReactionName, kSectionName};
    std::sort(expected.begin(), expected.end());
    if (names != expected) {
        fail("needs-reference-artifacts", "The legacy bundle inventory is not exact.");
    }
}

Domain readApprovedDomain(const std::filesystem::path& inputPath) {
    AbaqusInputReader reader;
    auto parsed = reader.read(inputPath);
    if (!parsed.hasValue()) {
        fail("needs-reference-artifacts", "The approved reference input cannot be parsed.");
    }
    AbaqusDomainMapper mapper;
    auto domain = mapper.map(parsed.value());
    if (!domain.hasValue()) {
        fail(
            "needs-reference-artifacts",
            "The approved reference input is not the required B33 model.");
    }
    return std::move(domain.value());
}

class Hdf5Handle {
public:
    using Closer = herr_t (*)(hid_t);

    Hdf5Handle() = default;
    Hdf5Handle(const hid_t value, Closer closer)
        : value_{value}, closer_{closer} {}
    Hdf5Handle(const Hdf5Handle&) = delete;
    Hdf5Handle& operator=(const Hdf5Handle&) = delete;
    Hdf5Handle(Hdf5Handle&& other) noexcept
        : value_{other.value_}, closer_{other.closer_} {
        other.value_ = -1;
        other.closer_ = nullptr;
    }
    Hdf5Handle& operator=(Hdf5Handle&& other) noexcept {
        if (this != &other) {
            reset();
            value_ = other.value_;
            closer_ = other.closer_;
            other.value_ = -1;
            other.closer_ = nullptr;
        }
        return *this;
    }
    ~Hdf5Handle() { reset(); }

    hid_t get() const noexcept { return value_; }

private:
    void reset() noexcept {
        if (value_ >= 0 && closer_ != nullptr) {
            (void)closer_(value_);
        }
        value_ = -1;
        closer_ = nullptr;
    }

    hid_t value_{-1};
    Closer closer_{nullptr};
};

class Hdf5ErrorSilencer {
public:
    Hdf5ErrorSilencer() {
        if (H5Eget_auto2(H5E_DEFAULT, &callback_, &clientData_) >= 0 &&
            H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr) >= 0) {
            active_ = true;
        }
    }
    Hdf5ErrorSilencer(const Hdf5ErrorSilencer&) = delete;
    Hdf5ErrorSilencer& operator=(const Hdf5ErrorSilencer&) = delete;
    ~Hdf5ErrorSilencer() {
        if (active_) {
            (void)H5Eset_auto2(H5E_DEFAULT, callback_, clientData_);
        }
    }

private:
    H5E_auto2_t callback_{nullptr};
    void* clientData_{nullptr};
    bool active_{false};
};

class Hdf5VlenReclaimer {
public:
    Hdf5VlenReclaimer(
        const hid_t memoryType,
        const hid_t dataSpace,
        void* const data) noexcept
        : memoryType_{memoryType}, dataSpace_{dataSpace}, data_{data} {}
    Hdf5VlenReclaimer(const Hdf5VlenReclaimer&) = delete;
    Hdf5VlenReclaimer& operator=(const Hdf5VlenReclaimer&) = delete;
    ~Hdf5VlenReclaimer() {
        if (active_) {
            (void)H5Dvlen_reclaim(
                memoryType_, dataSpace_, H5P_DEFAULT, data_);
        }
    }

    herr_t reclaim() noexcept {
        active_ = false;
        return H5Dvlen_reclaim(
            memoryType_, dataSpace_, H5P_DEFAULT, data_);
    }

private:
    hid_t memoryType_;
    hid_t dataSpace_;
    void* data_;
    bool active_{true};
};

hid_t requireId(const hid_t value, const char* message) {
    if (value < 0) {
        fail("schema-mismatch", message);
    }
    return value;
}

void requireHdf(const herr_t value, const char* message) {
    if (value < 0) {
        fail("schema-mismatch", message);
    }
}

Hdf5Handle openDataset(const hid_t file, const char* path) {
    return {requireId(H5Dopen2(file, path, H5P_DEFAULT),
                      "A required HDF5 dataset is missing."),
            H5Dclose};
}

std::vector<hsize_t> datasetDimensions(const hid_t dataset) {
    Hdf5Handle space{
        requireId(H5Dget_space(dataset), "Unable to inspect an HDF5 dataspace."),
        H5Sclose};
    const int rank = H5Sget_simple_extent_ndims(space.get());
    if (rank < 0) {
        fail("schema-mismatch", "Unable to inspect an HDF5 dataset rank.");
    }
    std::vector<hsize_t> dimensions(static_cast<std::size_t>(rank));
    if (rank > 0) {
        requireHdf(
            H5Sget_simple_extent_dims(space.get(), dimensions.data(), nullptr),
            "Unable to inspect HDF5 dataset dimensions.");
    }
    return dimensions;
}

std::string readStringAttribute(const hid_t object, const char* name) {
    Hdf5Handle attribute{
        requireId(H5Aopen(object, name, H5P_DEFAULT),
                  "A required HDF5 string attribute is missing."),
        H5Aclose};
    Hdf5Handle type{
        requireId(H5Aget_type(attribute.get()),
                  "Unable to inspect an HDF5 string attribute."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_STRING ||
        H5Tis_variable_str(type.get()) <= 0 ||
        H5Tget_cset(type.get()) != H5T_CSET_UTF8) {
        fail("schema-mismatch", "An HDF5 string attribute has the wrong type.");
    }
    char* raw = nullptr;
    requireHdf(
        H5Aread(attribute.get(), type.get(), &raw),
        "Unable to read an HDF5 string attribute.");
    if (raw == nullptr) {
        fail("schema-mismatch", "An HDF5 string attribute is null.");
    }
    const std::string value{raw};
    requireHdf(H5free_memory(raw), "Unable to release HDF5 string memory.");
    return value;
}

std::uint64_t readUint64Attribute(const hid_t object, const char* name) {
    Hdf5Handle attribute{
        requireId(H5Aopen(object, name, H5P_DEFAULT),
                  "A required HDF5 integer attribute is missing."),
        H5Aclose};
    Hdf5Handle type{
        requireId(H5Aget_type(attribute.get()),
                  "Unable to inspect an HDF5 integer attribute."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_INTEGER ||
        H5Tget_size(type.get()) != sizeof(std::uint64_t) ||
        H5Tget_sign(type.get()) != H5T_SGN_NONE) {
        fail("schema-mismatch", "An HDF5 integer attribute has the wrong type.");
    }
    std::uint64_t value = 0U;
    requireHdf(
        H5Aread(attribute.get(), H5T_NATIVE_UINT64, &value),
        "Unable to read an HDF5 integer attribute.");
    return value;
}

void requireStringAttribute(
    const hid_t object, const char* name, const char* expected) {
    if (readStringAttribute(object, name) != expected) {
        fail("schema-mismatch", "An HDF5 string attribute has the wrong value.");
    }
}

void requireResultAttributes(
    const hid_t dataset,
    const char* components,
    const char* units,
    const char* coordinateSystem,
    const char* location) {
    requireStringAttribute(dataset, "component_names", components);
    requireStringAttribute(dataset, "component_unit_dimensions", units);
    requireStringAttribute(dataset, "coordinate_system", coordinateSystem);
    requireStringAttribute(dataset, "location", location);
    requireStringAttribute(dataset, "step_name", kStepName);
    if (readUint64Attribute(dataset, "frame_index") != kFrameIndex) {
        fail("schema-mismatch", "An HDF5 result has the wrong frame identity.");
    }
}

std::vector<double> readDoubleDataset(
    const hid_t file,
    const char* path,
    const std::vector<hsize_t>& expectedDimensions,
    const char* components,
    const char* units,
    const char* coordinateSystem,
    const char* location) {
    auto dataset = openDataset(file, path);
    if (datasetDimensions(dataset.get()) != expectedDimensions) {
        fail("schema-mismatch", "An HDF5 result dataset has the wrong shape.");
    }
    Hdf5Handle type{
        requireId(H5Dget_type(dataset.get()),
                  "Unable to inspect an HDF5 result type."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_FLOAT ||
        H5Tget_size(type.get()) != sizeof(double) ||
        H5Tequal(type.get(), H5T_IEEE_F64LE) <= 0) {
        fail("schema-mismatch", "An HDF5 result dataset is not float64.");
    }
    requireResultAttributes(
        dataset.get(), components, units, coordinateSystem, location);
    std::size_t count = 1U;
    for (const hsize_t dimension : expectedDimensions) {
        if (dimension > (std::numeric_limits<std::size_t>::max)() / count) {
            fail("schema-mismatch", "An HDF5 result shape overflows size_t.");
        }
        count *= static_cast<std::size_t>(dimension);
    }
    std::vector<double> values(count);
    if (!values.empty()) {
        requireHdf(
            H5Dread(
                dataset.get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                H5P_DEFAULT, values.data()),
            "Unable to read an HDF5 result dataset.");
    }
    if (!std::all_of(values.begin(), values.end(), [](const double value) {
            return std::isfinite(value);
        })) {
        fail("schema-mismatch", "An HDF5 comparison value is nonfinite.");
    }
    return values;
}

void requireCompoundMembers(
    const hid_t dataset, const std::vector<std::string>& expected) {
    Hdf5Handle type{
        requireId(H5Dget_type(dataset), "Unable to inspect an HDF5 compound type."),
        H5Tclose};
    if (H5Tget_class(type.get()) != H5T_COMPOUND ||
        H5Tget_nmembers(type.get()) != static_cast<int>(expected.size())) {
        fail("schema-mismatch", "An HDF5 compound dataset has the wrong schema.");
    }
    for (std::size_t index = 0U; index < expected.size(); ++index) {
        char* raw = H5Tget_member_name(type.get(), static_cast<unsigned>(index));
        if (raw == nullptr) {
            fail("schema-mismatch", "Unable to inspect an HDF5 member name.");
        }
        const std::string actual{raw};
        requireHdf(H5free_memory(raw), "Unable to release an HDF5 member name.");
        if (actual != expected[index]) {
            fail("schema-mismatch", "An HDF5 compound member has the wrong name.");
        }
    }
}

Hdf5Handle makeVariableStringType() {
    Hdf5Handle type{
        requireId(H5Tcopy(H5T_C_S1), "Unable to create an HDF5 string type."),
        H5Tclose};
    requireHdf(
        H5Tset_size(type.get(), H5T_VARIABLE),
        "Unable to size an HDF5 string type.");
    requireHdf(
        H5Tset_cset(type.get(), H5T_CSET_UTF8),
        "Unable to configure an HDF5 string type.");
    return type;
}

struct NodeMemoryRow {
    std::uint64_t internalNodeId;
    char* instanceName;
    char* sourceLabel;
    double coordinates[3];
};

struct ElementMemoryRow {
    std::uint64_t internalElementId;
    char* instanceName;
    char* sourceLabel;
    std::uint64_t nodeInternalIds[2];
    double localAxes[9];
};

struct HdfNode {
    std::uint64_t internalNodeId;
    std::string instanceName;
    std::int64_t sourceNodeLabel;
    std::string sourceNodeLabelText;
    std::array<double, 3> coordinates;
};

struct HdfElement {
    std::uint64_t internalElementId;
    std::string instanceName;
    std::int64_t sourceElementLabel;
    std::string sourceElementLabelText;
    std::array<std::uint64_t, 2> nodeInternalIds;
    std::array<double, 9> localAxes;
};

std::vector<HdfNode> readNodeRows(const hid_t file) {
    auto dataset = openDataset(file, "/model/nodes");
    const auto dimensions = datasetDimensions(dataset.get());
    if (dimensions.size() != 1U || dimensions[0U] == 0U) {
        fail("schema-mismatch", "The HDF5 node table has the wrong shape.");
    }
    requireCompoundMembers(
        dataset.get(),
        {"internal_node_id", "instance_name", "source_label", "coordinates"});
    requireStringAttribute(dataset.get(), "coordinate_system", "global-cartesian");
    requireStringAttribute(dataset.get(), "units_label", "length");

    auto stringType = makeVariableStringType();
    const hsize_t coordinateDimensions[1] = {3U};
    Hdf5Handle coordinateType{
        requireId(
            H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinateDimensions),
            "Unable to create the node coordinate memory type."),
        H5Tclose};
    Hdf5Handle memoryType{
        requireId(H5Tcreate(H5T_COMPOUND, sizeof(NodeMemoryRow)),
                  "Unable to create the node memory type."),
        H5Tclose};
    requireHdf(
        H5Tinsert(memoryType.get(), "internal_node_id",
                  HOFFSET(NodeMemoryRow, internalNodeId), H5T_NATIVE_UINT64),
        "Unable to define the node ID memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "instance_name",
                  HOFFSET(NodeMemoryRow, instanceName), stringType.get()),
        "Unable to define the node instance memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "source_label",
                  HOFFSET(NodeMemoryRow, sourceLabel), stringType.get()),
        "Unable to define the node label memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "coordinates",
                  HOFFSET(NodeMemoryRow, coordinates), coordinateType.get()),
        "Unable to define the node coordinate memory field.");

    std::vector<NodeMemoryRow> raw(static_cast<std::size_t>(dimensions[0U]));
    Hdf5Handle space{
        requireId(H5Dget_space(dataset.get()),
                  "Unable to reopen the node dataspace."),
        H5Sclose};
    requireHdf(
        H5Dread(dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL,
                H5P_DEFAULT, raw.data()),
        "Unable to read the HDF5 node table.");
    Hdf5VlenReclaimer strings{memoryType.get(), space.get(), raw.data()};
    std::vector<HdfNode> rows;
    rows.reserve(raw.size());
    for (const auto& row : raw) {
        if (row.instanceName == nullptr || row.sourceLabel == nullptr) {
            fail("schema-mismatch", "An HDF5 node identity is null.");
        }
        const std::array<double, 3> coordinates = {
            row.coordinates[0U], row.coordinates[1U], row.coordinates[2U]};
        if (!std::all_of(
                coordinates.begin(), coordinates.end(), [](const double value) {
                    return std::isfinite(value);
                })) {
            fail("schema-mismatch", "An HDF5 node coordinate is nonfinite.");
        }
        rows.push_back({
            row.internalNodeId,
            row.instanceName,
            parsePositiveLabel(row.sourceLabel),
            row.sourceLabel,
            coordinates});
    }
    requireHdf(
        strings.reclaim(),
        "Unable to reclaim HDF5 node strings.");
    return rows;
}

std::vector<HdfElement> readElementRows(const hid_t file) {
    auto dataset = openDataset(file, "/model/elements");
    const auto dimensions = datasetDimensions(dataset.get());
    if (dimensions.size() != 1U || dimensions[0U] == 0U) {
        fail("schema-mismatch", "The HDF5 element table has the wrong shape.");
    }
    requireCompoundMembers(
        dataset.get(),
        {"internal_element_id", "instance_name", "source_label",
         "node_internal_ids", "local_axes"});
    requireStringAttribute(dataset.get(), "formulation", "B33-3D-Euler-Bernoulli");

    auto stringType = makeVariableStringType();
    const hsize_t connectivityDimensions[1] = {2U};
    const hsize_t axesDimensions[2] = {3U, 3U};
    Hdf5Handle connectivityType{
        requireId(
            H5Tarray_create2(H5T_NATIVE_UINT64, 1, connectivityDimensions),
            "Unable to create the connectivity memory type."),
        H5Tclose};
    Hdf5Handle axesType{
        requireId(
            H5Tarray_create2(H5T_NATIVE_DOUBLE, 2, axesDimensions),
            "Unable to create the local-axis memory type."),
        H5Tclose};
    Hdf5Handle memoryType{
        requireId(H5Tcreate(H5T_COMPOUND, sizeof(ElementMemoryRow)),
                  "Unable to create the element memory type."),
        H5Tclose};
    requireHdf(
        H5Tinsert(memoryType.get(), "internal_element_id",
                  HOFFSET(ElementMemoryRow, internalElementId),
                  H5T_NATIVE_UINT64),
        "Unable to define the element ID memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "instance_name",
                  HOFFSET(ElementMemoryRow, instanceName), stringType.get()),
        "Unable to define the element instance memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "source_label",
                  HOFFSET(ElementMemoryRow, sourceLabel), stringType.get()),
        "Unable to define the element label memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "node_internal_ids",
                  HOFFSET(ElementMemoryRow, nodeInternalIds),
                  connectivityType.get()),
        "Unable to define the connectivity memory field.");
    requireHdf(
        H5Tinsert(memoryType.get(), "local_axes",
                  HOFFSET(ElementMemoryRow, localAxes), axesType.get()),
        "Unable to define the local-axis memory field.");

    std::vector<ElementMemoryRow> raw(static_cast<std::size_t>(dimensions[0U]));
    Hdf5Handle space{
        requireId(H5Dget_space(dataset.get()),
                  "Unable to reopen the element dataspace."),
        H5Sclose};
    requireHdf(
        H5Dread(dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL,
                H5P_DEFAULT, raw.data()),
        "Unable to read the HDF5 element table.");
    Hdf5VlenReclaimer strings{memoryType.get(), space.get(), raw.data()};
    std::vector<HdfElement> rows;
    rows.reserve(raw.size());
    for (const auto& row : raw) {
        if (row.instanceName == nullptr || row.sourceLabel == nullptr) {
            fail("schema-mismatch", "An HDF5 element identity is null.");
        }
        HdfElement converted{};
        converted.internalElementId = row.internalElementId;
        converted.instanceName = row.instanceName;
        converted.sourceElementLabel = parsePositiveLabel(row.sourceLabel);
        converted.sourceElementLabelText = row.sourceLabel;
        converted.nodeInternalIds = {
            row.nodeInternalIds[0U], row.nodeInternalIds[1U]};
        std::copy(
            std::begin(row.localAxes), std::end(row.localAxes),
            converted.localAxes.begin());
        if (!std::all_of(
                converted.localAxes.begin(), converted.localAxes.end(),
                [](const double value) { return std::isfinite(value); })) {
            fail("schema-mismatch", "An HDF5 local axis is nonfinite.");
        }
        rows.push_back(std::move(converted));
    }
    requireHdf(
        strings.reclaim(),
        "Unable to reclaim HDF5 element strings.");
    return rows;
}

void requireFiniteStress(const hid_t file) {
    auto dataset = openDataset(file, kStressPath);
    const auto dimensions = datasetDimensions(dataset.get());
    if (dimensions.size() != 1U || dimensions[0U] == 0U) {
        fail("schema-mismatch", "The mandatory stress dataset has no rows.");
    }
    requireCompoundMembers(
        dataset.get(),
        {"internal_element_id", "gauss_point_index", "section_point_index",
         "x1", "x2", "source", "S11"});
    requireResultAttributes(
        dataset.get(), "S11", "force/length^2", "beam-local", "section-point");
    struct StressValue {
        double s11;
    };
    Hdf5Handle memoryType{
        requireId(H5Tcreate(H5T_COMPOUND, sizeof(StressValue)),
                  "Unable to create a stress memory type."),
        H5Tclose};
    requireHdf(
        H5Tinsert(memoryType.get(), "S11", HOFFSET(StressValue, s11),
                  H5T_NATIVE_DOUBLE),
        "Unable to define the stress memory field.");
    std::vector<StressValue> values(static_cast<std::size_t>(dimensions[0U]));
    requireHdf(
        H5Dread(dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL,
                H5P_DEFAULT, values.data()),
        "Unable to read the stress dataset.");
    if (!std::all_of(values.begin(), values.end(), [](const StressValue& value) {
            return std::isfinite(value.s11);
        })) {
        fail("schema-mismatch", "The mandatory stress dataset is nonfinite.");
    }
}

std::array<double, 9> expectedLocalAxes(
    const Domain& domain, const EulerBeam3DDefinition& element) {
    const auto& first = domain.nodes()[element.nodeIndices[0U]].coordinates;
    const auto& second = domain.nodes()[element.nodeIndices[1U]].coordinates;
    const auto& guide = domain.sections()[element.sectionIndex].firstAxis;
    const std::array<double, 3> delta = {
        second[0U] - first[0U],
        second[1U] - first[1U],
        second[2U] - first[2U]};
    const double length = std::hypot(delta[0U], delta[1U], delta[2U]);
    const std::array<double, 3> x = {
        delta[0U] / length, delta[1U] / length, delta[2U] / length};
    const double projection =
        guide[0U] * x[0U] + guide[1U] * x[1U] + guide[2U] * x[2U];
    const std::array<double, 3> yTrial = {
        guide[0U] - projection * x[0U],
        guide[1U] - projection * x[1U],
        guide[2U] - projection * x[2U]};
    const double yNorm = std::hypot(yTrial[0U], yTrial[1U], yTrial[2U]);
    const std::array<double, 3> y = {
        yTrial[0U] / yNorm, yTrial[1U] / yNorm, yTrial[2U] / yNorm};
    const std::array<double, 3> z = {
        x[1U] * y[2U] - x[2U] * y[1U],
        x[2U] * y[0U] - x[0U] * y[2U],
        x[0U] * y[1U] - x[1U] * y[0U]};
    return {
        x[0U], x[1U], x[2U],
        y[0U], y[1U], y[2U],
        z[0U], z[1U], z[2U]};
}

struct HdfProjection {
    std::vector<HdfNode> nodes;
    std::vector<HdfElement> elements;
    std::vector<double> displacement;
    std::vector<double> reaction;
    std::vector<double> sectionResultants;
};

HdfProjection readHdfProjection(
    const std::filesystem::path& results,
    const std::filesystem::path& input,
    const Domain& domain) {
    std::error_code error;
    if (!std::filesystem::is_regular_file(results, error) || error) {
        fail("needs-solver-results", "The authoritative FESA results.h5 is missing.");
    }
    Hdf5ErrorSilencer silence;
    if (H5Fis_hdf5(results.string().c_str()) <= 0) {
        fail("schema-mismatch", "The solver result is not an HDF5 file.");
    }
    Hdf5Handle file{
        requireId(H5Fopen(results.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT),
                  "Unable to open the solver HDF5 file read-only."),
        H5Fclose};
    Hdf5Handle metadata{
        requireId(H5Gopen2(file.get(), "/metadata", H5P_DEFAULT),
                  "The HDF5 metadata group is missing."),
        H5Gclose};
    if (readUint64Attribute(metadata.get(), "schema_version") != 0U ||
        readUint64Attribute(metadata.get(), "frame_index") != kFrameIndex) {
        fail("schema-mismatch", "The HDF5 schema or frame version is wrong.");
    }
    requireStringAttribute(
        metadata.get(), "feature_id", "linear-static-3d-euler-beam");
    requireStringAttribute(
        metadata.get(), "unit_system_label", "user-consistent-unspecified");
    requireStringAttribute(
        metadata.get(), "coordinate_convention",
        "global-cartesian; beam-local=(t,n1,t-cross-n1)");
    requireStringAttribute(
        metadata.get(), "element_formulation", "B33-3D-Euler-Bernoulli");
    requireStringAttribute(metadata.get(), "step_name", kStepName);
    const std::string sourceIdentity =
        readStringAttribute(metadata.get(), "source_input_identity");
    const std::string normalizedInput =
        std::filesystem::absolute(input).lexically_normal().generic_u8string();
    const std::string expectedIdentity =
        "path=" + normalizedInput +
        ";content_identity=" + domain.sourceContentIdentity();
    if (sourceIdentity != expectedIdentity) {
        fail("schema-mismatch", "The HDF5 source-input identity is inconsistent.");
    }

    HdfProjection projection{};
    projection.nodes = readNodeRows(file.get());
    projection.elements = readElementRows(file.get());
    if (projection.nodes.size() != domain.nodes().size() ||
        projection.elements.size() != domain.elements().size()) {
        fail("schema-mismatch", "HDF5 model identity counts do not match the input.");
    }
    for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
        const auto& actual = projection.nodes[node];
        const auto& expected = domain.nodes()[node];
        if (actual.internalNodeId != node ||
            actual.instanceName != expected.sourceId.instanceName ||
            actual.sourceNodeLabel != expected.sourceId.sourceLabel ||
            actual.sourceNodeLabelText != expected.sourceId.sourceLabelText ||
            actual.coordinates != expected.coordinates) {
            fail("schema-mismatch", "An HDF5 node identity does not match the input.");
        }
    }
    for (std::size_t element = 0U; element < projection.elements.size(); ++element) {
        const auto& actual = projection.elements[element];
        const auto& expected = domain.elements()[element];
        if (actual.internalElementId != element ||
            actual.instanceName != expected.sourceId.instanceName ||
            actual.sourceElementLabel != expected.sourceId.sourceLabel ||
            actual.sourceElementLabelText != expected.sourceId.sourceLabelText ||
            actual.nodeInternalIds[0U] != expected.nodeIndices[0U] ||
            actual.nodeInternalIds[1U] != expected.nodeIndices[1U]) {
            fail("schema-mismatch", "An HDF5 element identity does not match the input.");
        }
        const auto axes = expectedLocalAxes(domain, expected);
        for (std::size_t component = 0U; component < axes.size(); ++component) {
            if (std::abs(actual.localAxes[component] - axes[component]) > 1.0e-12) {
                fail("schema-mismatch", "An HDF5 local axis does not match the input.");
            }
        }
    }

    const hsize_t nodeCount = static_cast<hsize_t>(projection.nodes.size());
    const hsize_t elementCount = static_cast<hsize_t>(projection.elements.size());
    projection.displacement = readDoubleDataset(
        file.get(), kDisplacementPath, {nodeCount, 6U},
        "UX,UY,UZ,URX,URY,URZ",
        "length,length,length,radian,radian,radian",
        "global-cartesian", "nodal");
    projection.reaction = readDoubleDataset(
        file.get(), kReactionPath, {nodeCount, 6U},
        "RF1,RF2,RF3,RM1,RM2,RM3",
        "force,force,force,force*length,force*length,force*length",
        "global-cartesian", "nodal");
    projection.sectionResultants = readDoubleDataset(
        file.get(), kSectionPath, {elementCount, 2U, 4U},
        "N,T,My,Mz",
        "force,force*length,force*length,force*length",
        "beam-local", "endpoint-positive-local-x-section-cut");
    requireFiniteStress(file.get());
    return projection;
}

std::vector<const WideReferenceRow*> orderedRows(
    const ReferenceTable& table, const std::vector<HdfNode>& nodes) {
    if (table.rows.size() != nodes.size()) {
        fail("schema-mismatch", "The FESA/reference projected row sets differ.");
    }
    std::vector<const WideReferenceRow*> ordered;
    ordered.reserve(nodes.size());
    for (const auto& node : nodes) {
        const auto found = std::find_if(
            table.rows.begin(), table.rows.end(), [&](const WideReferenceRow& row) {
                return asciiLower(row.instanceName) ==
                           asciiLower(node.instanceName) &&
                       row.sourceNodeLabel == node.sourceNodeLabel;
            });
        if (found == table.rows.end() || found->instanceName != node.instanceName) {
            fail("schema-mismatch", "A reference row identity does not match HDF5.");
        }
        ordered.push_back(&*found);
    }
    return ordered;
}

double tableScale(const ReferenceTable& table, const std::size_t valueIndex) {
    double scale = 0.0;
    for (const auto& row : table.rows) {
        if (valueIndex >= row.values.size()) {
            fail("schema-mismatch", "A reference row has the wrong component arity.");
        }
        scale = (std::max)(scale, std::abs(row.values[valueIndex]));
    }
    return scale;
}

std::vector<NodeStationResultRow> normalizeStations(
    const Domain& domain,
    const HdfProjection& hdf,
    const ReferenceTable& sectionTable) {
    auto modelResult = AnalysisModel::create(domain);
    if (!modelResult.hasValue()) {
        fail("schema-mismatch", "The approved input cannot create an analysis view.");
    }
    const AnalysisModel model = std::move(modelResult.value());
    const std::array<double, 4> tolerances = {
        kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 0U),
        kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 3U),
        kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 1U),
        kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 2U)};
    std::vector<EndpointResultRow> endpoints;
    endpoints.reserve(hdf.elements.size() * 2U);
    for (std::size_t element = 0U; element < hdf.elements.size(); ++element) {
        for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
            const auto node = static_cast<std::size_t>(
                hdf.elements[element].nodeInternalIds[endpoint]);
            std::array<double, 4> values{};
            for (std::size_t component = 0U; component < values.size(); ++component) {
                values[component] =
                    hdf.sectionResultants[(element * 2U + endpoint) * 4U + component];
            }
            endpoints.push_back({
                static_cast<EntityIndex>(element),
                static_cast<int>(endpoint),
                domain.nodes()[node].sourceId,
                {},
                values});
        }
    }
    auto normalized = ResultRecovery::normalizeSectionResultantsToNodeStations(
        model, endpoints, tolerances);
    if (!normalized.hasValue()) {
        const auto& diagnostics = normalized.status().diagnostics();
        const std::string code = diagnostics.empty() ? std::string{} : diagnostics[0U].code;
        if (code == "node-station-tolerance-failure") {
            fail("tolerance-failure", "Interior endpoint section resultants disagree.");
        }
        fail("schema-mismatch", "A node station is not eligible for legacy projection.");
    }
    return std::move(normalized.value());
}

const NodeStationResultRow& findStation(
    const std::vector<NodeStationResultRow>& stations,
    const HdfNode& node) {
    const auto found = std::find_if(
        stations.begin(), stations.end(), [&](const NodeStationResultRow& row) {
            return row.node.instanceName == node.instanceName &&
                   row.node.sourceLabel == node.sourceNodeLabel;
        });
    if (found == stations.end()) {
        fail("schema-mismatch", "A projected HDF5 node station is missing.");
    }
    return *found;
}

CanonicalComparisonRow canonicalRow(
    const HdfNode& node,
    const ComparisonQuantity quantity,
    std::string component,
    const double value,
    std::string unit,
    std::string coordinateSystem,
    std::string datasetPath) {
    return {
        kModelId,
        kStepName,
        kFrameIndex,
        node.instanceName,
        node.sourceNodeLabel,
        quantity,
        std::move(component),
        value,
        std::move(unit),
        std::move(coordinateSystem),
        std::move(datasetPath)};
}

void appendNodalRows(
    ComparisonReport& report,
    const HdfProjection& hdf,
    const std::vector<const WideReferenceRow*>& reference,
    const ComparisonQuantity quantity,
    const std::array<std::string, 6>& components,
    const std::array<std::string, 6>& units,
    const std::vector<double>& fesaValues,
    const char* datasetPath) {
    for (std::size_t node = 0U; node < hdf.nodes.size(); ++node) {
        for (std::size_t component = 0U; component < components.size(); ++component) {
            auto fesa = canonicalRow(
                hdf.nodes[node], quantity, components[component],
                fesaValues[node * 6U + component], units[component],
                "global-cartesian", datasetPath);
            auto abaqus = canonicalRow(
                hdf.nodes[node], quantity, components[component],
                reference[node]->values[component], units[component],
                "global-cartesian", datasetPath);
            report.rows.push_back(
                {std::move(fesa), std::move(abaqus), 0.0, 0.0, false});
        }
    }
}

void appendSectionRows(
    ComparisonReport& report,
    const HdfProjection& hdf,
    const std::vector<const WideReferenceRow*>& reference,
    const std::vector<NodeStationResultRow>& stations) {
    const std::array<std::string, 4> components = {"N", "T", "My", "Mz"};
    const std::array<std::string, 4> units = {
        "force", "force*length", "force*length", "force*length"};
    const std::array<std::size_t, 4> referenceColumns = {0U, 3U, 1U, 2U};
    for (std::size_t node = 0U; node < hdf.nodes.size(); ++node) {
        const auto& station = findStation(stations, hdf.nodes[node]);
        for (std::size_t component = 0U; component < components.size(); ++component) {
            auto fesa = canonicalRow(
                hdf.nodes[node], ComparisonQuantity::sectionResultant,
                components[component], station.sectionResultant[component],
                units[component], "beam-local", kSectionPath);
            auto abaqus = canonicalRow(
                hdf.nodes[node], ComparisonQuantity::sectionResultant,
                components[component],
                reference[node]->values[referenceColumns[component]],
                units[component], "beam-local", kSectionPath);
            report.rows.push_back(
                {std::move(fesa), std::move(abaqus), 0.0, 0.0, false});
        }
    }
}

double absoluteFloor(
    const ComparisonQuantity quantity, const std::string&) {
    return quantity == ComparisonQuantity::displacement
               ? kKinematicFloor
               : kForceMomentFloor;
}

void evaluateGroup(
    ComparisonReport& report,
    const ComparisonQuantity quantity,
    const std::string& component) {
    std::vector<std::size_t> rowIndices;
    for (std::size_t index = 0U; index < report.rows.size(); ++index) {
        if (report.rows[index].reference.quantity == quantity &&
            report.rows[index].reference.component == component) {
            rowIndices.push_back(index);
        }
    }
    if (rowIndices.empty()) {
        fail("schema-mismatch", "A canonical comparison component has no rows.");
    }
    double referenceScale = 0.0;
    for (const std::size_t index : rowIndices) {
        referenceScale = (std::max)(
            referenceScale, std::abs(report.rows[index].reference.value));
    }
    const double tolerance =
        absoluteFloor(quantity, component) +
        kRelativeCoefficient * referenceScale;
    double maximumAbsolute = -1.0;
    double maximumNormalized = 0.0;
    std::size_t worstRow = rowIndices.front();
    long double squaredError = 0.0L;
    for (const std::size_t index : rowIndices) {
        auto& row = report.rows[index];
        row.absoluteError = std::abs(row.fesa.value - row.reference.value);
        if (!std::isfinite(row.absoluteError)) {
            fail("schema-mismatch", "A canonical row error is nonfinite.");
        }
        row.tolerance = tolerance;
        row.passed = row.absoluteError <= tolerance;
        report.passed = report.passed && row.passed;
        const double normalized = row.absoluteError / tolerance;
        if (row.absoluteError > maximumAbsolute) {
            maximumAbsolute = row.absoluteError;
            worstRow = index;
        }
        maximumNormalized = (std::max)(maximumNormalized, normalized);
        const long double error = static_cast<long double>(row.absoluteError);
        squaredError += error * error;
    }
    const double normError = std::sqrt(static_cast<double>(squaredError));
    const double rmsError = std::sqrt(
        static_cast<double>(squaredError /
                            static_cast<long double>(rowIndices.size())));
    if (!std::isfinite(normError) || !std::isfinite(rmsError)) {
        fail("schema-mismatch", "A component aggregate error is nonfinite.");
    }
    report.metrics.push_back({
        quantity,
        component,
        referenceScale,
        maximumAbsolute,
        maximumNormalized,
        rmsError,
        normError,
        worstRow});
}

PhysicsEvidence makePhysicsEvidence(
    const Domain& domain,
    const HdfProjection& hdf) {
    auto modelResult = AnalysisModel::create(domain);
    if (!modelResult.hasValue()) {
        fail("schema-mismatch", "The approved input cannot create physics evidence.");
    }
    const AnalysisModel model = std::move(modelResult.value());
    auto dofsResult = DofManager::create(model);
    if (!dofsResult.hasValue()) {
        fail("schema-mismatch", "The approved input cannot create a DOF map.");
    }
    const DofManager dofs = std::move(dofsResult.value());
    auto loadResult = LoadAssembler::assembleFullNodalLoad(model, dofs);
    if (!loadResult.hasValue()) {
        fail("schema-mismatch", "The approved input load cannot be assembled.");
    }
    const Vector load = std::move(loadResult.value());
    if (load.size() != hdf.reaction.size()) {
        fail("schema-mismatch", "The load and reaction spaces are inconsistent.");
    }

    PhysicsEvidence evidence{};
    long double residualSquared = 0.0L;
    for (const std::size_t freeDof : dofs.freeDofs()) {
        const long double value =
            static_cast<long double>(hdf.reaction[freeDof]);
        residualSquared += value * value;
    }
    evidence.freeResidualNorm =
        std::sqrt(static_cast<double>(residualSquared));
    for (std::size_t node = 0U; node < domain.nodes().size(); ++node) {
        const auto& coordinates = domain.nodes()[node].coordinates;
        const std::array<double, 3> applied = {
            load[node * 6U + 0U],
            load[node * 6U + 1U],
            load[node * 6U + 2U]};
        const std::array<double, 3> reaction = {
            hdf.reaction[node * 6U + 0U],
            hdf.reaction[node * 6U + 1U],
            hdf.reaction[node * 6U + 2U]};
        for (std::size_t component = 0U; component < 3U; ++component) {
            evidence.appliedForce[component] += applied[component];
            evidence.reactionForce[component] += reaction[component];
            evidence.appliedMomentAboutOrigin[component] +=
                load[node * 6U + 3U + component];
            evidence.reactionMomentAboutOrigin[component] +=
                hdf.reaction[node * 6U + 3U + component];
        }
        evidence.appliedMomentAboutOrigin[0U] +=
            coordinates[1U] * applied[2U] - coordinates[2U] * applied[1U];
        evidence.appliedMomentAboutOrigin[1U] +=
            coordinates[2U] * applied[0U] - coordinates[0U] * applied[2U];
        evidence.appliedMomentAboutOrigin[2U] +=
            coordinates[0U] * applied[1U] - coordinates[1U] * applied[0U];
        evidence.reactionMomentAboutOrigin[0U] +=
            coordinates[1U] * reaction[2U] - coordinates[2U] * reaction[1U];
        evidence.reactionMomentAboutOrigin[1U] +=
            coordinates[2U] * reaction[0U] - coordinates[0U] * reaction[2U];
        evidence.reactionMomentAboutOrigin[2U] +=
            coordinates[0U] * reaction[1U] - coordinates[1U] * reaction[0U];
    }
    evidence.endpointConsistencyPassed = true;
    return evidence;
}

const char* quantityName(const ComparisonQuantity quantity) {
    switch (quantity) {
    case ComparisonQuantity::displacement:
        return "displacement";
    case ComparisonQuantity::reaction:
        return "reaction";
    case ComparisonQuantity::sectionResultant:
        return "section_resultant";
    }
    return "unknown";
}

void writeJsonString(std::ostream& stream, const std::string& value) {
    static constexpr char digits[] = "0123456789abcdef";
    stream.put('"');
    for (const unsigned char character : value) {
        switch (character) {
        case '"':
            stream << "\\\"";
            break;
        case '\\':
            stream << "\\\\";
            break;
        case '\b':
            stream << "\\b";
            break;
        case '\f':
            stream << "\\f";
            break;
        case '\n':
            stream << "\\n";
            break;
        case '\r':
            stream << "\\r";
            break;
        case '\t':
            stream << "\\t";
            break;
        default:
            if (character < 0x20U) {
                stream << "\\u00" << digits[character >> 4U]
                       << digits[character & 0x0fU];
            } else {
                stream.put(static_cast<char>(character));
            }
            break;
        }
    }
    stream.put('"');
}

void writeCanonicalRow(
    std::ostream& stream, const CanonicalComparisonRow& row) {
    stream << "{\"model_id\":";
    writeJsonString(stream, row.modelId);
    stream << ",\"step_name\":";
    writeJsonString(stream, row.stepName);
    stream << ",\"frame_index\":" << row.frameIndex
           << ",\"instance_name\":";
    writeJsonString(stream, row.instanceName);
    stream << ",\"source_node_label\":" << row.sourceNodeLabel
           << ",\"quantity\":";
    writeJsonString(stream, quantityName(row.quantity));
    stream << ",\"component\":";
    writeJsonString(stream, row.component);
    stream << ",\"value\":" << row.value << ",\"unit_dimension\":";
    writeJsonString(stream, row.unitDimension);
    stream << ",\"coordinate_system\":";
    writeJsonString(stream, row.coordinateSystem);
    stream << ",\"hdf5_dataset_path\":";
    writeJsonString(stream, row.hdf5DatasetPath);
    stream << '}';
}

void writeArray(std::ostream& stream, const std::array<double, 3>& values) {
    stream << '[' << values[0U] << ',' << values[1U] << ',' << values[2U]
           << ']';
}

bool finiteReport(const ComparisonReport& report) {
    const auto finiteArray = [](const std::array<double, 3>& values) {
        return std::all_of(values.begin(), values.end(), [](const double value) {
            return std::isfinite(value);
        });
    };
    if (!std::isfinite(report.physicsEvidence.freeResidualNorm) ||
        !finiteArray(report.physicsEvidence.appliedForce) ||
        !finiteArray(report.physicsEvidence.reactionForce) ||
        !finiteArray(report.physicsEvidence.appliedMomentAboutOrigin) ||
        !finiteArray(report.physicsEvidence.reactionMomentAboutOrigin)) {
        return false;
    }
    for (const auto& row : report.rows) {
        if (!std::isfinite(row.fesa.value) ||
            !std::isfinite(row.reference.value) ||
            !std::isfinite(row.absoluteError) ||
            !std::isfinite(row.tolerance)) {
            return false;
        }
    }
    return std::all_of(
        report.metrics.begin(), report.metrics.end(),
        [](const ComponentMetrics& metric) {
            return std::isfinite(metric.referenceScale) &&
                   std::isfinite(metric.maximumAbsoluteError) &&
                   std::isfinite(metric.maximumNormalizedError) &&
                   std::isfinite(metric.rmsError) &&
                   std::isfinite(metric.normError);
        });
}

} // namespace

Result<ComparisonReport> ReferenceComparison::compare(
    const std::filesystem::path& resultsHdf5,
    const std::filesystem::path& legacyReferenceDirectory) {
    try {
        requireExactArtifactInventory(legacyReferenceDirectory);
        const auto input = legacyReferenceDirectory / kInputName;
        Domain domain = readApprovedDomain(input);
        const ReferenceTable displacement = readReferenceCsv(
            legacyReferenceDirectory / kDisplacementName,
            {"Frame", "Part Instance Name", "Node Label", "U-U1", "U-U2",
             "U-U3", "UR-UR1", "UR-UR2", "UR-UR3"});
        const ReferenceTable reaction = readReferenceCsv(
            legacyReferenceDirectory / kReactionName,
            {"Frame", "Part Instance Name", "Node Label", "RF-RF1", "RF-RF2",
             "RF-RF3", "RM-RM1", "RM-RM2", "RM-RM3"});
        const ReferenceTable section = readReferenceCsv(
            legacyReferenceDirectory / kSectionName,
            {"Frame", "Part Instance Name", "Node Label", "SF-SF1", "SM-SM1",
             "SM-SM2", "SM-SM3"});
        HdfProjection hdf = readHdfProjection(resultsHdf5, input, domain);
        const auto displacementRows = orderedRows(displacement, hdf.nodes);
        const auto reactionRows = orderedRows(reaction, hdf.nodes);
        const auto sectionRows = orderedRows(section, hdf.nodes);
        const auto stations = normalizeStations(domain, hdf, section);
        if (stations.size() != hdf.nodes.size()) {
            fail("schema-mismatch", "The HDF5 node-station row set is incomplete.");
        }

        ComparisonReport report{};
        report.passed = true;
        appendNodalRows(
            report,
            hdf,
            displacementRows,
            ComparisonQuantity::displacement,
            {"UX", "UY", "UZ", "URX", "URY", "URZ"},
            {"length", "length", "length", "radian", "radian", "radian"},
            hdf.displacement,
            kDisplacementPath);
        appendNodalRows(
            report,
            hdf,
            reactionRows,
            ComparisonQuantity::reaction,
            {"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"},
            {"force", "force", "force", "force*length", "force*length",
             "force*length"},
            hdf.reaction,
            kReactionPath);
        appendSectionRows(report, hdf, sectionRows, stations);

        for (const std::string& component :
             {"UX", "UY", "UZ", "URX", "URY", "URZ"}) {
            evaluateGroup(report, ComparisonQuantity::displacement, component);
        }
        for (const std::string& component :
             {"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"}) {
            evaluateGroup(report, ComparisonQuantity::reaction, component);
        }
        for (const std::string& component : {"N", "T", "My", "Mz"}) {
            evaluateGroup(report, ComparisonQuantity::sectionResultant, component);
        }
        report.physicsEvidence = makePhysicsEvidence(domain, hdf);
        report.stressComparisonApplicable = false;
        report.stressComparisonReason =
            "Abaqus beam stress comparison is N/A; analytical/unit and HDF5 "
            "schema tests provide stress evidence.";
        return Result<ComparisonReport>::success(std::move(report));
    } catch (const ComparisonFailure& failure) {
        return Result<ComparisonReport>::failure(
            comparisonFailureStatus(failure.code(), failure.what()));
    } catch (const std::exception& failure) {
        return Result<ComparisonReport>::failure(comparisonFailureStatus(
            "schema-mismatch", failure.what()));
    }
}

Status ReferenceComparison::writeDeterministicJson(
    const ComparisonReport& report,
    const std::filesystem::path& outputJson) {
    if (outputJson.empty() || outputJson.filename().empty() ||
        !finiteReport(report)) {
        return Status::failure(
            FailureCategory::output,
            {{Severity::error,
              "comparison-json-write-failure",
              {},
              "",
              kModelId,
              "The deterministic comparison report or output path is invalid."}});
    }
    std::ofstream stream{outputJson, std::ios::binary | std::ios::trunc};
    if (!stream) {
        return Status::failure(
            FailureCategory::output,
            {{Severity::error,
              "comparison-json-write-failure",
              {},
              "",
              kModelId,
              "The deterministic comparison JSON cannot be opened."}});
    }
    stream.imbue(std::locale::classic());
    stream << std::setprecision(std::numeric_limits<double>::max_digits10);
    stream << "{\"rows\":[";
    for (std::size_t index = 0U; index < report.rows.size(); ++index) {
        if (index != 0U) {
            stream << ',';
        }
        const auto& row = report.rows[index];
        stream << "{\"fesa\":";
        writeCanonicalRow(stream, row.fesa);
        stream << ",\"reference\":";
        writeCanonicalRow(stream, row.reference);
        stream << ",\"absolute_error\":" << row.absoluteError
               << ",\"tolerance\":" << row.tolerance
               << ",\"passed\":" << (row.passed ? "true" : "false") << '}';
    }
    stream << "],\"metrics\":[";
    for (std::size_t index = 0U; index < report.metrics.size(); ++index) {
        if (index != 0U) {
            stream << ',';
        }
        const auto& metric = report.metrics[index];
        stream << "{\"quantity\":";
        writeJsonString(stream, quantityName(metric.quantity));
        stream << ",\"component\":";
        writeJsonString(stream, metric.component);
        stream << ",\"reference_scale\":" << metric.referenceScale
               << ",\"maximum_absolute_error\":"
               << metric.maximumAbsoluteError
               << ",\"maximum_normalized_error\":"
               << metric.maximumNormalizedError
               << ",\"rms_error\":" << metric.rmsError
               << ",\"norm_error\":" << metric.normError
               << ",\"worst_row\":" << metric.worstRow << '}';
    }
    stream << "],\"physics_evidence\":{\"free_residual_norm\":"
           << report.physicsEvidence.freeResidualNorm
           << ",\"applied_force\":";
    writeArray(stream, report.physicsEvidence.appliedForce);
    stream << ",\"reaction_force\":";
    writeArray(stream, report.physicsEvidence.reactionForce);
    stream << ",\"applied_moment_about_origin\":";
    writeArray(stream, report.physicsEvidence.appliedMomentAboutOrigin);
    stream << ",\"reaction_moment_about_origin\":";
    writeArray(stream, report.physicsEvidence.reactionMomentAboutOrigin);
    stream << ",\"endpoint_consistency_passed\":"
           << (report.physicsEvidence.endpointConsistencyPassed ? "true" : "false")
           << "},\"stress_comparison_applicable\":"
           << (report.stressComparisonApplicable ? "true" : "false")
           << ",\"stress_comparison_reason\":";
    writeJsonString(stream, report.stressComparisonReason);
    stream << ",\"passed\":" << (report.passed ? "true" : "false") << "}\n";
    if (!stream) {
        return Status::failure(
            FailureCategory::output,
            {{Severity::error,
              "comparison-json-write-failure",
              {},
              "",
              kModelId,
              "The deterministic comparison JSON write failed."}});
    }
    return Status::ok();
}

} // namespace fesa::test