From 613bdb9db53cce8f33e96769d9e2153b7b44e215 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 12 Aug 2026 22:18:45 +0900 Subject: [PATCH] feat(linear-static-mitc4-shell): step 13 - shell-reference-comparison --- tests/CMakeLists.txt | 3 + .../reference/mitc4_reference_cases_test.cpp | 194 ++++ .../reference/mitc4_reference_comparison.cpp | 936 ++++++++++++++++++ .../reference/mitc4_reference_comparison.hpp | 81 ++ .../mitc4_reference_comparison_test.cpp | 792 +++++++++++++++ 5 files changed, 2006 insertions(+) create mode 100644 tests/reference/mitc4_reference_cases_test.cpp create mode 100644 tests/reference/mitc4_reference_comparison.cpp create mode 100644 tests/reference/mitc4_reference_comparison.hpp create mode 100644 tests/reference/mitc4_reference_comparison_test.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6cf4378..b404eda 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -59,6 +59,9 @@ add_executable( reference/reference_comparison.cpp reference/reference_comparison_test.cpp reference/b33_reference_comparison_test.cpp + reference/mitc4_reference_comparison.cpp + reference/mitc4_reference_comparison_test.cpp + reference/mitc4_reference_cases_test.cpp ) target_link_libraries( diff --git a/tests/reference/mitc4_reference_cases_test.cpp b/tests/reference/mitc4_reference_cases_test.cpp new file mode 100644 index 0000000..fd38c9f --- /dev/null +++ b/tests/reference/mitc4_reference_cases_test.cpp @@ -0,0 +1,194 @@ +#include "mitc4_reference_comparison.hpp" + +#include "fesa/app/fesa_application.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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* kInternalFormulation = "FESA-MITC4"; +constexpr const char* kIntegrationRule = + "2x2x2-gauss; mitc4-edge-midpoint-shear"; +constexpr std::size_t kNodeCount = 49U; +constexpr std::size_t kComponentCount = 6U; + +std::string readBytes(const std::filesystem::path& path) { + std::ifstream stream{path, std::ios::binary}; + if (!stream) { + throw std::runtime_error{"Unable to read declared reference artifact."}; + } + return {std::istreambuf_iterator{stream}, + std::istreambuf_iterator{}}; +} + +struct ArtifactSnapshot { + std::string bytes; + std::filesystem::file_time_type lastWriteTime; +}; + +ArtifactSnapshot snapshot(const std::filesystem::path& path) { + return {readBytes(path), std::filesystem::last_write_time(path)}; +} + +void expectUnchanged( + const std::filesystem::path& path, const ArtifactSnapshot& before) { + EXPECT_EQ(readBytes(path), before.bytes) << path.string(); + EXPECT_EQ(std::filesystem::last_write_time(path), before.lastWriteTime) + << path.string(); +} + +struct CaseEvidence { + fesa::test::Mitc4ComparisonReport report; + std::filesystem::path comparisonJson; +}; + +CaseEvidence runCase( + const std::string& caseId, + const std::string& sourceElementType, + const std::filesystem::path& referenceDirectory, + const std::filesystem::path& input, + const std::filesystem::path& csv, + const std::string& outputName) { + const auto inputBefore = snapshot(input); + const auto csvBefore = snapshot(csv); + const std::filesystem::path outputDirectory = + std::filesystem::path{FESA_TEST_BINARY_DIR} / "reference" / outputName; + std::error_code error; + std::filesystem::remove_all(outputDirectory, error); + error.clear(); + if (!std::filesystem::create_directories(outputDirectory, error) || error) { + throw std::runtime_error{"Unable to create MITC4 evidence directory."}; + } + const auto results = outputDirectory / "results.h5"; + const auto comparison = outputDirectory / "comparison.json"; + + fesa::FesaApplication application; + EXPECT_EQ( + application.run({input.string(), "--output", results.string()}), 0); + EXPECT_TRUE(std::filesystem::is_regular_file(results)); + auto comparisonResult = fesa::test::Mitc4ReferenceComparison::compare( + {caseId, sourceElementType, input, csv, results}); + if (!comparisonResult.hasValue()) { + ADD_FAILURE() << "MITC4 comparison precheck failed for " << caseId; + return {{}, comparison}; + } + EXPECT_TRUE( + fesa::test::Mitc4ReferenceComparison::writeDeterministicJson( + comparisonResult.value(), comparison) + .isOk()); + EXPECT_TRUE(std::filesystem::is_regular_file(comparison)); + + std::vector generated; + for (const auto& entry : + std::filesystem::directory_iterator{outputDirectory}) { + generated.push_back(entry.path().filename().string()); + } + std::sort(generated.begin(), generated.end()); + EXPECT_EQ( + generated, + (std::vector{"comparison.json", "results.h5"})); + EXPECT_TRUE(std::filesystem::is_directory(referenceDirectory)); + expectUnchanged(input, inputBefore); + expectUnchanged(csv, csvBefore); + return {std::move(comparisonResult.value()), comparison}; +} + +void expectCommonMetadata( + const fesa::test::Mitc4ComparisonReport& report, + const std::string& caseId, + const std::string& sourceElementType) { + EXPECT_EQ(report.caseId, caseId); + EXPECT_EQ(report.sourceElementType, sourceElementType); + EXPECT_EQ(report.internalFormulation, kInternalFormulation); + EXPECT_EQ(report.integrationRule, kIntegrationRule); +} + +void expectComparisonCoverage( + const fesa::test::Mitc4ComparisonReport& report) { + ASSERT_EQ(report.rows.size(), kNodeCount * kComponentCount); + ASSERT_EQ(report.metrics.size(), kComponentCount); + ASSERT_EQ(report.vectorMetrics.size(), kNodeCount); + EXPECT_TRUE(report.passed); + const std::size_t blockingRows = static_cast(std::count_if( + report.rows.begin(), report.rows.end(), + [](const fesa::test::Mitc4RowDecision& row) { + return row.blocking; + })); + const std::size_t rotationRows = report.rows.size() - blockingRows; + EXPECT_EQ(blockingRows, kNodeCount * 3U); + EXPECT_EQ(rotationRows, kNodeCount * 3U); + EXPECT_TRUE(std::all_of( + report.rows.begin(), report.rows.end(), + [](const fesa::test::Mitc4RowDecision& row) { + return !row.blocking || row.withinTolerance; + })); + EXPECT_EQ( + report.warnings.size(), + static_cast(std::count_if( + report.rows.begin(), report.rows.end(), + [](const fesa::test::Mitc4RowDecision& row) { + return !row.blocking && !row.withinTolerance; + }))); +} + +// MITC4-E2E-S4-001 +TEST(Mitc4S4Reference, PreservesS4AndWritesCommonMitc4Metadata) { + const std::filesystem::path root{FESA_TEST_SOURCE_DIR}; + const auto directory = root / "reference" / "shell"; + const auto evidence = runCase( + "shell-s4", "S4", directory, directory / "shell.inp", + directory / "shell displacements.csv", "mitc4-shell-s4-metadata"); + expectCommonMetadata(evidence.report, "shell-s4", "S4"); +} + +// MITC4-E2E-S4-002 +TEST(Mitc4S4Reference, PassesBlockingUAndReportsEveryUrRow) { + const std::filesystem::path root{FESA_TEST_SOURCE_DIR}; + const auto directory = root / "reference" / "shell"; + const auto evidence = runCase( + "shell-s4", "S4", directory, directory / "shell.inp", + directory / "shell displacements.csv", "mitc4-shell-s4-comparison"); + expectCommonMetadata(evidence.report, "shell-s4", "S4"); + expectComparisonCoverage(evidence.report); +} + +// MITC4-E2E-S4R-001 +TEST(Mitc4S4RReference, PreservesS4RAndWritesCommonMitc4Metadata) { + const std::filesystem::path root{FESA_TEST_SOURCE_DIR}; + const auto directory = root / "reference" / "shellR"; + const auto evidence = runCase( + "shell-s4r", "S4R", directory, directory / "shellR.inp", + directory / "shellR displacements.csv", "mitc4-shell-s4r-metadata"); + expectCommonMetadata(evidence.report, "shell-s4r", "S4R"); +} + +// MITC4-E2E-S4R-002 +TEST(Mitc4S4RReference, PassesBlockingUAndReportsEveryUrRow) { + const std::filesystem::path root{FESA_TEST_SOURCE_DIR}; + const auto directory = root / "reference" / "shellR"; + const auto evidence = runCase( + "shell-s4r", "S4R", directory, directory / "shellR.inp", + directory / "shellR displacements.csv", "mitc4-shell-s4r-comparison"); + expectCommonMetadata(evidence.report, "shell-s4r", "S4R"); + expectComparisonCoverage(evidence.report); +} + +} // namespace diff --git a/tests/reference/mitc4_reference_comparison.cpp b/tests/reference/mitc4_reference_comparison.cpp new file mode 100644 index 0000000..796a28a --- /dev/null +++ b/tests/reference/mitc4_reference_comparison.cpp @@ -0,0 +1,936 @@ +#include "mitc4_reference_comparison.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa::test { +namespace { + +constexpr const char* kDisplacementPath = + "/steps/Step-1/frames/0/nodal/displacement"; +constexpr const char* kInternalFormulation = "FESA-MITC4"; +constexpr const char* kIntegrationRule = + "2x2x2-gauss; mitc4-edge-midpoint-shear"; +constexpr double kAbsoluteFloor = 1.0e-9; +constexpr double kRelativeCoefficient = 1.0e-6; +constexpr std::array kComponents{ + "U1", "U2", "U3", "UR1", "UR2", "UR3"}; +const std::vector kExpectedHeader{ + "Part Instance Name", "Node Label", "U-U1", "U-U2", "U-U3", + "UR-UR1", "UR-UR2", "UR-UR3"}; + +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 failureStatus( + const std::string& caseId, + const std::string& code, + const std::string& message) { + return Status::failure( + FailureCategory::model, + {{Severity::error, code, {}, "", caseId, 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(character)); + }); + const auto end = std::find_if_not( + value.rbegin(), value.rend(), [&](const char character) { + return isSpace(static_cast(character)); + }).base(); + return begin < end ? std::string{begin, end} : std::string{}; +} + +std::vector splitCsvLine(const std::string& line) { + std::vector 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 source-node label is invalid."); + } + return value; +} + +double parseFiniteDouble(const std::string& field) { + if (field.empty()) { + fail("schema-mismatch", "A displacement CSV 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 displacement CSV numeric field is invalid or nonfinite."); + } + return value; +} + +struct IdentityKey { + std::string instanceName; + std::int64_t sourceNodeLabel; + + bool operator<(const IdentityKey& other) const noexcept { + if (instanceName != other.instanceName) { + return instanceName < other.instanceName; + } + return sourceNodeLabel < other.sourceNodeLabel; + } +}; + +struct WideRow { + IdentityKey identity; + std::array values; +}; + +std::vector readReferenceCsv(const std::filesystem::path& path) { + std::ifstream stream{path}; + if (!stream) { + fail( + "needs-reference-artifacts", + "The declared displacement CSV is missing or unreadable."); + } + std::string line; + if (!std::getline(stream, line)) { + fail("schema-mismatch", "The declared displacement CSV is empty."); + } + if (!line.empty() && line.back() == '\r') { + line.pop_back(); + } + if (splitCsvLine(line) != kExpectedHeader) { + fail( + "schema-mismatch", + "The displacement CSV header does not match the six-component contract."); + } + + std::vector rows; + std::map identities; + while (std::getline(stream, line)) { + if (!line.empty() && line.back() == '\r') { + line.pop_back(); + } + if (line.empty()) { + fail("schema-mismatch", "Blank displacement CSV rows are not allowed."); + } + const auto fields = splitCsvLine(line); + if (fields.size() != kExpectedHeader.size() || fields[0U].empty()) { + fail("schema-mismatch", "A displacement CSV row has invalid schema."); + } + WideRow row{}; + row.identity.instanceName = fields[0U]; + row.identity.sourceNodeLabel = parsePositiveLabel(fields[1U]); + for (std::size_t component = 0U; + component < row.values.size(); + ++component) { + row.values[component] = parseFiniteDouble(fields[component + 2U]); + } + if (!identities.emplace(row.identity, rows.size()).second) { + fail("schema-mismatch", "A displacement CSV row identity is duplicated."); + } + rows.push_back(std::move(row)); + } + if (rows.empty()) { + fail("schema-mismatch", "The displacement CSV contains no data rows."); + } + return rows; +} + +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() { + if (value_ >= 0 && closer_ != nullptr) { + (void)closer_(value_); + } + } + + hid_t get() const noexcept { return value_; } + +private: + 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_); + } + } + + void reclaim() { + active_ = false; + if (H5Dvlen_reclaim(memoryType_, dataSpace_, H5P_DEFAULT, data_) < 0) { + fail("schema-mismatch", "Unable to reclaim HDF5 variable strings."); + } + } + +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 makeUtf8StringType() { + Hdf5Handle type{ + requireId(H5Tcopy(H5T_C_S1), "Unable to copy an HDF5 string type."), + H5Tclose}; + requireHdf( + H5Tset_size(type.get(), H5T_VARIABLE), + "Unable to define an HDF5 variable string type."); + requireHdf( + H5Tset_cset(type.get(), H5T_CSET_UTF8), + "Unable to define an HDF5 UTF-8 string type."); + return type; +} + +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 dimensions(const hid_t dataset) { + Hdf5Handle space{ + requireId(H5Dget_space(dataset), "Unable to inspect HDF5 dimensions."), + H5Sclose}; + const int rank = H5Sget_simple_extent_ndims(space.get()); + if (rank < 0) { + fail("schema-mismatch", "Unable to inspect HDF5 rank."); + } + std::vector result(static_cast(rank)); + if (rank > 0) { + requireHdf( + H5Sget_simple_extent_dims(space.get(), result.data(), nullptr), + "Unable to inspect HDF5 extents."); + } + return result; +} + +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 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 result{raw}; + requireHdf(H5free_memory(raw), "Unable to free HDF5 attribute memory."); + return result; +} + +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 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 result = 0U; + requireHdf( + H5Aread(attribute.get(), H5T_NATIVE_UINT64, &result), + "Unable to read an HDF5 integer attribute."); + return result; +} + +void requireStringAttribute( + const hid_t object, const char* name, const std::string& expected) { + if (readStringAttribute(object, name) != expected) { + fail("schema-mismatch", "An HDF5 string attribute has the wrong value."); + } +} + +void requireExactCompoundMembers( + const hid_t dataset, const std::vector& 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(expected.size())) { + fail("schema-mismatch", "An HDF5 compound schema has the wrong member count."); + } + for (std::size_t index = 0U; index < expected.size(); ++index) { + char* name = H5Tget_member_name(type.get(), static_cast(index)); + if (name == nullptr) { + fail("schema-mismatch", "Unable to inspect an HDF5 compound member."); + } + const std::string actual{name}; + requireHdf(H5free_memory(name), "Unable to free HDF5 member memory."); + if (actual != expected[index]) { + fail("schema-mismatch", "An HDF5 compound member is out of contract order."); + } + } +} + +struct NodeReadRow { + std::uint64_t internalNodeId; + char* instanceName; + char* sourceLabel; + double coordinates[3]; +}; + +std::vector readNodes(const hid_t file) { + auto dataset = openDataset(file, "/model/nodes"); + const auto shape = dimensions(dataset.get()); + if (shape.size() != 1U || shape[0U] == 0U || + shape[0U] > static_cast((std::numeric_limits::max)())) { + fail("schema-mismatch", "The HDF5 node dataset has an invalid shape."); + } + requireExactCompoundMembers( + 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 = makeUtf8StringType(); + const hsize_t coordinateDimensions[] = {3U}; + Hdf5Handle coordinates{ + requireId( + H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinateDimensions), + "Unable to create node coordinate memory type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireId( + H5Tcreate(H5T_COMPOUND, sizeof(NodeReadRow)), + "Unable to create node memory type."), + H5Tclose}; + requireHdf( + H5Tinsert( + memoryType.get(), "internal_node_id", + HOFFSET(NodeReadRow, internalNodeId), H5T_NATIVE_UINT64), + "Unable to define node ID memory field."); + requireHdf( + H5Tinsert( + memoryType.get(), "instance_name", HOFFSET(NodeReadRow, instanceName), + stringType.get()), + "Unable to define node instance memory field."); + requireHdf( + H5Tinsert( + memoryType.get(), "source_label", HOFFSET(NodeReadRow, sourceLabel), + stringType.get()), + "Unable to define node label memory field."); + requireHdf( + H5Tinsert( + memoryType.get(), "coordinates", HOFFSET(NodeReadRow, coordinates), + coordinates.get()), + "Unable to define node coordinate memory field."); + + const std::size_t count = static_cast(shape[0U]); + std::vector raw(count); + requireHdf( + H5Dread( + dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, + raw.data()), + "Unable to read HDF5 node identities."); + Hdf5Handle space{ + requireId(H5Dget_space(dataset.get()), "Unable to inspect node space."), + H5Sclose}; + Hdf5VlenReclaimer reclaimer{memoryType.get(), space.get(), raw.data()}; + std::vector rows; + rows.reserve(count); + std::map identities; + for (std::size_t index = 0U; index < count; ++index) { + if (raw[index].internalNodeId != index || raw[index].instanceName == nullptr || + raw[index].sourceLabel == nullptr || raw[index].instanceName[0] == '\0' || + !std::all_of( + std::begin(raw[index].coordinates), + std::end(raw[index].coordinates), + [](const double value) { return std::isfinite(value); })) { + fail("schema-mismatch", "An HDF5 node row has invalid identity or values."); + } + WideRow row{}; + row.identity.instanceName = raw[index].instanceName; + row.identity.sourceNodeLabel = parsePositiveLabel(raw[index].sourceLabel); + if (!identities.emplace(row.identity, rows.size()).second) { + fail("schema-mismatch", "An HDF5 node identity is duplicated."); + } + rows.push_back(std::move(row)); + } + reclaimer.reclaim(); + return rows; +} + +struct ElementIdentityReadRow { + char* sourceElementType; + char* internalFormulation; +}; + +void requireElementIdentity( + const hid_t file, const std::string& expectedSourceType) { + auto dataset = openDataset(file, "/model/elements"); + const auto shape = dimensions(dataset.get()); + if (shape.size() != 1U || shape[0U] == 0U || + shape[0U] > static_cast((std::numeric_limits::max)())) { + fail("schema-mismatch", "The HDF5 element dataset has an invalid shape."); + } + requireExactCompoundMembers( + dataset.get(), + {"internal_element_id", "instance_name", "source_label", + "source_element_type", "internal_formulation", "node_internal_ids", + "shell_section_internal_id", "material_internal_id"}); + requireStringAttribute(dataset.get(), "formulation", kInternalFormulation); + + auto stringType = makeUtf8StringType(); + Hdf5Handle memoryType{ + requireId( + H5Tcreate(H5T_COMPOUND, sizeof(ElementIdentityReadRow)), + "Unable to create element identity memory type."), + H5Tclose}; + requireHdf( + H5Tinsert( + memoryType.get(), "source_element_type", + HOFFSET(ElementIdentityReadRow, sourceElementType), stringType.get()), + "Unable to define source element type memory field."); + requireHdf( + H5Tinsert( + memoryType.get(), "internal_formulation", + HOFFSET(ElementIdentityReadRow, internalFormulation), stringType.get()), + "Unable to define formulation memory field."); + const std::size_t count = static_cast(shape[0U]); + std::vector rows(count); + requireHdf( + H5Dread( + dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL, H5P_DEFAULT, + rows.data()), + "Unable to read HDF5 element identity."); + Hdf5Handle space{ + requireId(H5Dget_space(dataset.get()), "Unable to inspect element space."), + H5Sclose}; + Hdf5VlenReclaimer reclaimer{memoryType.get(), space.get(), rows.data()}; + for (const auto& row : rows) { + if (row.sourceElementType == nullptr || row.internalFormulation == nullptr || + row.sourceElementType != expectedSourceType || + row.internalFormulation != std::string{kInternalFormulation}) { + fail( + "schema-mismatch", + "The HDF5 source element type or internal formulation is invalid."); + } + } + reclaimer.reclaim(); +} + +std::vector readDisplacement( + const hid_t file, const std::size_t nodeCount) { + auto dataset = openDataset(file, kDisplacementPath); + if (dimensions(dataset.get()) != + std::vector{static_cast(nodeCount), 6U}) { + fail("schema-mismatch", "The HDF5 displacement dataset has the wrong shape."); + } + Hdf5Handle type{ + requireId(H5Dget_type(dataset.get()), "Unable to inspect displacement 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", "The HDF5 displacement dataset is not float64 LE."); + } + requireStringAttribute(dataset.get(), "component_names", "UX,UY,UZ,URX,URY,URZ"); + requireStringAttribute( + dataset.get(), "component_unit_dimensions", + "length,length,length,radian,radian,radian"); + requireStringAttribute(dataset.get(), "coordinate_system", "global-cartesian"); + requireStringAttribute(dataset.get(), "location", "nodal"); + requireStringAttribute(dataset.get(), "step_name", "Step-1"); + if (readUint64Attribute(dataset.get(), "frame_index") != 0U) { + fail("schema-mismatch", "The HDF5 displacement frame identity is invalid."); + } + std::vector values(nodeCount * kComponents.size()); + if (!values.empty()) { + requireHdf( + H5Dread( + dataset.get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), + "Unable to read HDF5 displacement values."); + } + if (!std::all_of(values.begin(), values.end(), [](const double value) { + return std::isfinite(value); + })) { + fail("schema-mismatch", "An HDF5 displacement value is nonfinite."); + } + return values; +} + +struct Hdf5Projection { + std::vector rows; + std::string sourceElementType; + std::string internalFormulation; + std::string integrationRule; +}; + +Hdf5Projection readHdf5(const Mitc4ReferenceCase& referenceCase) { + Hdf5ErrorSilencer silencer; + Hdf5Handle file{ + requireId( + H5Fopen( + referenceCase.resultsHdf5Path.string().c_str(), H5F_ACC_RDONLY, + H5P_DEFAULT), + "The authoritative HDF5 output cannot be opened."), + 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") != 0U) { + fail("schema-mismatch", "The HDF5 schema or frame version is invalid."); + } + requireStringAttribute(metadata.get(), "feature_id", "linear-static-mitc4-shell"); + requireStringAttribute(metadata.get(), "step_name", "Step-1"); + requireStringAttribute( + metadata.get(), "internal_formulation", kInternalFormulation); + requireStringAttribute(metadata.get(), "integration_rule", kIntegrationRule); + const std::string normalizedInput = + std::filesystem::absolute(referenceCase.inputPath) + .lexically_normal() + .generic_u8string(); + const std::string sourceIdentity = + readStringAttribute(metadata.get(), "source_input_identity"); + if (sourceIdentity.find("path=" + normalizedInput + ";content_identity=") != 0U) { + fail( + "schema-mismatch", + "The HDF5 source input identity does not match the declared input."); + } + auto rows = readNodes(file.get()); + requireElementIdentity(file.get(), referenceCase.expectedSourceElementType); + const auto values = readDisplacement(file.get(), rows.size()); + for (std::size_t row = 0U; row < rows.size(); ++row) { + std::copy_n( + values.begin() + static_cast(row * kComponents.size()), + kComponents.size(), rows[row].values.begin()); + } + return { + std::move(rows), referenceCase.expectedSourceElementType, + kInternalFormulation, kIntegrationRule}; +} + +void requireArtifacts(const Mitc4ReferenceCase& referenceCase) { + if (referenceCase.caseId.empty() || + (referenceCase.expectedSourceElementType != "S4" && + referenceCase.expectedSourceElementType != "S4R")) { + fail("schema-mismatch", "The MITC4 reference case identity is invalid."); + } + std::error_code error; + for (const auto* path : { + &referenceCase.inputPath, + &referenceCase.displacementCsvPath, + &referenceCase.resultsHdf5Path}) { + if (!std::filesystem::is_regular_file(*path, error) || error) { + fail( + "needs-reference-artifacts", + "A declared MITC4 input, displacement CSV, or HDF5 file is missing."); + } + } +} + +std::string finiteText(const double value) { + std::ostringstream stream; + stream.imbue(std::locale::classic()); + stream << std::setprecision(std::numeric_limits::max_digits10) + << value; + return stream.str(); +} + +std::string jsonEscape(const std::string& value) { + std::ostringstream stream; + 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" << std::hex << std::setw(2) + << std::setfill('0') << static_cast(character) + << std::dec << std::setfill(' '); + } else { + stream << static_cast(character); + } + break; + } + } + return stream.str(); +} + +} // namespace + +Result Mitc4ReferenceComparison::compare( + const Mitc4ReferenceCase& referenceCase) { + try { + requireArtifacts(referenceCase); + const auto referenceRows = + readReferenceCsv(referenceCase.displacementCsvPath); + auto hdf5 = readHdf5(referenceCase); + + std::map referenceByIdentity; + for (const auto& row : referenceRows) { + referenceByIdentity.emplace(row.identity, &row); + } + if (referenceRows.size() != hdf5.rows.size()) { + fail( + "schema-mismatch", + "The HDF5 and CSV source-node row counts are not equal."); + } + for (const auto& row : hdf5.rows) { + if (referenceByIdentity.find(row.identity) == referenceByIdentity.end()) { + fail( + "schema-mismatch", + "The HDF5 and CSV source-node identities are not equal."); + } + } + + std::array scales{}; + for (const auto& row : referenceRows) { + for (std::size_t component = 0U; + component < scales.size(); + ++component) { + scales[component] = + (std::max)(scales[component], std::abs(row.values[component])); + } + } + + Mitc4ComparisonReport report{}; + report.caseId = referenceCase.caseId; + report.sourceElementType = std::move(hdf5.sourceElementType); + report.internalFormulation = std::move(hdf5.internalFormulation); + report.integrationRule = std::move(hdf5.integrationRule); + report.passed = true; + report.rows.reserve(hdf5.rows.size() * kComponents.size()); + report.vectorMetrics.reserve(hdf5.rows.size()); + std::array errorNorms{}; + std::array maximumErrors{}; + std::array maximumNormalized{}; + std::array worstRows{}; + double globalWorstNormalized = -1.0; + + for (const auto& hdf5Row : hdf5.rows) { + const auto reference = referenceByIdentity.find(hdf5Row.identity); + if (reference == referenceByIdentity.end()) { + fail("schema-mismatch", "A projected reference row is missing."); + } + std::array errors{}; + for (std::size_t component = 0U; + component < kComponents.size(); + ++component) { + const double tolerance = + kAbsoluteFloor + kRelativeCoefficient * scales[component]; + const double absoluteError = std::abs( + hdf5Row.values[component] - + reference->second->values[component]); + const double normalizedError = absoluteError / tolerance; + if (!std::isfinite(tolerance) || !(tolerance > 0.0) || + !std::isfinite(absoluteError) || + !std::isfinite(normalizedError)) { + fail( + "schema-mismatch", + "A finite row produced a nonfinite comparison metric."); + } + const bool blocking = component < 3U; + const bool withinTolerance = absoluteError <= tolerance; + const std::size_t rowIndex = report.rows.size(); + report.rows.push_back({ + referenceCase.caseId, + hdf5Row.identity.instanceName, + hdf5Row.identity.sourceNodeLabel, + kComponents[component], + hdf5Row.values[component], + reference->second->values[component], + absoluteError, + tolerance, + normalizedError, + blocking, + withinTolerance}); + errors[component] = absoluteError; + errorNorms[component] = + std::hypot(errorNorms[component], absoluteError); + if (normalizedError > maximumNormalized[component]) { + maximumNormalized[component] = normalizedError; + maximumErrors[component] = absoluteError; + worstRows[component] = rowIndex; + } + if (normalizedError > globalWorstNormalized) { + globalWorstNormalized = normalizedError; + report.worstRow = rowIndex; + } + if (blocking && !withinTolerance) { + report.passed = false; + } else if (!blocking && !withinTolerance) { + const std::string message = + "case=" + referenceCase.caseId + ";instance=" + + hdf5Row.identity.instanceName + ";node=" + + std::to_string(hdf5Row.identity.sourceNodeLabel) + + ";component=" + kComponents[component] + + ";absolute_error=" + finiteText(absoluteError) + + ";tolerance=" + finiteText(tolerance); + report.warnings.push_back({ + "rotation-reference-exceedance", rowIndex, message}); + } + } + report.vectorMetrics.push_back({ + hdf5Row.identity.instanceName, + hdf5Row.identity.sourceNodeLabel, + std::hypot(errors[0U], errors[1U], errors[2U]), + std::hypot(errors[3U], errors[4U], errors[5U])}); + } + + report.metrics.reserve(kComponents.size()); + const double rowCount = static_cast(hdf5.rows.size()); + for (std::size_t component = 0U; + component < kComponents.size(); + ++component) { + report.metrics.push_back({ + kComponents[component], + scales[component], + kAbsoluteFloor + kRelativeCoefficient * scales[component], + maximumErrors[component], + maximumNormalized[component], + errorNorms[component] / std::sqrt(rowCount), + errorNorms[component], + worstRows[component]}); + } + return Result::success(std::move(report)); + } catch (const ComparisonFailure& exception) { + return Result::failure(failureStatus( + referenceCase.caseId, exception.code(), exception.what())); + } catch (const std::exception& exception) { + return Result::failure(failureStatus( + referenceCase.caseId, "comparison-failure", exception.what())); + } +} + +Status Mitc4ReferenceComparison::writeDeterministicJson( + const Mitc4ComparisonReport& report, + const std::filesystem::path& outputJson) { + try { + std::ofstream stream{outputJson, std::ios::binary | std::ios::trunc}; + if (!stream) { + return failureStatus( + report.caseId, + "comparison-report-write-failed", + "The deterministic MITC4 JSON report cannot be opened."); + } + stream.imbue(std::locale::classic()); + stream << std::setprecision(std::numeric_limits::max_digits10); + stream << "{\"case_id\":\"" << jsonEscape(report.caseId) + << "\",\"source_element_type\":\"" + << jsonEscape(report.sourceElementType) + << "\",\"internal_formulation\":\"" + << jsonEscape(report.internalFormulation) + << "\",\"integration_rule\":\"" + << jsonEscape(report.integrationRule) << "\",\"rows\":["; + for (std::size_t index = 0U; index < report.rows.size(); ++index) { + if (index != 0U) { + stream << ','; + } + const auto& row = report.rows[index]; + stream << "{\"case_id\":\"" << jsonEscape(row.caseId) + << "\",\"instance_name\":\"" + << jsonEscape(row.instanceName) + << "\",\"source_node_label\":" << row.sourceNodeLabel + << ",\"component\":\"" << jsonEscape(row.component) + << "\",\"fesa_value\":" << row.fesaValue + << ",\"reference_value\":" << row.referenceValue + << ",\"absolute_error\":" << row.absoluteError + << ",\"tolerance\":" << row.tolerance + << ",\"normalized_error\":" << row.normalizedError + << ",\"blocking\":" << (row.blocking ? "true" : "false") + << ",\"within_tolerance\":" + << (row.withinTolerance ? "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 << "{\"component\":\"" << jsonEscape(metric.component) + << "\",\"reference_scale\":" << metric.referenceScale + << ",\"tolerance\":" << metric.tolerance + << ",\"maximum_absolute_error\":" + << metric.maximumAbsoluteError + << ",\"maximum_normalized_error\":" + << metric.maximumNormalizedError + << ",\"rms_error\":" << metric.rmsError + << ",\"vector_norm_error\":" << metric.vectorNormError + << ",\"worst_row\":" << metric.worstRow << '}'; + } + stream << "],\"vector_metrics\":["; + for (std::size_t index = 0U; index < report.vectorMetrics.size(); ++index) { + if (index != 0U) { + stream << ','; + } + const auto& metric = report.vectorMetrics[index]; + stream << "{\"instance_name\":\"" + << jsonEscape(metric.instanceName) + << "\",\"source_node_label\":" << metric.sourceNodeLabel + << ",\"displacement_norm_error\":" + << metric.displacementNormError + << ",\"rotation_norm_error\":" + << metric.rotationNormError << '}'; + } + stream << "],\"warnings\":["; + for (std::size_t index = 0U; index < report.warnings.size(); ++index) { + if (index != 0U) { + stream << ','; + } + const auto& warning = report.warnings[index]; + stream << "{\"code\":\"" << jsonEscape(warning.code) + << "\",\"row\":" << warning.row + << ",\"message\":\"" << jsonEscape(warning.message) + << "\"}"; + } + stream << "],\"worst_row\":" << report.worstRow + << ",\"passed\":" << (report.passed ? "true" : "false") + << "}\n"; + stream.flush(); + if (!stream) { + return failureStatus( + report.caseId, + "comparison-report-write-failed", + "The deterministic MITC4 JSON report could not be completed."); + } + return Status::ok(); + } catch (const std::exception& exception) { + return failureStatus( + report.caseId, "comparison-report-write-failed", exception.what()); + } +} + +} // namespace fesa::test diff --git a/tests/reference/mitc4_reference_comparison.hpp b/tests/reference/mitc4_reference_comparison.hpp new file mode 100644 index 0000000..9a77bf3 --- /dev/null +++ b/tests/reference/mitc4_reference_comparison.hpp @@ -0,0 +1,81 @@ +#pragma once + +#include "fesa/core/status.hpp" + +#include +#include +#include +#include +#include + +namespace fesa::test { + +struct Mitc4ReferenceCase { + std::string caseId; + std::string expectedSourceElementType; + std::filesystem::path inputPath; + std::filesystem::path displacementCsvPath; + std::filesystem::path resultsHdf5Path; +}; + +struct Mitc4RowDecision { + std::string caseId; + std::string instanceName; + std::int64_t sourceNodeLabel; + std::string component; + double fesaValue; + double referenceValue; + double absoluteError; + double tolerance; + double normalizedError; + bool blocking; + bool withinTolerance; +}; + +struct Mitc4ComponentMetrics { + std::string component; + double referenceScale; + double tolerance; + double maximumAbsoluteError; + double maximumNormalizedError; + double rmsError; + double vectorNormError; + std::size_t worstRow; +}; + +struct Mitc4VectorMetrics { + std::string instanceName; + std::int64_t sourceNodeLabel; + double displacementNormError; + double rotationNormError; +}; + +struct Mitc4Warning { + std::string code; + std::size_t row; + std::string message; +}; + +struct Mitc4ComparisonReport { + std::string caseId; + std::string sourceElementType; + std::string internalFormulation; + std::string integrationRule; + std::vector rows; + std::vector metrics; + std::vector vectorMetrics; + std::vector warnings; + std::size_t worstRow; + bool passed; +}; + +class Mitc4ReferenceComparison { +public: + static Result compare( + const Mitc4ReferenceCase& referenceCase); + static Status writeDeterministicJson( + const Mitc4ComparisonReport& report, + const std::filesystem::path& outputJson); +}; + +} // namespace fesa::test diff --git a/tests/reference/mitc4_reference_comparison_test.cpp b/tests/reference/mitc4_reference_comparison_test.cpp new file mode 100644 index 0000000..038c759 --- /dev/null +++ b/tests/reference/mitc4_reference_comparison_test.cpp @@ -0,0 +1,792 @@ +#include "mitc4_reference_comparison.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 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 instanceName{kInstanceName}; + std::int64_t sourceNodeLabel{}; + std::array values{}; +}; + +std::vector 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 sourceElementType{"S4"}; + std::string internalFormulation{kInternalFormulation}; + std::string integrationRule{kIntegrationRule}; + std::string displacementComponents{"UX,UY,UZ,URX,URY,URZ"}; +}; + +struct NodeWriteRow { + std::uint64_t internalNodeId; + const char* instanceName; + const char* sourceLabel; + double coordinates[3]; +}; + +struct ElementWriteRow { + std::uint64_t internalElementId; + const char* instanceName; + const char* sourceLabel; + const char* sourceElementType; + const char* internalFormulation; + std::uint64_t nodeInternalIds[4]; + std::uint64_t shellSectionInternalId; + std::uint64_t materialInternalId; +}; + +void writeNodes( + const hid_t file, const std::vector& values) { + std::vector labels; + labels.reserve(values.size()); + for (const auto& row : values) { + labels.push_back(std::to_string(row.sourceNodeLabel)); + } + std::vector rows; + rows.reserve(values.size()); + for (std::size_t index = 0U; index < values.size(); ++index) { + rows.push_back({ + static_cast(index), + values[index].instanceName.c_str(), + labels[index].c_str(), + {static_cast(index), 0.0, 0.0}}); + } + + auto stringType = makeUtf8StringType(); + const hsize_t coordinateDimensions[] = {3U}; + Hdf5Handle fileCoordinates{ + requireId( + H5Tarray_create2(H5T_IEEE_F64LE, 1, coordinateDimensions), + "Unable to create coordinate file type."), + H5Tclose}; + Hdf5Handle memoryCoordinates{ + requireId( + H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinateDimensions), + "Unable to create coordinate memory type."), + H5Tclose}; + Hdf5Handle fileType{ + requireId( + H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), + "Unable to create node file type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireId( + H5Tcreate(H5T_COMPOUND, sizeof(NodeWriteRow)), + "Unable to create node memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, + const hid_t integerType, + const hid_t coordinateType) { + requireHdf( + H5Tinsert( + type, "internal_node_id", + HOFFSET(NodeWriteRow, internalNodeId), integerType), + "Unable to define node ID."); + requireHdf( + H5Tinsert( + type, "instance_name", HOFFSET(NodeWriteRow, instanceName), + stringType.get()), + "Unable to define node instance."); + requireHdf( + H5Tinsert( + type, "source_label", HOFFSET(NodeWriteRow, sourceLabel), + stringType.get()), + "Unable to define node label."); + requireHdf( + H5Tinsert( + type, "coordinates", HOFFSET(NodeWriteRow, coordinates), + coordinateType), + "Unable to define node coordinates."); + }; + insert(fileType.get(), H5T_STD_U64LE, fileCoordinates.get()); + insert(memoryType.get(), H5T_NATIVE_UINT64, memoryCoordinates.get()); + const hsize_t dimensions[] = {static_cast(rows.size())}; + Hdf5Handle space{ + requireId( + H5Screate_simple(1, dimensions, nullptr), + "Unable to create node space."), + H5Sclose}; + Hdf5Handle dataset{ + requireId( + H5Dcreate2( + file, "/model/nodes", fileType.get(), space.get(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + "Unable to create node dataset."), + H5Dclose}; + requireHdf( + H5Dwrite( + dataset.get(), memoryType.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.sourceElementType.c_str(), + options.internalFormulation.c_str(), + {0U, 1U, 0U, 1U}, + 0U, + 0U}; + auto stringType = makeUtf8StringType(); + const hsize_t nodeDimensions[] = {4U}; + Hdf5Handle fileNodes{ + requireId( + H5Tarray_create2(H5T_STD_U64LE, 1, nodeDimensions), + "Unable to create element node file type."), + H5Tclose}; + Hdf5Handle memoryNodes{ + requireId( + H5Tarray_create2(H5T_NATIVE_UINT64, 1, nodeDimensions), + "Unable to create element node memory type."), + H5Tclose}; + Hdf5Handle fileType{ + requireId( + H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), + "Unable to create element file type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireId( + H5Tcreate(H5T_COMPOUND, sizeof(ElementWriteRow)), + "Unable to create element memory type."), + H5Tclose}; + const auto insert = [&](const hid_t type, + const hid_t integerType, + const hid_t nodesType) { + requireHdf( + H5Tinsert( + type, "internal_element_id", + HOFFSET(ElementWriteRow, internalElementId), integerType), + "Unable to define element ID."); + requireHdf( + H5Tinsert( + type, "instance_name", HOFFSET(ElementWriteRow, instanceName), + stringType.get()), + "Unable to define element instance."); + requireHdf( + H5Tinsert( + type, "source_label", HOFFSET(ElementWriteRow, sourceLabel), + stringType.get()), + "Unable to define element label."); + requireHdf( + H5Tinsert( + type, "source_element_type", + HOFFSET(ElementWriteRow, sourceElementType), stringType.get()), + "Unable to define source element type."); + requireHdf( + H5Tinsert( + type, "internal_formulation", + HOFFSET(ElementWriteRow, internalFormulation), stringType.get()), + "Unable to define internal formulation."); + requireHdf( + H5Tinsert( + type, "node_internal_ids", + HOFFSET(ElementWriteRow, nodeInternalIds), nodesType), + "Unable to define element connectivity."); + requireHdf( + H5Tinsert( + type, "shell_section_internal_id", + HOFFSET(ElementWriteRow, shellSectionInternalId), integerType), + "Unable to define section ID."); + requireHdf( + H5Tinsert( + type, "material_internal_id", + HOFFSET(ElementWriteRow, materialInternalId), integerType), + "Unable to define material ID."); + }; + insert(fileType.get(), H5T_STD_U64LE, fileNodes.get()); + insert(memoryType.get(), H5T_NATIVE_UINT64, memoryNodes.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", fileType.get(), space.get(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), + "Unable to create element dataset."), + H5Dclose}; + requireHdf( + H5Dwrite( + dataset.get(), memoryType.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& values, + const Hdf5Options& options) { + std::vector 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(values.size()), + static_cast(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.displacementComponents); + 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& inputPath, + const std::vector& 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(inputPath) + .lexically_normal() + .generic_u8string() + + ";content_identity=test"); + writeStringAttribute( + metadata.get(), "internal_formulation", options.internalFormulation); + writeStringAttribute( + metadata.get(), "integration_rule", options.integrationRule); + 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& 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.instanceName << ',' << row.sourceNodeLabel; + 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 inputStream{input_, std::ios::binary | std::ios::trunc}; + inputStream << "*Element, type=S4\n1,1,2,1,2\n"; + inputStream.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& result, + const std::string& code) { + ASSERT_FALSE(result.hasValue()); + ASSERT_FALSE(result.status().diagnostics().empty()); + EXPECT_EQ(result.status().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.sourceNodeLabel == label && row.component == component; + }); + return found == report.rows.end() ? nullptr : &*found; +} + +const fesa::test::Mitc4ComponentMetrics* findMetric( + const fesa::test::Mitc4ComparisonReport& report, + const std::string& component) { + const auto found = std::find_if( + report.metrics.begin(), report.metrics.end(), + [&](const fesa::test::Mitc4ComponentMetrics& metric) { + return metric.component == component; + }); + 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{stream}, + std::istreambuf_iterator{}}; +} + +// MITC4-REF-001 +TEST(Mitc4ReferenceComparison, MapsTrimmedHeaderAndSixComponentsBySourceIdentity) { + ContractFixture fixture{"mapping"}; + auto csvRows = defaultRows(); + std::reverse(csvRows.begin(), csvRows.end()); + writeCsv( + fixture.csv(), csvRows, + " 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].sourceNodeLabel, 1); + EXPECT_EQ(report.rows[0U].component, "U1"); + EXPECT_EQ(report.rows[5U].component, "UR3"); + EXPECT_EQ(report.rows[6U].sourceNodeLabel, 2); + EXPECT_TRUE(std::all_of( + report.rows.begin(), report.rows.end(), + [](const fesa::test::Mitc4RowDecision& row) { + return row.caseId == "shell-contract" && + row.instanceName == kInstanceName && row.withinTolerance; + })); +} + +// 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::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::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].instanceName = "WRONG-INSTANCE"; + writeCsv(fixture.csv(), rows); + expectFailureCode( + fesa::test::Mitc4ReferenceComparison::compare( + fixture.referenceCase()), + "schema-mismatch"); + } + { + ContractFixture fixture{"hdf5-schema"}; + Hdf5Options options; + options.displacementComponents = "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, AppliesAbaqusComponentScaleWithoutClampOrRowDenominator) { + 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 fesaValues = reference; + fesaValues[0U].values[0U] = 1.5e-6; + fesaValues[1U].values[0U] += 2.002e-6; + fesaValues[0U].values[1U] = 0.999e-9; + fesaValues[1U].values[1U] = 1.001e-9; + writeHdf5(fixture.results(), fixture.input(), fesaValues); + + auto result = fesa::test::Mitc4ReferenceComparison::compare( + fixture.referenceCase()); + ASSERT_TRUE(result.hasValue()); + const auto& report = result.value(); + EXPECT_FALSE(report.passed); + const auto* u1Zero = findRow(report, 1, "U1"); + const auto* u1Scaled = findRow(report, 2, "U1"); + const auto* u2Near = findRow(report, 1, "U2"); + const auto* u2Over = findRow(report, 2, "U2"); + ASSERT_NE(u1Zero, nullptr); + ASSERT_NE(u1Scaled, nullptr); + ASSERT_NE(u2Near, nullptr); + ASSERT_NE(u2Over, nullptr); + EXPECT_DOUBLE_EQ(u1Zero->tolerance, 2.001e-6); + EXPECT_TRUE(u1Zero->withinTolerance); + EXPECT_FALSE(u1Scaled->withinTolerance); + EXPECT_DOUBLE_EQ(u2Near->tolerance, 1.0e-9); + EXPECT_TRUE(u2Near->withinTolerance); + EXPECT_FALSE(u2Over->withinTolerance); + const auto* u1Metric = findMetric(report, "U1"); + ASSERT_NE(u1Metric, nullptr); + EXPECT_DOUBLE_EQ(u1Metric->referenceScale, 2.0); +} + +// MITC4-REF-004 +TEST(Mitc4ReferenceComparison, RotationExceedanceWarnsWithoutBlockingTranslationVerdict) { + ContractFixture fixture{"warning"}; + auto fesaValues = defaultRows(); + fesaValues[0U].values[3U] = 1.0; + writeHdf5(fixture.results(), fixture.input(), fesaValues); + + 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-exceedance"); + const auto* row = findRow(report, 1, "UR1"); + ASSERT_NE(row, nullptr); + EXPECT_FALSE(row->blocking); + EXPECT_FALSE(row->withinTolerance); + EXPECT_EQ(report.warnings[0U].row, + static_cast(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); +} + +// MITC4-REF-005 +TEST(Mitc4ReferenceComparison, ReportsMetricsVectorsWorstRowAndJsonDeterministically) { + ContractFixture fixture{"report"}; + auto fesaValues = defaultRows(); + fesaValues[0U].values[0U] = 0.5e-9; + fesaValues[0U].values[1U] = -0.25e-9; + fesaValues[0U].values[3U] = 0.75e-9; + writeHdf5(fixture.results(), fixture.input(), fesaValues); + + 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(), 6U); + ASSERT_EQ(report.vectorMetrics.size(), 2U); + EXPECT_NEAR( + report.vectorMetrics[0U].displacementNormError, + std::sqrt(0.3125) * 1.0e-9, + 1.0e-21); + EXPECT_DOUBLE_EQ( + report.vectorMetrics[0U].rotationNormError, 0.75e-9); + ASSERT_LT(report.worstRow, report.rows.size()); + EXPECT_EQ(report.rows[report.worstRow].component, "UR1"); + for (const auto& metric : report.metrics) { + EXPECT_TRUE(std::isfinite(metric.referenceScale)); + EXPECT_TRUE(std::isfinite(metric.maximumAbsoluteError)); + EXPECT_TRUE(std::isfinite(metric.maximumNormalizedError)); + EXPECT_TRUE(std::isfinite(metric.rmsError)); + EXPECT_TRUE(std::isfinite(metric.vectorNormError)); + EXPECT_LT(metric.worstRow, report.rows.size()); + } + + const auto jsonA = fixture.root() / "comparison-a.json"; + const auto jsonB = fixture.root() / "comparison-b.json"; + ASSERT_TRUE( + fesa::test::Mitc4ReferenceComparison::writeDeterministicJson( + report, jsonA) + .isOk()); + ASSERT_TRUE( + fesa::test::Mitc4ReferenceComparison::writeDeterministicJson( + report, jsonB) + .isOk()); + const std::string first = readBytes(jsonA); + EXPECT_EQ(first, readBytes(jsonB)); + for (const char* key : { + "\"rows\"", "\"metrics\"", "\"vector_metrics\"", + "\"warnings\"", "\"worst_row\"", "\"passed\""}) { + EXPECT_NE(first.find(key), std::string::npos) << key; + } +} + +// MITC4-REF-006 +TEST(Mitc4ReferenceComparison, RequiresOnlyDeclaredInputCsvAndHdf5) { + ContractFixture fixture{"minimal-artifacts"}; + std::vector 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{"case.inp", "displacements.csv", "results.h5"})); + auto result = fesa::test::Mitc4ReferenceComparison::compare( + fixture.referenceCase()); + ASSERT_TRUE(result.hasValue()); + EXPECT_TRUE(result.value().passed); +} + +} // namespace