FESADev/tests/reference/mitc4_reference_comparison_test.cpp

#include "mitc4_reference_comparison.h"

#include <gtest/gtest.h>
#include <hdf5.h>

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

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

namespace {

constexpr const char* kInstanceName = "PART-1-1";
constexpr const char* kInternalFormulation = "FESA-MITC4";
constexpr const char* kIntegrationRule =
    "2x2x2-gauss; mitc4-edge-midpoint-shear";
constexpr const char* kDisplacementPath =
    "/steps/Step-1/frames/0/nodal/displacement";
constexpr std::array<const char*, 6> kComponents{"U1",  "U2",  "U3",
                                                 "UR1", "UR2", "UR3"};

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

  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() {
    if (value_ >= 0 && closer_ != nullptr) {
      (void)closer_(value_);
    }
  }

  hid_t Get() const noexcept { return value_; }

 private:
  hid_t value_;
  Closer closer_;
};

hid_t RequireId(const hid_t value, const char* message) {
  if (value < 0) {
    throw std::runtime_error{message};
  }
  return value;
}

void RequireHdf(const herr_t value, const char* message) {
  if (value < 0) {
    throw std::runtime_error{message};
  }
}

Hdf5Handle MakeUtf8StringType() {
  Hdf5Handle type{RequireId(H5Tcopy(H5T_C_S1), "Unable to copy string type."),
                  H5Tclose};
  RequireHdf(H5Tset_size(type.Get(), H5T_VARIABLE),
             "Unable to create variable string type.");
  RequireHdf(H5Tset_cset(type.Get(), H5T_CSET_UTF8),
             "Unable to create UTF-8 string type.");
  return type;
}

void WriteStringAttribute(const hid_t object, const char* name,
                          const std::string& value) {
  auto type = MakeUtf8StringType();
  Hdf5Handle space{
      RequireId(H5Screate(H5S_SCALAR), "Unable to create attribute space."),
      H5Sclose};
  Hdf5Handle attribute{
      RequireId(H5Acreate2(object, name, type.Get(), space.Get(), H5P_DEFAULT,
                           H5P_DEFAULT),
                "Unable to create string attribute."),
      H5Aclose};
  const char* raw = value.c_str();
  RequireHdf(H5Awrite(attribute.Get(), type.Get(), &raw),
             "Unable to write string attribute.");
}

void WriteUint64Attribute(const hid_t object, const char* name,
                          const std::uint64_t value) {
  Hdf5Handle space{
      RequireId(H5Screate(H5S_SCALAR), "Unable to create attribute space."),
      H5Sclose};
  Hdf5Handle attribute{
      RequireId(H5Acreate2(object, name, H5T_STD_U64LE, space.Get(),
                           H5P_DEFAULT, H5P_DEFAULT),
                "Unable to create integer attribute."),
      H5Aclose};
  RequireHdf(H5Awrite(attribute.Get(), H5T_NATIVE_UINT64, &value),
             "Unable to write integer attribute.");
}

Hdf5Handle CreateGroup(const hid_t parent, const char* path) {
  return {
      RequireId(H5Gcreate2(parent, path, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT),
                "Unable to create HDF5 group."),
      H5Gclose};
}

struct ComparisonValueRow {
  std::string instance_name{kInstanceName};
  std::int64_t source_node_label{};
  std::array<double, 6> values{};
};

std::vector<ComparisonValueRow> DefaultRows() {
  return {{kInstanceName, 1, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}},
          {kInstanceName, 2, {2.0, -4.0, 8.0, 0.1, -0.2, 0.3}}};
}

struct Hdf5Options {
  std::string source_element_type{"S4"};
  std::string internal_formulation{kInternalFormulation};
  std::string integration_rule{kIntegrationRule};
  std::string displacement_components{"UX,UY,UZ,URX,URY,URZ"};
};

struct NodeWriteRow {
  std::uint64_t internal_node_id;
  const char* instance_name;
  const char* source_label;
  double coordinates[3];
};

struct ElementWriteRow {
  std::uint64_t internal_element_id;
  const char* instance_name;
  const char* source_label;
  const char* source_element_type;
  const char* internal_formulation;
  std::uint64_t node_internal_ids[4];
  std::uint64_t shell_section_internal_id;
  std::uint64_t material_internal_id;
};

void WriteNodes(const hid_t file,
                const std::vector<ComparisonValueRow>& values) {
  std::vector<std::string> labels;
  labels.reserve(values.size());
  for (const auto& row : values) {
    labels.push_back(std::to_string(row.source_node_label));
  }
  std::vector<NodeWriteRow> rows;
  rows.reserve(values.size());
  for (std::size_t index = 0U; index < values.size(); ++index) {
    rows.push_back({static_cast<std::uint64_t>(index),
                    values[index].instance_name.c_str(),
                    labels[index].c_str(),
                    {static_cast<double>(index), 0.0, 0.0}});
  }

  auto string_type = MakeUtf8StringType();
  const hsize_t coordinate_dimensions[] = {3U};
  Hdf5Handle file_coordinates{
      RequireId(H5Tarray_create2(H5T_IEEE_F64LE, 1, coordinate_dimensions),
                "Unable to create coordinate file type."),
      H5Tclose};
  Hdf5Handle memory_coordinates{
      RequireId(H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinate_dimensions),
                "Unable to create coordinate memory type."),
      H5Tclose};
  Hdf5Handle file_type{RequireId(H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)),
                                 "Unable to create node file type."),
                       H5Tclose};
  Hdf5Handle memory_type{
      RequireId(H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)),
                "Unable to create node memory type."),
      H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t coordinate_type) {
    RequireHdf(H5Tinsert(type, "internal_node_id",
                         HOFFSET(NodeWriteRow, internal_node_id), integer_type),
               "Unable to define node ID.");
    RequireHdf(
        H5Tinsert(type, "instance_name", HOFFSET(NodeWriteRow, instance_name),
                  string_type.Get()),
        "Unable to define node instance.");
    RequireHdf(
        H5Tinsert(type, "source_label", HOFFSET(NodeWriteRow, source_label),
                  string_type.Get()),
        "Unable to define node label.");
    RequireHdf(H5Tinsert(type, "coordinates",
                         HOFFSET(NodeWriteRow, coordinates), coordinate_type),
               "Unable to define node coordinates.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, file_coordinates.Get());
  insert(memory_type.Get(), H5T_NATIVE_UINT64, memory_coordinates.Get());
  const hsize_t dimensions[] = {static_cast<hsize_t>(rows.size())};
  Hdf5Handle space{RequireId(H5Screate_simple(1, dimensions, nullptr),
                             "Unable to create node space."),
                   H5Sclose};
  Hdf5Handle dataset{
      RequireId(H5Dcreate2(file, "/model/nodes", file_type.Get(), space.Get(),
                           H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT),
                "Unable to create node dataset."),
      H5Dclose};
  RequireHdf(H5Dwrite(dataset.Get(), memory_type.Get(), H5S_ALL, H5S_ALL,
                      H5P_DEFAULT, rows.data()),
             "Unable to write node dataset.");
  WriteStringAttribute(dataset.Get(), "coordinate_system", "global-cartesian");
  WriteStringAttribute(dataset.Get(), "units_label", "length");
}

void WriteElements(const hid_t file, const Hdf5Options& options) {
  ElementWriteRow row{0U,
                      kInstanceName,
                      "1",
                      options.source_element_type.c_str(),
                      options.internal_formulation.c_str(),
                      {0U, 1U, 0U, 1U},
                      0U,
                      0U};
  auto string_type = MakeUtf8StringType();
  const hsize_t node_dimensions[] = {4U};
  Hdf5Handle file_nodes{
      RequireId(H5Tarray_create2(H5T_STD_U64LE, 1, node_dimensions),
                "Unable to create element node file type."),
      H5Tclose};
  Hdf5Handle memory_nodes{
      RequireId(H5Tarray_create2(H5T_NATIVE_UINT64, 1, node_dimensions),
                "Unable to create element node memory type."),
      H5Tclose};
  Hdf5Handle file_type{
      RequireId(H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)),
                "Unable to create element file type."),
      H5Tclose};
  Hdf5Handle memory_type{
      RequireId(H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)),
                "Unable to create element memory type."),
      H5Tclose};
  const auto insert = [&](const hid_t type, const hid_t integer_type,
                          const hid_t nodes_type) {
    RequireHdf(
        H5Tinsert(type, "internal_element_id",
                  HOFFSET(ElementWriteRow, internal_element_id), integer_type),
        "Unable to define element ID.");
    RequireHdf(
        H5Tinsert(type, "instance_name",
                  HOFFSET(ElementWriteRow, instance_name), string_type.Get()),
        "Unable to define element instance.");
    RequireHdf(
        H5Tinsert(type, "source_label", HOFFSET(ElementWriteRow, source_label),
                  string_type.Get()),
        "Unable to define element label.");
    RequireHdf(H5Tinsert(type, "source_element_type",
                         HOFFSET(ElementWriteRow, source_element_type),
                         string_type.Get()),
               "Unable to define source element type.");
    RequireHdf(H5Tinsert(type, "internal_formulation",
                         HOFFSET(ElementWriteRow, internal_formulation),
                         string_type.Get()),
               "Unable to define internal formulation.");
    RequireHdf(
        H5Tinsert(type, "node_internal_ids",
                  HOFFSET(ElementWriteRow, node_internal_ids), nodes_type),
        "Unable to define element connectivity.");
    RequireHdf(H5Tinsert(type, "shell_section_internal_id",
                         HOFFSET(ElementWriteRow, shell_section_internal_id),
                         integer_type),
               "Unable to define section ID.");
    RequireHdf(
        H5Tinsert(type, "material_internal_id",
                  HOFFSET(ElementWriteRow, material_internal_id), integer_type),
        "Unable to define material ID.");
  };
  insert(file_type.Get(), H5T_STD_U64LE, file_nodes.Get());
  insert(memory_type.Get(), H5T_NATIVE_UINT64, memory_nodes.Get());
  const hsize_t dimensions[] = {1U};
  Hdf5Handle space{RequireId(H5Screate_simple(1, dimensions, nullptr),
                             "Unable to create element space."),
                   H5Sclose};
  Hdf5Handle dataset{
      RequireId(H5Dcreate2(file, "/model/elements", file_type.Get(),
                           space.Get(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT),
                "Unable to create element dataset."),
      H5Dclose};
  RequireHdf(H5Dwrite(dataset.Get(), memory_type.Get(), H5S_ALL, H5S_ALL,
                      H5P_DEFAULT, &row),
             "Unable to write element dataset.");
  WriteStringAttribute(dataset.Get(), "formulation", kInternalFormulation);
}

void WriteDisplacement(const hid_t file,
                       const std::vector<ComparisonValueRow>& values,
                       const Hdf5Options& options) {
  std::vector<double> flattened;
  flattened.reserve(values.size() * kComponents.size());
  for (const auto& row : values) {
    flattened.insert(flattened.end(), row.values.begin(), row.values.end());
  }
  const hsize_t dimensions[] = {static_cast<hsize_t>(values.size()),
                                static_cast<hsize_t>(kComponents.size())};
  Hdf5Handle space{RequireId(H5Screate_simple(2, dimensions, nullptr),
                             "Unable to create displacement space."),
                   H5Sclose};
  Hdf5Handle dataset{
      RequireId(H5Dcreate2(file, kDisplacementPath, H5T_IEEE_F64LE, space.Get(),
                           H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT),
                "Unable to create displacement dataset."),
      H5Dclose};
  RequireHdf(H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                      H5P_DEFAULT, flattened.data()),
             "Unable to write displacement dataset.");
  WriteStringAttribute(dataset.Get(), "component_names",
                       options.displacement_components);
  WriteStringAttribute(dataset.Get(), "component_unit_dimensions",
                       "length,length,length,radian,radian,radian");
  WriteStringAttribute(dataset.Get(), "coordinate_system", "global-cartesian");
  WriteStringAttribute(dataset.Get(), "location", "nodal");
  WriteStringAttribute(dataset.Get(), "step_name", "Step-1");
  WriteUint64Attribute(dataset.Get(), "frame_index", 0U);
}

void WriteHdf5(const std::filesystem::path& path,
               const std::filesystem::path& input_path,
               const std::vector<ComparisonValueRow>& values,
               const Hdf5Options& options = {}) {
  Hdf5Handle file{RequireId(H5Fcreate(path.string().c_str(), H5F_ACC_TRUNC,
                                      H5P_DEFAULT, H5P_DEFAULT),
                            "Unable to create comparison HDF5 fixture."),
                  H5Fclose};
  auto metadata = CreateGroup(file.Get(), "/metadata");
  auto model = CreateGroup(file.Get(), "/model");
  auto steps = CreateGroup(file.Get(), "/steps");
  auto step = CreateGroup(file.Get(), "/steps/Step-1");
  auto frames = CreateGroup(file.Get(), "/steps/Step-1/frames");
  auto frame = CreateGroup(file.Get(), "/steps/Step-1/frames/0");
  auto nodal = CreateGroup(file.Get(), "/steps/Step-1/frames/0/nodal");
  (void)model;
  (void)steps;
  (void)step;
  (void)frames;
  (void)frame;
  (void)nodal;

  WriteUint64Attribute(metadata.Get(), "schema_version", 0U);
  WriteStringAttribute(metadata.Get(), "feature_id",
                       "linear-static-mitc4-shell");
  WriteStringAttribute(metadata.Get(), "source_input_identity",
                       "path=" +
                           std::filesystem::absolute(input_path)
                               .lexically_normal()
                               .generic_u8string() +
                           ";content_identity=test");
  WriteStringAttribute(metadata.Get(), "internal_formulation",
                       options.internal_formulation);
  WriteStringAttribute(metadata.Get(), "integration_rule",
                       options.integration_rule);
  WriteStringAttribute(metadata.Get(), "step_name", "Step-1");
  WriteUint64Attribute(metadata.Get(), "frame_index", 0U);
  WriteNodes(file.Get(), values);
  WriteElements(file.Get(), options);
  WriteDisplacement(file.Get(), values, options);
}

void WriteCsv(
    const std::filesystem::path& path,
    const std::vector<ComparisonValueRow>& rows,
    const std::string& header =
        "Part Instance Name,Node Label,U-U1,U-U2,U-U3,UR-UR1,UR-UR2,UR-UR3") {
  std::ofstream stream{path, std::ios::binary | std::ios::trunc};
  if (!stream) {
    throw std::runtime_error{"Unable to create comparison CSV fixture."};
  }
  stream << header << '\n' << std::setprecision(17);
  for (const auto& row : rows) {
    stream << row.instance_name << ',' << row.source_node_label;
    for (const double value : row.values) {
      stream << ',' << value;
    }
    stream << '\n';
  }
}

class ContractFixture {
 public:
  explicit ContractFixture(const std::string& name)
      : root_{std::filesystem::path{FESA_TEST_BINARY_DIR} / "reference" /
              ("mitc4-comparator-" + name)},
        input_{root_ / "case.inp"},
        csv_{root_ / "displacements.csv"},
        results_{root_ / "results.h5"} {
    std::error_code error;
    std::filesystem::remove_all(root_, error);
    error.clear();
    if (!std::filesystem::create_directories(root_, error) || error) {
      throw std::runtime_error{
          "Unable to create comparison fixture directory."};
    }
    std::ofstream input_stream{input_, std::ios::binary | std::ios::trunc};
    input_stream << "*Element, type=S4\n1,1,2,1,2\n";
    input_stream.close();
    WriteCsv(csv_, DefaultRows());
    WriteHdf5(results_, input_, DefaultRows());
  }

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

  fesa::test::Mitc4ReferenceCase ReferenceCase() const {
    return {"shell-contract", "S4", input_, csv_, results_};
  }

  const std::filesystem::path& Root() const noexcept { return root_; }
  const std::filesystem::path& Input() const noexcept { return input_; }
  const std::filesystem::path& Csv() const noexcept { return csv_; }
  const std::filesystem::path& Results() const noexcept { return results_; }

 private:
  std::filesystem::path root_;
  std::filesystem::path input_;
  std::filesystem::path csv_;
  std::filesystem::path results_;
};

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

const fesa::test::Mitc4RowDecision* FindRow(
    const fesa::test::Mitc4ComparisonReport& report, const std::int64_t label,
    const std::string& component) {
  const auto found = std::find_if(report.rows.begin(), report.rows.end(),
                                  [&](const fesa::test::Mitc4RowDecision& row) {
                                    return row.source_node_label == label &&
                                           row.component == component;
                                  });
  return found == report.rows.end() ? nullptr : &*found;
}

const fesa::test::Mitc4ComponentMetrics* FindMetric(
    const fesa::test::Mitc4ComparisonReport& report,
    const std::string& family_identity) {
  const auto found =
      std::find_if(report.metrics.begin(), report.metrics.end(),
                   [&](const fesa::test::Mitc4ComponentMetrics& metric) {
                     return metric.family_identity == family_identity;
                   });
  return found == report.metrics.end() ? nullptr : &*found;
}

std::string ReadBytes(const std::filesystem::path& path) {
  std::ifstream stream{path, std::ios::binary};
  return {std::istreambuf_iterator<char>{stream},
          std::istreambuf_iterator<char>{}};
}

// MITC4-REF-001
TEST(Mitc4ReferenceComparison,
     MapsTrimmedHeaderAndSixComponentsBySourceIdentity) {
  ContractFixture fixture{"mapping"};
  auto csv_rows = DefaultRows();
  std::reverse(csv_rows.begin(), csv_rows.end());
  for (auto& row : csv_rows) {
    row.instance_name = "part-1-1";
  }
  WriteCsv(fixture.Csv(), csv_rows,
           " Part Instance Name , Node Label , U-U1 , U-U2 , U-U3 , "
           "UR-UR1 , UR-UR2 , UR-UR3 ");

  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  const auto& report = result.Value();
  ASSERT_TRUE(report.passed);
  ASSERT_EQ(report.rows.size(), 12U);
  EXPECT_EQ(report.rows[0U].source_node_label, 1);
  EXPECT_EQ(report.rows[0U].component, "U1");
  EXPECT_EQ(report.rows[5U].component, "UR3");
  EXPECT_EQ(report.rows[6U].source_node_label, 2);
  ASSERT_EQ(report.metrics.size(), 2U);
  EXPECT_EQ(report.metrics[0U].row_count, 6U);
  EXPECT_EQ(report.metrics[1U].row_count, 6U);
  EXPECT_TRUE(std::all_of(report.rows.begin(), report.rows.end(),
                          [](const fesa::test::Mitc4RowDecision& row) {
                            return row.case_id == "shell-contract" &&
                                   row.instance_name == kInstanceName &&
                                   row.within_tolerance;
                          }));
}

// MITC4-REF-002
TEST(Mitc4ReferenceComparison, RejectsInvalidInventoryBeforeNumericComparison) {
  {
    ContractFixture fixture{"missing-input"};
    ASSERT_TRUE(std::filesystem::remove(fixture.Input()));
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "needs-reference-artifacts");
  }
  {
    ContractFixture fixture{"header"};
    WriteCsv(fixture.Csv(), DefaultRows(),
             "Part Instance Name,Node Label,U1,U-U2,U-U3,UR-UR1,UR-UR2,UR-UR3");
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"missing-row"};
    auto rows = DefaultRows();
    rows.pop_back();
    WriteCsv(fixture.Csv(), rows);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"extra-row"};
    auto rows = DefaultRows();
    rows.push_back({kInstanceName, 3, {}});
    WriteCsv(fixture.Csv(), rows);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"duplicate-row"};
    auto rows = DefaultRows();
    rows.push_back(rows.front());
    WriteCsv(fixture.Csv(), rows);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"nonfinite-csv"};
    auto rows = DefaultRows();
    rows[0U].values[0U] = std::numeric_limits<double>::quiet_NaN();
    WriteCsv(fixture.Csv(), rows);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"nonfinite-hdf5"};
    auto rows = DefaultRows();
    rows[0U].values[0U] = std::numeric_limits<double>::infinity();
    WriteHdf5(fixture.Results(), fixture.Input(), rows);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"source-identity"};
    auto rows = DefaultRows();
    rows[0U].instance_name = "WRONG-INSTANCE";
    WriteCsv(fixture.Csv(), rows);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
  {
    ContractFixture fixture{"hdf5-schema"};
    Hdf5Options options;
    options.displacement_components = "U1,U2,U3,UR1,UR2,UR3";
    WriteHdf5(fixture.Results(), fixture.Input(), DefaultRows(), options);
    ExpectFailureCode(
        fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()),
        "schema-mismatch");
  }
}

// MITC4-REF-003
TEST(Mitc4ReferenceComparison,
     SharesTranslationScaleAcrossUComponentsAndUsesNearZeroBranches) {
  ContractFixture fixture{"tolerance"};
  auto reference = DefaultRows();
  reference[0U].values[0U] = 0.0;
  reference[1U].values[0U] = 2.0;
  reference[0U].values[1U] = 0.0;
  reference[1U].values[1U] = 0.0;
  WriteCsv(fixture.Csv(), reference);
  auto fesa_values = reference;
  fesa_values[0U].values[0U] = 0.999 * 8.0e-2;
  fesa_values[1U].values[0U] += 1.001e-1;
  fesa_values[0U].values[1U] = 0.999 * 8.0e-2;
  fesa_values[1U].values[1U] = 1.001 * 8.0e-2;
  WriteHdf5(fixture.Results(), fixture.Input(), fesa_values);

  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  const auto& report = result.Value();
  EXPECT_FALSE(report.passed);
  const auto* u1_zero = FindRow(report, 1, "U1");
  const auto* u1_scaled = FindRow(report, 2, "U1");
  const auto* u2_near = FindRow(report, 1, "U2");
  const auto* u2_over = FindRow(report, 2, "U2");
  ASSERT_NE(u1_zero, nullptr);
  ASSERT_NE(u1_scaled, nullptr);
  ASSERT_NE(u2_near, nullptr);
  ASSERT_NE(u2_over, nullptr);
  EXPECT_DOUBLE_EQ(u1_zero->tolerance, 8.0e-2);
  EXPECT_TRUE(u1_zero->within_tolerance);
  EXPECT_FALSE(u1_zero->relative_error_applicable);
  EXPECT_FALSE(u1_scaled->within_tolerance);
  EXPECT_NEAR(u1_scaled->normalized_error, 0.05005, 1.0e-12);
  EXPECT_DOUBLE_EQ(u2_near->tolerance, 8.0e-2);
  EXPECT_TRUE(u2_near->within_tolerance);
  EXPECT_FALSE(u2_near->relative_error_applicable);
  EXPECT_FALSE(u2_over->within_tolerance);
  EXPECT_DOUBLE_EQ(u2_over->tolerance, 8.0e-2);
  const auto* translation = FindMetric(report, "displacement-translation");
  ASSERT_NE(translation, nullptr);
  EXPECT_DOUBLE_EQ(translation->reference_scale, 8.0);
  EXPECT_EQ(translation->components,
            (std::vector<std::string>{"U1", "U2", "U3"}));
}

// MITC4-REF-004
TEST(Mitc4ReferenceComparison,
     RotationExceedanceWarnsWithoutBlockingTranslationVerdict) {
  ContractFixture fixture{"warning"};
  auto fesa_values = DefaultRows();
  fesa_values[0U].values[3U] = 1.0;
  WriteHdf5(fixture.Results(), fixture.Input(), fesa_values);

  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  const auto& report = result.Value();
  EXPECT_TRUE(report.passed);
  ASSERT_EQ(report.warnings.size(), 2U);
  EXPECT_EQ(report.warnings[0U].code, "rotation-reference-exceedance");
  EXPECT_EQ(report.warnings[1U].code, "rotation-reference-rms-exceedance");
  const auto* row = FindRow(report, 1, "UR1");
  ASSERT_NE(row, nullptr);
  EXPECT_FALSE(row->blocking);
  EXPECT_FALSE(row->within_tolerance);
  EXPECT_EQ(report.warnings[0U].row,
            static_cast<std::size_t>(row - report.rows.data()));
  EXPECT_NE(report.warnings[0U].message.find("shell-contract"),
            std::string::npos);
  EXPECT_NE(report.warnings[0U].message.find("UR1"), std::string::npos);
}

TEST(Mitc4ReferenceComparison,
     RotationRmsExceedanceWarnsWhenEveryRotationRowPasses) {
  ContractFixture fixture{"rotation-rms"};
  auto fesa_values = DefaultRows();
  fesa_values[1U].values[3U] = 0.1049;
  fesa_values[1U].values[4U] = -0.2098;
  fesa_values[1U].values[5U] = 0.3147;
  WriteHdf5(fixture.Results(), fixture.Input(), fesa_values);

  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  const auto& report = result.Value();
  EXPECT_TRUE(report.passed);
  ASSERT_EQ(report.warnings.size(), 1U);
  EXPECT_EQ(report.warnings[0U].code, "rotation-reference-rms-exceedance");
  EXPECT_TRUE(std::all_of(report.rows.begin(), report.rows.end(),
                          [](const fesa::test::Mitc4RowDecision& row) {
                            return row.within_tolerance;
                          }));
}

TEST(Mitc4ReferenceComparison,
     RotationRowAndRmsExceedancesProduceSeparateWarnings) {
  ContractFixture fixture{"rotation-row-and-rms"};
  auto fesa_values = DefaultRows();
  fesa_values[0U].values[3U] = 1.0;
  WriteHdf5(fixture.Results(), fixture.Input(), fesa_values);

  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  const auto& report = result.Value();
  EXPECT_TRUE(report.passed);
  ASSERT_EQ(report.warnings.size(), 2U);
  EXPECT_EQ(report.warnings[0U].code, "rotation-reference-exceedance");
  EXPECT_EQ(report.warnings[1U].code, "rotation-reference-rms-exceedance");
  EXPECT_LT(report.warnings[0U].row, report.rows.size());
  EXPECT_LT(report.warnings[1U].row, report.rows.size());
}

// MITC4-REF-005
TEST(Mitc4ReferenceComparison,
     ReportsMetricsVectorsWorstRowAndJsonDeterministically) {
  ContractFixture fixture{"report"};
  auto fesa_values = DefaultRows();
  fesa_values[0U].values[0U] = 0.5e-9;
  fesa_values[0U].values[1U] = -0.25e-9;
  fesa_values[0U].values[3U] = 0.75e-9;
  WriteHdf5(fixture.Results(), fixture.Input(), fesa_values);

  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  const auto& report = result.Value();
  ASSERT_TRUE(report.passed);
  ASSERT_EQ(report.metrics.size(), 2U);
  ASSERT_EQ(report.vector_metrics.size(), 2U);
  EXPECT_NEAR(report.vector_metrics[0U].displacement_norm_error,
              std::sqrt(0.3125) * 1.0e-9, 1.0e-21);
  EXPECT_DOUBLE_EQ(report.vector_metrics[0U].rotation_norm_error, 0.75e-9);
  ASSERT_LT(report.worst_row, report.rows.size());
  EXPECT_EQ(report.rows[report.worst_row].component, "UR1");
  for (const auto& metric : report.metrics) {
    EXPECT_TRUE(std::isfinite(metric.reference_scale));
    EXPECT_TRUE(std::isfinite(metric.near_zero_band));
    EXPECT_TRUE(std::isfinite(metric.maximum_absolute_error));
    EXPECT_TRUE(std::isfinite(metric.relative_rms));
    EXPECT_FALSE(metric.components.empty());
    EXPECT_LT(metric.worst_row, report.rows.size());
  }

  const auto json_a = fixture.Root() / "comparison-a.json";
  const auto json_b = fixture.Root() / "comparison-b.json";
  ASSERT_TRUE(fesa::test::Mitc4ReferenceComparison::WriteDeterministicJson(
                  report, json_a)
                  .IsOk());
  ASSERT_TRUE(fesa::test::Mitc4ReferenceComparison::WriteDeterministicJson(
                  report, json_b)
                  .IsOk());
  const std::string first = ReadBytes(json_a);
  EXPECT_EQ(first, ReadBytes(json_b));
  for (const char* key :
       {"\"rows\"", "\"metrics\"", "\"vector_metrics\"", "\"warnings\"",
        "\"family_identity\"", "\"near_zero_band\"", "\"relative_rms\"",
        "\"row_count\"", "\"worst_row\"", "\"passed\""}) {
    EXPECT_NE(first.find(key), std::string::npos) << key;
  }
}

// MITC4-REF-006
TEST(Mitc4ReferenceComparison, RequiresOnlyDeclaredInputCsvAndHdf5) {
  ContractFixture fixture{"minimal-artifacts"};
  std::vector<std::string> names;
  for (const auto& entry :
       std::filesystem::directory_iterator{fixture.Root()}) {
    names.push_back(entry.path().filename().string());
  }
  std::sort(names.begin(), names.end());
  EXPECT_EQ(names, (std::vector<std::string>{"case.inp", "displacements.csv",
                                             "results.h5"}));
  auto result =
      fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase());
  ASSERT_TRUE(result.HasValue());
  EXPECT_TRUE(result.Value().passed);
}

}  // namespace