feat(linear-static-mitc4-shell): step 13 - shell-reference-comparison

This commit is contained in:
KOKO\Mimi
2026-08-12 22:18:45 +09:00
parent 492c89f283
commit 613bdb9db5
5 changed files with 2006 additions and 0 deletions
+3
View File
@@ -59,6 +59,9 @@ add_executable(
reference/reference_comparison.cpp reference/reference_comparison.cpp
reference/reference_comparison_test.cpp reference/reference_comparison_test.cpp
reference/b33_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( target_link_libraries(
@@ -0,0 +1,194 @@
#include "mitc4_reference_comparison.hpp"
#include "fesa/app/fesa_application.hpp"
#include <gtest/gtest.h>
#include <algorithm>
#include <cstddef>
#include <filesystem>
#include <fstream>
#include <iterator>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
#ifndef FESA_TEST_SOURCE_DIR
#error FESA_TEST_SOURCE_DIR must identify the repository root.
#endif
#ifndef FESA_TEST_BINARY_DIR
#error FESA_TEST_BINARY_DIR must identify the CMake binary root.
#endif
namespace {
constexpr const char* 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<char>{stream},
std::istreambuf_iterator<char>{}};
}
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<std::string> 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<std::string>{"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::size_t>(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::size_t>(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
@@ -0,0 +1,936 @@
#include "mitc4_reference_comparison.hpp"
#include <hdf5.h>
#include <algorithm>
#include <array>
#include <cerrno>
#include <cctype>
#include <charconv>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <locale>
#include <map>
#include <sstream>
#include <stdexcept>
#include <string>
#include <system_error>
#include <utility>
#include <vector>
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<const char*, 6> kComponents{
"U1", "U2", "U3", "UR1", "UR2", "UR3"};
const std::vector<std::string> 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<unsigned char>(character));
});
const auto end = std::find_if_not(
value.rbegin(), value.rend(), [&](const char character) {
return isSpace(static_cast<unsigned char>(character));
}).base();
return begin < end ? std::string{begin, end} : std::string{};
}
std::vector<std::string> splitCsvLine(const std::string& line) {
std::vector<std::string> fields;
std::size_t start = 0U;
while (true) {
const std::size_t comma = line.find(',', start);
fields.push_back(trim(line.substr(start, comma - start)));
if (comma == std::string::npos) {
break;
}
start = comma + 1U;
}
return fields;
}
std::int64_t parsePositiveLabel(const std::string& field) {
std::int64_t value = 0;
const char* const begin = field.data();
const char* const end = begin + field.size();
const auto parsed = std::from_chars(begin, end, value);
if (parsed.ec != std::errc{} || parsed.ptr != end || value <= 0) {
fail("schema-mismatch", "A 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<double, 6> values;
};
std::vector<WideRow> 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<WideRow> rows;
std::map<IdentityKey, std::size_t> 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<hsize_t> 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<hsize_t> result(static_cast<std::size_t>(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<const char*>& expected) {
Hdf5Handle type{
requireId(H5Dget_type(dataset), "Unable to inspect an HDF5 compound type."),
H5Tclose};
if (H5Tget_class(type.get()) != H5T_COMPOUND ||
H5Tget_nmembers(type.get()) != static_cast<int>(expected.size())) {
fail("schema-mismatch", "An HDF5 compound 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<unsigned>(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<WideRow> 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<hsize_t>((std::numeric_limits<std::size_t>::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<std::size_t>(shape[0U]);
std::vector<NodeReadRow> 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<WideRow> rows;
rows.reserve(count);
std::map<IdentityKey, std::size_t> 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<hsize_t>((std::numeric_limits<std::size_t>::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<std::size_t>(shape[0U]);
std::vector<ElementIdentityReadRow> 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<double> readDisplacement(
const hid_t file, const std::size_t nodeCount) {
auto dataset = openDataset(file, kDisplacementPath);
if (dimensions(dataset.get()) !=
std::vector<hsize_t>{static_cast<hsize_t>(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<double> 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<WideRow> 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<std::ptrdiff_t>(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<double>::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<unsigned>(character)
<< std::dec << std::setfill(' ');
} else {
stream << static_cast<char>(character);
}
break;
}
}
return stream.str();
}
} // namespace
Result<Mitc4ComparisonReport> Mitc4ReferenceComparison::compare(
const Mitc4ReferenceCase& referenceCase) {
try {
requireArtifacts(referenceCase);
const auto referenceRows =
readReferenceCsv(referenceCase.displacementCsvPath);
auto hdf5 = readHdf5(referenceCase);
std::map<IdentityKey, const WideRow*> 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<double, 6> 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<double, 6> errorNorms{};
std::array<double, 6> maximumErrors{};
std::array<double, 6> maximumNormalized{};
std::array<std::size_t, 6> 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<double, 6> 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<double>(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<Mitc4ComparisonReport>::success(std::move(report));
} catch (const ComparisonFailure& exception) {
return Result<Mitc4ComparisonReport>::failure(failureStatus(
referenceCase.caseId, exception.code(), exception.what()));
} catch (const std::exception& exception) {
return Result<Mitc4ComparisonReport>::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<double>::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
@@ -0,0 +1,81 @@
#pragma once
#include "fesa/core/status.hpp"
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <string>
#include <vector>
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<Mitc4RowDecision> rows;
std::vector<Mitc4ComponentMetrics> metrics;
std::vector<Mitc4VectorMetrics> vectorMetrics;
std::vector<Mitc4Warning> warnings;
std::size_t worstRow;
bool passed;
};
class Mitc4ReferenceComparison {
public:
static Result<Mitc4ComparisonReport> compare(
const Mitc4ReferenceCase& referenceCase);
static Status writeDeterministicJson(
const Mitc4ComparisonReport& report,
const std::filesystem::path& outputJson);
};
} // namespace fesa::test
@@ -0,0 +1,792 @@
#include "mitc4_reference_comparison.hpp"
#include <hdf5.h>
#include <gtest/gtest.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 instanceName{kInstanceName};
std::int64_t sourceNodeLabel{};
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 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<ComparisonValueRow>& values) {
std::vector<std::string> labels;
labels.reserve(values.size());
for (const auto& row : values) {
labels.push_back(std::to_string(row.sourceNodeLabel));
}
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].instanceName.c_str(),
labels[index].c_str(),
{static_cast<double>(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<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", 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<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.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<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(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<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.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<fesa::test::Mitc4ComparisonReport>& 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<char>{stream},
std::istreambuf_iterator<char>{}};
}
// 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<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].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<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);
}
// 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<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