baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
681794289411ebab68d03f554c70129f35a6c845
FESADev/tests/reference/nonlinear_static_reference_comparison.cpp
T

1723 lines
72 KiB
C++
Raw Blame History

#include "nonlinear_static_reference_comparison.h"
#include <hdf5.h>
#include <algorithm>
#include <array>
#include <atomic>
#include <charconv>
#include <chrono>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <map>
#include <optional>
#include <sstream>
#include <stdexcept>
#include <string>
#include <system_error>
#include <utility>
#include <vector>
#include "fesa/app/fesa_application.h"
#include "fesa/elements/finite_rotation.h"
#include "fesa/math/vector3.h"
#include "mitc4_reference_comparison.h"
#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 fesa::test {
namespace {
constexpr std::array<const char*, 6> kDisplacementComponents{
"U1", "U2", "U3", "UR1", "UR2", "UR3"};
constexpr std::array<const char*, 6> kReactionComponents{
"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"};
constexpr double kInternalGlobalBalanceRelativeTolerance = 1.0e-8;
constexpr const char* kNonlinearFrame =
"Increment 10: Step Time = 1.000";
constexpr std::array<std::int64_t, 24> kCoarseConstrainedNodeLabels{
1, 3, 4, 5, 6, 7, 8, 9, 14, 15, 16, 17,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33};
[[noreturn]] void Fail(const std::string& code, const std::string& message);
struct CaseContract {
std::size_t node_count;
bool nonlinear;
std::string family_prefix;
};
const CaseContract& ContractFor(const NonlinearReferenceCase& reference_case) {
static const CaseContract kCoarse{49U, true, "coarse."};
static const CaseContract kFine{729U, true, "fine."};
static const CaseContract kPressure{729U, true, "pressure."};
static const CaseContract kLinear{49U, false, "linear."};
if (reference_case.case_id == "shellnl-s4-final") {
return kCoarse;
}
if (reference_case.case_id == "shellnl-meshfine-s4-final") {
return kFine;
}
if (reference_case.case_id == "shellnl-meshfine-pressure-s4-follower") {
return kPressure;
}
if (reference_case.case_id == "shell-s4-linear-limit") {
return kLinear;
}
Fail("needs-reference-artifacts", "The reference case identifier is unknown.");
}
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 NonlinearReferenceCase& reference_case,
const std::string& code, const std::string& message) {
return Status::Failure(FailureCategory::kModel,
{{Severity::kError, code, {}, "",
reference_case.case_id, message}});
}
std::filesystem::path Normalized(const std::filesystem::path& path) {
std::error_code error;
const std::filesystem::path absolute = std::filesystem::absolute(path, error);
return (error ? path : absolute).lexically_normal();
}
bool IsDescendantOf(const std::filesystem::path& path,
const std::filesystem::path& parent) {
const auto normalized_path = Normalized(path);
const auto normalized_parent = Normalized(parent);
auto path_part = normalized_path.begin();
auto parent_part = normalized_parent.begin();
for (; parent_part != normalized_parent.end(); ++parent_part, ++path_part) {
if (path_part == normalized_path.end() || *path_part != *parent_part) {
return false;
}
}
return true;
}
std::string GenericPath(const std::filesystem::path& path) {
return Normalized(path).generic_u8string();
}
std::string Trim(const std::string& value) {
const auto is_space = [](const unsigned char character) {
return character == ' ' || character == '\t' || character == '\r' ||
character == '\n' || character == '\f' || character == '\v';
};
const auto begin = std::find_if_not(
value.begin(), value.end(),
[&is_space](const char character) {
return is_space(static_cast<unsigned char>(character));
});
const auto end = std::find_if_not(
value.rbegin(), value.rend(),
[&is_space](const char character) {
return is_space(static_cast<unsigned char>(character));
})
.base();
return begin < end ? std::string{begin, end} : std::string{};
}
std::string UppercaseAscii(std::string value) {
std::transform(value.begin(), value.end(), value.begin(),
[](const char character) {
return character >= 'a' && character <= 'z'
? static_cast<char>(character - 'a' + 'A')
: character;
});
return value;
}
std::string ReadBytes(const std::filesystem::path& path) {
std::ifstream stream{path, std::ios::binary};
if (!stream) {
Fail("needs-reference-artifacts",
"A declared reference artifact is missing or unreadable.");
}
return {std::istreambuf_iterator<char>{stream},
std::istreambuf_iterator<char>{}};
}
std::uint32_t RotateRight(const std::uint32_t value,
const std::uint32_t amount) {
return (value >> amount) | (value << (32U - amount));
}
std::string Sha256(const std::string& input) {
constexpr std::array<std::uint32_t, 64> kRoundConstants{
0x428a2f98U, 0x71374491U, 0xb5c0fbcfU, 0xe9b5dba5U, 0x3956c25bU,
0x59f111f1U, 0x923f82a4U, 0xab1c5ed5U, 0xd807aa98U, 0x12835b01U,
0x243185beU, 0x550c7dc3U, 0x72be5d74U, 0x80deb1feU, 0x9bdc06a7U,
0xc19bf174U, 0xe49b69c1U, 0xefbe4786U, 0x0fc19dc6U, 0x240ca1ccU,
0x2de92c6fU, 0x4a7484aaU, 0x5cb0a9dcU, 0x76f988daU, 0x983e5152U,
0xa831c66dU, 0xb00327c8U, 0xbf597fc7U, 0xc6e00bf3U, 0xd5a79147U,
0x06ca6351U, 0x14292967U, 0x27b70a85U, 0x2e1b2138U, 0x4d2c6dfcU,
0x53380d13U, 0x650a7354U, 0x766a0abbU, 0x81c2c92eU, 0x92722c85U,
0xa2bfe8a1U, 0xa81a664bU, 0xc24b8b70U, 0xc76c51a3U, 0xd192e819U,
0xd6990624U, 0xf40e3585U, 0x106aa070U, 0x19a4c116U, 0x1e376c08U,
0x2748774cU, 0x34b0bcb5U, 0x391c0cb3U, 0x4ed8aa4aU, 0x5b9cca4fU,
0x682e6ff3U, 0x748f82eeU, 0x78a5636fU, 0x84c87814U, 0x8cc70208U,
0x90befffaU, 0xa4506cebU, 0xbef9a3f7U, 0xc67178f2U};
std::vector<std::uint8_t> bytes;
bytes.reserve(input.size() + 72U);
for (const unsigned char byte : input) {
bytes.push_back(byte);
}
const std::uint64_t bit_length =
static_cast<std::uint64_t>(bytes.size()) * 8U;
bytes.push_back(0x80U);
while ((bytes.size() % 64U) != 56U) {
bytes.push_back(0U);
}
for (int shift = 56; shift >= 0; shift -= 8) {
bytes.push_back(static_cast<std::uint8_t>(bit_length >> shift));
}
std::array<std::uint32_t, 8> state{0x6a09e667U, 0xbb67ae85U, 0x3c6ef372U,
0xa54ff53aU, 0x510e527fU, 0x9b05688cU,
0x1f83d9abU, 0x5be0cd19U};
for (std::size_t offset = 0U; offset < bytes.size(); offset += 64U) {
std::array<std::uint32_t, 64> words{};
for (std::size_t index = 0U; index < 16U; ++index) {
const std::size_t byte = offset + index * 4U;
words[index] = (static_cast<std::uint32_t>(bytes[byte]) << 24U) |
(static_cast<std::uint32_t>(bytes[byte + 1U]) << 16U) |
(static_cast<std::uint32_t>(bytes[byte + 2U]) << 8U) |
static_cast<std::uint32_t>(bytes[byte + 3U]);
}
for (std::size_t index = 16U; index < words.size(); ++index) {
const std::uint32_t sigma0 = RotateRight(words[index - 15U], 7U) ^
RotateRight(words[index - 15U], 18U) ^
(words[index - 15U] >> 3U);
const std::uint32_t sigma1 = RotateRight(words[index - 2U], 17U) ^
RotateRight(words[index - 2U], 19U) ^
(words[index - 2U] >> 10U);
words[index] = words[index - 16U] + sigma0 + words[index - 7U] +
sigma1;
}
std::uint32_t a = state[0U];
std::uint32_t b = state[1U];
std::uint32_t c = state[2U];
std::uint32_t d = state[3U];
std::uint32_t e = state[4U];
std::uint32_t f = state[5U];
std::uint32_t g = state[6U];
std::uint32_t h = state[7U];
for (std::size_t index = 0U; index < words.size(); ++index) {
const std::uint32_t sigma1 = RotateRight(e, 6U) ^ RotateRight(e, 11U) ^
RotateRight(e, 25U);
const std::uint32_t choose = (e & f) ^ (~e & g);
const std::uint32_t temporary1 =
h + sigma1 + choose + kRoundConstants[index] + words[index];
const std::uint32_t sigma0 = RotateRight(a, 2U) ^ RotateRight(a, 13U) ^
RotateRight(a, 22U);
const std::uint32_t majority = (a & b) ^ (a & c) ^ (b & c);
const std::uint32_t temporary2 = sigma0 + majority;
h = g;
g = f;
f = e;
e = d + temporary1;
d = c;
c = b;
b = a;
a = temporary1 + temporary2;
}
state[0U] += a;
state[1U] += b;
state[2U] += c;
state[3U] += d;
state[4U] += e;
state[5U] += f;
state[6U] += g;
state[7U] += h;
}
std::ostringstream result;
result << std::uppercase << std::hex << std::setfill('0');
for (const std::uint32_t word : state) {
result << std::setw(8) << word;
}
return result.str();
}
std::string Fnv1a64(const std::string& input) {
std::uint64_t hash = 14695981039346656037ULL;
for (const unsigned char byte : input) {
hash ^= static_cast<std::uint64_t>(byte);
hash *= 1099511628211ULL;
}
std::ostringstream result;
result << "fnv1a64:" << std::hex << std::setfill('0') << std::setw(16)
<< hash;
return result.str();
}
struct DeclaredArtifact {
std::filesystem::path path;
const char* sha256;
std::uintmax_t expected_byte_size{0U};
NonlinearReferenceArtifactDisposition disposition{
NonlinearReferenceArtifactDisposition::kRequired};
};
std::vector<DeclaredArtifact> DeclaredArtifacts(
const NonlinearReferenceCase& reference_case) {
const std::filesystem::path root{FESA_TEST_SOURCE_DIR};
if (reference_case.case_id == "shellnl-meshfine-pressure-s4-follower") {
const auto directory = root / "reference" / "shellNL-meshfine-pressure";
const std::vector<DeclaredArtifact> artifacts{
{directory / "shellNL-meshfine-pressure.inp",
"53718CCC9B41C23FCC1B509A51D66F1347E87ACA43A74F1D6CEC4D2085F4B90E",
56665U, reference_case.input_artifact_disposition},
{directory / "shellNL-meshfine-pressure displacements.csv",
"2FE53F98EAC9F04779100721C3EDF8049EB1DBFE733EBCBD819AEEED03B81198",
111182U, reference_case.displacement_artifact_disposition},
{directory / "shellNL-meshfine-pressure reactions.csv", nullptr, 0U,
reference_case.reaction_artifact_disposition},
{directory / "shellNL-meshfine-pressure stresses.csv", nullptr, 0U,
reference_case.stress_artifact_disposition}};
if (Normalized(reference_case.original_input_path) !=
Normalized(artifacts[0U].path) ||
Normalized(reference_case.displacement_csv_path) !=
Normalized(artifacts[1U].path) ||
Normalized(reference_case.reaction_csv_path) !=
Normalized(artifacts[2U].path) ||
Normalized(reference_case.stress_inventory_path) !=
Normalized(artifacts[3U].path) ||
reference_case.input_artifact_disposition !=
NonlinearReferenceArtifactDisposition::kRequired ||
reference_case.displacement_artifact_disposition !=
NonlinearReferenceArtifactDisposition::kRequired ||
reference_case.reaction_artifact_disposition !=
NonlinearReferenceArtifactDisposition::kReviewOnly ||
reference_case.stress_artifact_disposition !=
NonlinearReferenceArtifactDisposition::kExcluded) {
Fail("needs-reference-artifacts",
"The pressure reference case does not use its declared artifact "
"dispositions.");
}
return artifacts;
}
if (reference_case.case_id == "shellnl-meshfine-s4-final") {
const auto directory = root / "reference" / "shellNL-meshfine";
const std::vector<DeclaredArtifact> artifacts{
{directory / "shellNL-meshfine.inp",
"40DF1CE77311C8D799564A21A60F3BDEBB8148D52CAE22B7292FDFAD32018220"},
{directory / "shellNL-meshfine displacements.csv",
"740819FF8AED0E279961581F7CDBFA12784CFBA3A488C35A652E80103F64433F"},
{directory / "shellNL-meshfine reactions.csv",
"AEB01C98C596792B0E5F2B86A81945CAB0BEB0C6969CFC6738C6F81DF1F94363"},
{directory / "shellNL-meshfine stresses.csv",
"6D6213AB5DB9D6B1F5FBD641367D08BACC805967B6C5A95BCE0E1A5395377782"}};
if (Normalized(reference_case.original_input_path) !=
Normalized(artifacts[0U].path) ||
Normalized(reference_case.displacement_csv_path) !=
Normalized(artifacts[1U].path) ||
Normalized(reference_case.reaction_csv_path) !=
Normalized(artifacts[2U].path) ||
Normalized(reference_case.stress_inventory_path) !=
Normalized(artifacts[3U].path)) {
Fail("needs-reference-artifacts",
"The fine nonlinear reference case does not use declared paths.");
}
return artifacts;
}
if (reference_case.case_id == "shellnl-s4-final") {
const auto directory = root / "reference" / "shellNL";
const std::vector<DeclaredArtifact> artifacts{
{directory / "shellNL.inp",
"4BF3ECB0F7F1E8462C573E993DBC97AC5D920C91CF6F23AD4C6891D77FFBADA2"},
{directory / "shellNL displacements.csv",
"4E0F99CA9372B1ADF7FAF102F43BC0B887DC743FDEDA378E4342D2CF8E0C0786"},
{directory / "shellNL reactions.csv",
"001CB1541A9E94424E5E650509D46BC8724970DF4DE47184D9DB09C7938AB069"},
{directory / "shellNL stresses.csv",
"5500E5114B1C44887528E30CEC054A1EB0F68BE3F4B8CE81FA59FC74A881483F"}};
if (Normalized(reference_case.original_input_path) !=
Normalized(artifacts[0U].path) ||
Normalized(reference_case.displacement_csv_path) !=
Normalized(artifacts[1U].path) ||
Normalized(reference_case.reaction_csv_path) !=
Normalized(artifacts[2U].path) ||
Normalized(reference_case.stress_inventory_path) !=
Normalized(artifacts[3U].path)) {
Fail("needs-reference-artifacts",
"The nonlinear reference case does not use its declared paths.");
}
return artifacts;
}
if (reference_case.case_id == "shell-s4-linear-limit") {
const auto directory = root / "reference" / "shell";
const std::vector<DeclaredArtifact> artifacts{
{directory / "shell.inp",
"4005851E1AB22FD3A16AC17A8D5DA3E051233F69F37419079F3553AD134ECFCF"},
{directory / "shell displacements.csv",
"C81D94E0B4A849F87AA0F79C83A79B94D5661AC79E44ED826919AB432C87746B"},
{directory / "shell reactions.csv",
"0E34EEC27E700C3C0C744C1DAFD8A3B95EFEF34B8616FE96665876586FE82D41"},
{directory / "shell stresses.csv",
"2B75A5B0DEB045A283C9BB7E9B1153B79F7537FF487214BF4A23316B59B660A7"}};
if (Normalized(reference_case.original_input_path) !=
Normalized(artifacts[0U].path) ||
Normalized(reference_case.displacement_csv_path) !=
Normalized(artifacts[1U].path) ||
Normalized(reference_case.reaction_csv_path) !=
Normalized(artifacts[2U].path) ||
Normalized(reference_case.stress_inventory_path) !=
Normalized(artifacts[3U].path)) {
Fail("needs-reference-artifacts",
"The linear control case does not use its declared paths.");
}
return artifacts;
}
Fail("needs-reference-artifacts", "The reference case identifier is unknown.");
}
ArtifactSnapshot Snapshot(const DeclaredArtifact& artifact) {
std::error_code error;
if (!std::filesystem::is_regular_file(artifact.path, error) || error) {
Fail("needs-reference-artifacts",
"A declared reference artifact is missing or unreadable.");
}
const std::string bytes = ReadBytes(artifact.path);
const std::string hash = Sha256(bytes);
if (artifact.sha256 == nullptr || hash != artifact.sha256) {
Fail("reference-artifact-hash-mismatch",
"A declared reference artifact does not match its approved hash.");
}
const std::uintmax_t size = std::filesystem::file_size(artifact.path, error);
if (error) {
Fail("needs-reference-artifacts",
"A declared reference artifact size cannot be read.");
}
if (artifact.expected_byte_size != 0U &&
size != artifact.expected_byte_size) {
Fail("reference-artifact-size-mismatch",
"A declared reference artifact does not match its approved size.");
}
const auto write_time = std::filesystem::last_write_time(artifact.path, error);
if (error) {
Fail("needs-reference-artifacts",
"A declared reference artifact timestamp cannot be read.");
}
return {Normalized(artifact.path), hash, size, write_time};
}
std::vector<std::string> SplitCsv(const std::string& line) {
std::vector<std::string> fields;
std::string field;
bool quoted = false;
for (std::size_t index = 0U; index < line.size(); ++index) {
const char character = line[index];
if (character == '"') {
if (quoted && index + 1U < line.size() && line[index + 1U] == '"') {
field.push_back(character);
++index;
} else {
quoted = !quoted;
}
} else if (character == ',' && !quoted) {
fields.push_back(Trim(field));
field.clear();
} else {
field.push_back(character);
}
}
if (quoted) {
Fail("schema-mismatch", "A reference CSV row has an unclosed quote.");
}
fields.push_back(Trim(field));
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 reference CSV node label is invalid.");
}
return value;
}
double ParseFiniteDouble(const std::string& field) {
char* end = nullptr;
const double value = std::strtod(field.c_str(), &end);
if (field.empty() || end == field.c_str() || end == nullptr ||
*end != '\0' || !std::isfinite(value)) {
Fail("schema-mismatch", "A reference CSV value is nonfinite or invalid.");
}
return value;
}
struct CsvWideRow {
std::string instance_name;
std::int64_t source_node_label;
std::array<double, 6> values;
};
std::vector<CsvWideRow> ReadCsvWideRows(
const std::filesystem::path& path, const bool requires_frame,
const std::array<const char*, 6>& components,
const std::size_t expected_node_count) {
std::ifstream stream{path, std::ios::binary};
if (!stream) {
Fail("needs-reference-artifacts",
"A declared comparison CSV is missing or unreadable.");
}
std::string line;
if (!std::getline(stream, line)) {
Fail("schema-mismatch", "A declared comparison CSV is empty.");
}
if (!line.empty() && line.back() == '\r') {
line.pop_back();
}
std::vector<std::string> expected_header;
if (requires_frame) {
expected_header.push_back("Frame");
}
expected_header.push_back("Part Instance Name");
expected_header.push_back("Node Label");
for (const char* component : components) {
expected_header.push_back(std::string{components == kDisplacementComponents
? (std::string{"U-"} + component)
: (std::string{"RF-"} + component)});
}
if (components == kDisplacementComponents) {
expected_header[requires_frame ? 3U : 2U] = "U-U1";
expected_header[requires_frame ? 4U : 3U] = "U-U2";
expected_header[requires_frame ? 5U : 4U] = "U-U3";
expected_header[requires_frame ? 6U : 5U] = "UR-UR1";
expected_header[requires_frame ? 7U : 6U] = "UR-UR2";
expected_header[requires_frame ? 8U : 7U] = "UR-UR3";
} else {
expected_header[requires_frame ? 3U : 2U] = "RF-RF1";
expected_header[requires_frame ? 4U : 3U] = "RF-RF2";
expected_header[requires_frame ? 5U : 4U] = "RF-RF3";
expected_header[requires_frame ? 6U : 5U] = "RM-RM1";
expected_header[requires_frame ? 7U : 6U] = "RM-RM2";
expected_header[requires_frame ? 8U : 7U] = "RM-RM3";
}
if (SplitCsv(line) != expected_header) {
Fail("schema-mismatch", "A comparison CSV header is out of contract.");
}
std::vector<CsvWideRow> rows;
std::map<std::pair<std::string, std::int64_t>, std::size_t> identities;
while (std::getline(stream, line)) {
if (!line.empty() && line.back() == '\r') {
line.pop_back();
}
if (line.empty()) {
Fail("schema-mismatch", "A comparison CSV has a blank row.");
}
const auto fields = SplitCsv(line);
if (fields.size() != expected_header.size()) {
Fail("schema-mismatch", "A comparison CSV row has the wrong width.");
}
const std::size_t shift = requires_frame ? 1U : 0U;
if (requires_frame && fields[0U] != kNonlinearFrame) {
Fail("schema-mismatch", "A nonlinear CSV row has the wrong Frame value.");
}
if (fields[shift].empty()) {
Fail("schema-mismatch", "A comparison CSV instance is empty.");
}
CsvWideRow row{fields[shift], ParsePositiveLabel(fields[shift + 1U]), {}};
for (std::size_t component = 0U; component < row.values.size();
++component) {
row.values[component] = ParseFiniteDouble(fields[shift + 2U + component]);
}
const auto identity =
std::make_pair(UppercaseAscii(row.instance_name), row.source_node_label);
if (!identities.emplace(identity, rows.size()).second) {
Fail("schema-mismatch", "A comparison CSV wide-row identity is duplicated.");
}
rows.push_back(std::move(row));
}
if (rows.size() != expected_node_count) {
Fail("schema-mismatch", "A comparison CSV node inventory is out of contract.");
}
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_, &client_data_) >= 0 &&
H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr) >= 0) {
active_ = true;
}
}
~Hdf5ErrorSilencer() {
if (active_) {
(void)H5Eset_auto2(H5E_DEFAULT, callback_, client_data_);
}
}
private:
H5E_auto2_t callback_{nullptr};
void* client_data_{nullptr};
bool active_{false};
};
class Hdf5VlenReclaimer {
public:
Hdf5VlenReclaimer(const hid_t memory_type, const hid_t data_space,
void* const data)
: memory_type_{memory_type}, data_space_{data_space}, data_{data} {}
~Hdf5VlenReclaimer() {
if (active_) {
(void)H5Dvlen_reclaim(memory_type_, data_space_, H5P_DEFAULT, data_);
}
}
void Reclaim() {
active_ = false;
if (H5Dvlen_reclaim(memory_type_, data_space_, H5P_DEFAULT, data_) < 0) {
Fail("schema-mismatch", "HDF5 variable-length strings cannot be freed.");
}
}
private:
hid_t memory_type_;
hid_t data_space_;
void* data_;
bool active_{true};
};
hid_t RequireId(const hid_t value, const char* message) {
if (value < 0) {
Fail("schema-mismatch", message);
}
return value;
}
void RequireHdf(const herr_t value, const char* message) {
if (value < 0) {
Fail("schema-mismatch", message);
}
}
Hdf5Handle OpenDataset(const hid_t file, const std::string& path) {
return {RequireId(H5Dopen2(file, path.c_str(), H5P_DEFAULT),
"A required HDF5 dataset is missing."),
H5Dclose};
}
std::vector<hsize_t> Dimensions(const hid_t dataset) {
Hdf5Handle space{RequireId(H5Dget_space(dataset),
"HDF5 dimensions cannot be inspected."),
H5Sclose};
const int rank = H5Sget_simple_extent_ndims(space.Get());
if (rank < 0) {
Fail("schema-mismatch", "HDF5 rank cannot be inspected.");
}
std::vector<hsize_t> dimensions(static_cast<std::size_t>(rank));
if (rank > 0) {
RequireHdf(H5Sget_simple_extent_dims(space.Get(), dimensions.data(), nullptr),
"HDF5 extents cannot be inspected.");
}
return dimensions;
}
Hdf5Handle MakeUtf8StringType() {
Hdf5Handle type{RequireId(H5Tcopy(H5T_C_S1),
"An HDF5 string type cannot be created."),
H5Tclose};
RequireHdf(H5Tset_size(type.Get(), H5T_VARIABLE),
"An HDF5 variable string type cannot be configured.");
RequireHdf(H5Tset_cset(type.Get(), H5T_CSET_UTF8),
"An HDF5 UTF-8 string type cannot be configured.");
return type;
}
std::string ReadStringAttribute(const hid_t object, const char* name) {
Hdf5Handle attribute{RequireId(H5Aopen(object, name, H5P_DEFAULT),
"A required HDF5 attribute is missing."),
H5Aclose};
Hdf5Handle type{RequireId(H5Aget_type(attribute.Get()),
"An HDF5 attribute type cannot be read."),
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 is out of contract.");
}
char* raw = nullptr;
RequireHdf(H5Aread(attribute.Get(), type.Get(), &raw),
"An HDF5 string attribute cannot be read.");
if (raw == nullptr) {
Fail("schema-mismatch", "An HDF5 string attribute is null.");
}
const std::string result{raw};
RequireHdf(H5free_memory(raw), "HDF5 attribute memory cannot be freed.");
return result;
}
std::uint64_t ReadUint64Attribute(const hid_t object, const char* name) {
Hdf5Handle attribute{RequireId(H5Aopen(object, name, H5P_DEFAULT),
"A required HDF5 attribute is missing."),
H5Aclose};
std::uint64_t value = 0U;
RequireHdf(H5Aread(attribute.Get(), H5T_NATIVE_UINT64, &value),
"An HDF5 integer attribute cannot be read.");
return value;
}
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.");
}
}
struct NodeReadRow {
std::uint64_t internal_node_id;
char* instance_name;
char* source_label;
double coordinates[3];
};
struct Hdf5Node {
std::string instance_name;
std::int64_t source_node_label;
std::array<double, 3> coordinates;
};
std::vector<Hdf5Node> ReadNodes(const hid_t file,
const std::size_t expected_node_count) {
auto dataset = OpenDataset(file, "/model/nodes");
const auto shape = Dimensions(dataset.Get());
if (shape.size() != 1U || shape[0U] != expected_node_count) {
Fail("schema-mismatch", "The HDF5 node inventory is out of contract.");
}
auto string_type = MakeUtf8StringType();
const hsize_t coordinate_dimensions[] = {3U};
Hdf5Handle coordinates{RequireId(H5Tarray_create2(
H5T_NATIVE_DOUBLE, 1, coordinate_dimensions),
"An HDF5 coordinate type cannot be created."),
H5Tclose};
Hdf5Handle memory_type{RequireId(H5Tcreate(H5T_COMPOUND, sizeof(NodeReadRow)),
"An HDF5 node type cannot be created."),
H5Tclose};
RequireHdf(H5Tinsert(memory_type.Get(), "internal_node_id",
HOFFSET(NodeReadRow, internal_node_id),
H5T_NATIVE_UINT64),
"An HDF5 node identifier cannot be read.");
RequireHdf(H5Tinsert(memory_type.Get(), "instance_name",
HOFFSET(NodeReadRow, instance_name), string_type.Get()),
"An HDF5 node instance cannot be read.");
RequireHdf(H5Tinsert(memory_type.Get(), "source_label",
HOFFSET(NodeReadRow, source_label), string_type.Get()),
"An HDF5 node label cannot be read.");
RequireHdf(H5Tinsert(memory_type.Get(), "coordinates",
HOFFSET(NodeReadRow, coordinates), coordinates.Get()),
"An HDF5 node coordinates cannot be read.");
std::vector<NodeReadRow> raw(static_cast<std::size_t>(shape[0U]));
RequireHdf(H5Dread(dataset.Get(), memory_type.Get(), H5S_ALL, H5S_ALL,
H5P_DEFAULT, raw.data()),
"HDF5 node identities cannot be read.");
Hdf5Handle space{RequireId(H5Dget_space(dataset.Get()),
"HDF5 node space cannot be read."),
H5Sclose};
Hdf5VlenReclaimer reclaimer{memory_type.Get(), space.Get(), raw.data()};
std::map<std::pair<std::string, std::int64_t>, std::size_t> identities;
std::vector<Hdf5Node> nodes;
nodes.reserve(raw.size());
for (std::size_t index = 0U; index < raw.size(); ++index) {
if (raw[index].internal_node_id != index || raw[index].instance_name == nullptr ||
raw[index].source_label == nullptr ||
!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 identity is invalid.");
}
Hdf5Node node{raw[index].instance_name,
ParsePositiveLabel(raw[index].source_label),
{raw[index].coordinates[0U], raw[index].coordinates[1U],
raw[index].coordinates[2U]}};
const auto identity =
std::make_pair(UppercaseAscii(node.instance_name), node.source_node_label);
if (!identities.emplace(identity, nodes.size()).second) {
Fail("schema-mismatch", "An HDF5 node identity is duplicated.");
}
nodes.push_back(std::move(node));
}
reclaimer.Reclaim();
return nodes;
}
struct AcceptedFrameReadRow {
std::uint64_t frame_index;
std::uint64_t increment_index;
double load_factor;
std::uint64_t iteration_count;
std::uint8_t is_factor_one;
};
std::uint64_t ReadFactorOneFrame(const hid_t file) {
auto dataset = OpenDataset(file, "/steps/Step-1/accepted_frames");
const auto shape = Dimensions(dataset.Get());
if (shape.size() != 1U || shape[0U] == 0U) {
Fail("schema-mismatch", "The accepted-frame table is invalid.");
}
Hdf5Handle type{RequireId(H5Tcreate(H5T_COMPOUND, sizeof(AcceptedFrameReadRow)),
"An accepted-frame type cannot be created."),
H5Tclose};
RequireHdf(H5Tinsert(type.Get(), "frame_index",
HOFFSET(AcceptedFrameReadRow, frame_index),
H5T_NATIVE_UINT64),
"An accepted-frame index cannot be read.");
RequireHdf(H5Tinsert(type.Get(), "increment_index",
HOFFSET(AcceptedFrameReadRow, increment_index),
H5T_NATIVE_UINT64),
"An accepted-frame increment cannot be read.");
RequireHdf(H5Tinsert(type.Get(), "load_factor",
HOFFSET(AcceptedFrameReadRow, load_factor),
H5T_NATIVE_DOUBLE),
"An accepted-frame load factor cannot be read.");
RequireHdf(H5Tinsert(type.Get(), "iteration_count",
HOFFSET(AcceptedFrameReadRow, iteration_count),
H5T_NATIVE_UINT64),
"An accepted-frame iteration count cannot be read.");
RequireHdf(H5Tinsert(type.Get(), "is_factor_one",
HOFFSET(AcceptedFrameReadRow, is_factor_one),
H5T_NATIVE_UINT8),
"An accepted-frame factor-one flag cannot be read.");
std::vector<AcceptedFrameReadRow> rows(static_cast<std::size_t>(shape[0U]));
RequireHdf(H5Dread(dataset.Get(), type.Get(), H5S_ALL, H5S_ALL, H5P_DEFAULT,
rows.data()),
"The accepted-frame table cannot be read.");
std::vector<std::uint64_t> factor_one;
for (const auto& row : rows) {
if (!std::isfinite(row.load_factor) || row.increment_index == 0U ||
row.iteration_count == 0U || row.is_factor_one > 1U) {
Fail("schema-mismatch", "An accepted-frame row is invalid.");
}
if (row.is_factor_one == 1U) {
if (row.load_factor != 1.0) {
Fail("schema-mismatch", "A factor-one frame must have load factor one.");
}
factor_one.push_back(row.frame_index);
}
}
if (factor_one.size() != 1U) {
Fail("schema-mismatch", "The HDF5 result needs one factor-one frame.");
}
return factor_one.front();
}
std::vector<double> ReadNodalDataset(const hid_t file, const std::string& path,
const std::string& expected_components,
const std::size_t node_count) {
auto dataset = OpenDataset(file, path);
if (Dimensions(dataset.Get()) != std::vector<hsize_t>{node_count, 6U}) {
Fail("schema-mismatch", "A nodal HDF5 dataset has the wrong shape.");
}
RequireStringAttribute(dataset.Get(), "component_names", expected_components);
std::vector<double> values(node_count * 6U);
RequireHdf(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
H5P_DEFAULT, values.data()),
"A nodal HDF5 dataset cannot be read.");
if (!std::all_of(values.begin(), values.end(),
[](const double value) { return std::isfinite(value); })) {
Fail("schema-mismatch", "A nodal HDF5 value is nonfinite.");
}
return values;
}
std::vector<std::uint8_t> ReadConstraintMask(const hid_t file,
const std::size_t node_count) {
auto dataset = OpenDataset(file, "/model/nodal_constraint_mask");
if (Dimensions(dataset.Get()) != std::vector<hsize_t>{node_count, 6U}) {
Fail("schema-mismatch", "The HDF5 constraint mask has the wrong shape.");
}
std::vector<std::uint8_t> values(node_count * 6U);
RequireHdf(H5Dread(dataset.Get(), H5T_NATIVE_UINT8, H5S_ALL, H5S_ALL,
H5P_DEFAULT, values.data()),
"The HDF5 constraint mask cannot be read.");
if (!std::all_of(values.begin(), values.end(),
[](const std::uint8_t value) { return value <= 1U; })) {
Fail("schema-mismatch", "The HDF5 constraint mask contains an invalid flag.");
}
return values;
}
struct Hdf5Projection {
std::filesystem::path copied_input_path;
std::vector<Hdf5Node> nodes;
std::vector<double> displacement;
std::vector<double> internal_force;
std::vector<double> external_force;
std::vector<double> residual;
std::vector<double> reaction;
std::vector<std::uint8_t> constraint_mask;
std::array<double, 6> global_equilibrium;
};
std::filesystem::path CopiedInputPathFor(
const std::filesystem::path& hdf5_results_path) {
const std::string stem = hdf5_results_path.stem().string();
constexpr const char* kSuffix = "_results";
const std::size_t suffix_size = std::char_traits<char>::length(kSuffix);
if (stem.size() <= suffix_size ||
stem.compare(stem.size() - suffix_size, suffix_size, kSuffix) != 0) {
Fail("schema-mismatch", "The generated HDF5 file does not have a sibling stem.");
}
return hdf5_results_path.parent_path() /
(stem.substr(0U, stem.size() - suffix_size) + ".inp");
}
std::array<double, 6> ReadGlobalBalance(const hid_t file,
const std::string& path) {
auto dataset = OpenDataset(file, path);
if (Dimensions(dataset.Get()) != std::vector<hsize_t>{6U}) {
Fail("schema-mismatch", "The global balance dataset has the wrong shape.");
}
RequireStringAttribute(dataset.Get(), "component_names",
"FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3");
RequireStringAttribute(dataset.Get(), "location", "global-balance-error");
std::array<double, 6> values{};
RequireHdf(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
H5P_DEFAULT, values.data()),
"The global balance dataset cannot be read.");
if (!std::all_of(values.begin(), values.end(),
[](const double value) { return std::isfinite(value); })) {
Fail("schema-mismatch", "The global balance dataset contains nonfinite data.");
}
return values;
}
Hdf5Projection ReadNonlinearHdf5(
const NonlinearReferenceCase& reference_case,
const std::filesystem::path& hdf5_results_path) {
Hdf5ErrorSilencer silence;
Hdf5Handle file{RequireId(H5Fopen(hdf5_results_path.string().c_str(),
H5F_ACC_RDONLY, H5P_DEFAULT),
"The generated HDF5 result 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") != 1U) {
Fail("schema-mismatch", "The nonlinear HDF5 schema version is invalid.");
}
RequireStringAttribute(metadata.Get(), "feature_id",
"nonlinear-static-newton-raphson");
RequireStringAttribute(metadata.Get(), "analysis_procedure",
"nonlinear-static-full-newton");
const CaseContract& contract = ContractFor(reference_case);
const std::filesystem::path copied_input = CopiedInputPathFor(hdf5_results_path);
const std::string expected_identity =
"path=" + GenericPath(copied_input) + ";content_identity=" +
Fnv1a64(ReadBytes(reference_case.original_input_path));
if (ReadStringAttribute(metadata.Get(), "source_input_identity") !=
expected_identity) {
Fail("schema-mismatch",
"The HDF5 source path or source content identity is out of contract.");
}
const auto nodes = ReadNodes(file.Get(), contract.node_count);
const std::uint64_t frame_index = ReadFactorOneFrame(file.Get());
const std::string root = "/steps/Step-1/frames/" +
std::to_string(frame_index) + "/nodal/";
return {copied_input,
nodes,
ReadNodalDataset(file.Get(), root + "displacement",
"UX,UY,UZ,URX,URY,URZ", contract.node_count),
ReadNodalDataset(file.Get(), root + "generalized_internal_force",
"FI1,FI2,FI3,MI_PHI1,MI_PHI2,MI_PHI3",
contract.node_count),
ReadNodalDataset(file.Get(), root + "generalized_external_force",
"FE1,FE2,FE3,ME_PHI1,ME_PHI2,ME_PHI3",
contract.node_count),
ReadNodalDataset(file.Get(), root + "generalized_residual",
"RU1,RU2,RU3,RPHI1,RPHI2,RPHI3",
contract.node_count),
ReadNodalDataset(file.Get(), root + "reaction",
"RF1,RF2,RF3,RM1,RM2,RM3", contract.node_count),
ReadConstraintMask(file.Get(), contract.node_count),
ReadGlobalBalance(file.Get(),
"/steps/Step-1/frames/" +
std::to_string(frame_index) +
"/global/equilibrium")};
}
void RequireExactNodeLabels(const std::vector<CsvWideRow>& rows,
const std::size_t expected_node_count) {
std::vector<std::int64_t> labels;
labels.reserve(rows.size());
for (const auto& row : rows) {
labels.push_back(row.source_node_label);
}
std::sort(labels.begin(), labels.end());
if (labels.size() != expected_node_count) {
Fail("schema-mismatch", "The CSV node inventory is out of contract.");
}
for (std::size_t index = 0U; index < labels.size(); ++index) {
if (labels[index] != static_cast<std::int64_t>(index + 1U)) {
Fail("schema-mismatch", "The CSV node labels are out of contract.");
}
}
}
void RequireExactNodeLabels(const std::vector<Hdf5Node>& nodes,
const std::size_t expected_node_count) {
std::vector<std::int64_t> labels;
labels.reserve(nodes.size());
for (const auto& node : nodes) {
labels.push_back(node.source_node_label);
}
std::sort(labels.begin(), labels.end());
if (labels.size() != expected_node_count) {
Fail("schema-mismatch", "The HDF5 node inventory is out of contract.");
}
for (std::size_t index = 0U; index < labels.size(); ++index) {
if (labels[index] != static_cast<std::int64_t>(index + 1U)) {
Fail("schema-mismatch", "The HDF5 node labels are out of contract.");
}
}
}
template <typename Row>
void RequirePressureInstance(const std::vector<Row>& rows) {
const std::string expected = "PART-1-1";
if (!std::all_of(rows.begin(), rows.end(), [&expected](const Row& row) {
return UppercaseAscii(row.instance_name) == expected;
})) {
Fail("schema-mismatch",
"The pressure reference instance identity is out of contract.");
}
}
std::string Key(const std::string& instance_name,
const std::int64_t source_node_label,
const std::string& component) {
return UppercaseAscii(instance_name) + "\x1f" +
std::to_string(source_node_label) + "\x1f" + component;
}
bool IsSelected(const std::int64_t label,
const std::vector<std::int64_t>& selected_labels) {
return selected_labels.empty() ||
std::find(selected_labels.begin(), selected_labels.end(), label) !=
selected_labels.end();
}
std::map<std::string, double> ProjectCsv(
const std::vector<CsvWideRow>& rows,
const std::array<const char*, 6>& components,
const std::vector<std::size_t>& selected_components,
const std::vector<std::int64_t>& selected_labels) {
std::map<std::string, double> projected;
for (const auto& row : rows) {
if (!IsSelected(row.source_node_label, selected_labels)) {
continue;
}
for (const std::size_t component : selected_components) {
const std::string key =
Key(row.instance_name, row.source_node_label, components[component]);
if (!projected.emplace(key, row.values[component]).second) {
Fail("schema-mismatch", "A CSV component-row identity is duplicated.");
}
}
}
return projected;
}
std::map<std::string, double> ProjectHdf5(
const Hdf5Projection& projection, const std::vector<double>& values,
const std::array<const char*, 6>& components,
const std::vector<std::size_t>& selected_components,
const std::vector<std::int64_t>& selected_labels) {
std::map<std::string, double> projected;
if (values.size() != projection.nodes.size() * 6U) {
Fail("schema-mismatch", "An HDF5 nodal vector has the wrong row count.");
}
for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
if (!IsSelected(projection.nodes[node].source_node_label, selected_labels)) {
continue;
}
for (const std::size_t component : selected_components) {
const std::string key = Key(projection.nodes[node].instance_name,
projection.nodes[node].source_node_label,
components[component]);
const double value = values[node * 6U + component];
if (!std::isfinite(value) || !projected.emplace(key, value).second) {
Fail("schema-mismatch", "An HDF5 component-row identity is invalid.");
}
}
}
return projected;
}
ReferenceToleranceFamilyDecision EvaluateFamily(
const std::string& identity, const std::vector<std::string>& components,
const bool blocking, const std::map<std::string, double>& hdf5,
const std::map<std::string, double>& reference) {
if (hdf5.empty() || hdf5.size() != reference.size()) {
Fail("schema-mismatch", "HDF5 and CSV required row inventories differ.");
}
ReferenceToleranceFamily family{identity, components, blocking, {}};
family.rows.reserve(hdf5.size());
auto hdf5_row = hdf5.begin();
auto reference_row = reference.begin();
for (; hdf5_row != hdf5.end(); ++hdf5_row, ++reference_row) {
if (hdf5_row->first != reference_row->first ||
!std::isfinite(hdf5_row->second) || !std::isfinite(reference_row->second)) {
Fail("schema-mismatch", "HDF5 and CSV required row keys differ or are nonfinite.");
}
family.rows.push_back({hdf5_row->second, reference_row->second});
}
const auto evaluation = ReferenceTolerancePolicy::Evaluate({family});
if (!evaluation.valid || evaluation.families.size() != 1U) {
Fail("schema-mismatch", "The common tolerance policy rejected a finite family.");
}
return evaluation.families.front();
}
void AppendFamily(NonlinearReferenceComparisonReport& report,
const std::string& identity,
const std::vector<std::string>& components,
const NonlinearReferenceDisposition disposition,
const std::map<std::string, double>& hdf5,
const std::map<std::string, double>& reference) {
const bool blocking = disposition == NonlinearReferenceDisposition::kBlocking;
const auto tolerance =
EvaluateFamily(identity, components, blocking, hdf5, reference);
if (tolerance.row_count != hdf5.size()) {
Fail("schema-mismatch", "The tolerance family row count is inconsistent.");
}
report.families.push_back(
{identity, components, disposition, hdf5.size(), tolerance});
if (blocking && !tolerance.passed) {
report.external_blocking_passed = false;
}
}
std::vector<std::int64_t> ExpectedConstrainedLabels(
const NonlinearReferenceCase& reference_case) {
if (reference_case.case_id == "shellnl-s4-final") {
return {kCoarseConstrainedNodeLabels.begin(), kCoarseConstrainedNodeLabels.end()};
}
if (reference_case.case_id == "shellnl-meshfine-s4-final") {
std::vector<std::int64_t> labels{1};
for (std::int64_t label = 3; label <= 9; ++label) {
labels.push_back(label);
}
for (std::int64_t label = 34; label <= 57; ++label) {
labels.push_back(label);
}
for (std::int64_t label = 82; label <= 153; ++label) {
labels.push_back(label);
}
return labels;
}
Fail("schema-mismatch", "A linear case does not define nonlinear constraints.");
}
std::vector<std::int64_t> RequiredConstrainedLabels(
const NonlinearReferenceCase& reference_case,
const Hdf5Projection& projection) {
std::vector<std::int64_t> constrained;
for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
const auto begin = projection.constraint_mask.begin() +
static_cast<std::ptrdiff_t>(node * 6U);
if (std::all_of(begin, begin + 6U,
[](const std::uint8_t value) { return value == 1U; })) {
constrained.push_back(projection.nodes[node].source_node_label);
}
}
std::sort(constrained.begin(), constrained.end());
if (constrained != ExpectedConstrainedLabels(reference_case)) {
Fail("schema-mismatch", "The HDF5 constrained-node selection is out of contract.");
}
return constrained;
}
fesa::Vector3 VectorAt(const std::vector<double>& values,
const std::size_t offset) {
return {values[offset], values[offset + 1U], values[offset + 2U]};
}
fesa::Vector3 SpatialMoment(const fesa::Vector3& phi,
const fesa::Vector3& generalized) {
const Matrix transform = InverseTransposeLeftJacobianSO3(phi);
return {transform(0U, 0U) * generalized.X() +
transform(0U, 1U) * generalized.Y() +
transform(0U, 2U) * generalized.Z(),
transform(1U, 0U) * generalized.X() +
transform(1U, 1U) * generalized.Y() +
transform(1U, 2U) * generalized.Z(),
transform(2U, 0U) * generalized.X() +
transform(2U, 1U) * generalized.Y() +
transform(2U, 2U) * generalized.Z()};
}
fesa::Vector3 MaskedResidual(const Hdf5Projection& projection,
const std::size_t offset) {
return {projection.constraint_mask[offset] == 1U
? projection.residual[offset]
: 0.0,
projection.constraint_mask[offset + 1U] == 1U
? projection.residual[offset + 1U]
: 0.0,
projection.constraint_mask[offset + 2U] == 1U
? projection.residual[offset + 2U]
: 0.0};
}
double InfinityNorm(const Hdf5Projection& projection,
const std::vector<double>& values,
const bool rotations) {
double norm = 0.0;
const std::size_t first_component = rotations ? 3U : 0U;
for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
const std::size_t offset = node * 6U;
for (std::size_t component = 0U; component < 3U; ++component) {
const std::size_t dof = offset + first_component + component;
if (projection.constraint_mask[dof] == 0U) {
norm = (std::max)(norm, std::abs(values[dof]));
}
}
}
return norm;
}
std::size_t FreeRowCount(const Hdf5Projection& projection,
const bool rotations) {
std::size_t count = 0U;
const std::size_t first_component = rotations ? 3U : 0U;
for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
const std::size_t offset = node * 6U;
for (std::size_t component = 0U; component < 3U; ++component) {
if (projection.constraint_mask[offset + first_component + component] == 0U) {
++count;
}
}
}
return count;
}
void AppendInternal(NonlinearReferenceComparisonReport* report,
const NonlinearInternalEquilibriumReport& row) {
report->internal_equilibrium.push_back(row);
if (!row.passed) {
report->internal_blocking_passed = false;
}
}
void AppendFreeResidual(NonlinearReferenceComparisonReport* report,
const std::string& identity,
const Hdf5Projection& projection,
const bool rotations) {
const double residual = InfinityNorm(projection, projection.residual, rotations);
const double internal =
InfinityNorm(projection, projection.internal_force, rotations);
const double external =
InfinityNorm(projection, projection.external_force, rotations);
if (external == 0.0) {
AppendInternal(report, {identity, FreeRowCount(projection, rotations), residual,
0.0, std::nullopt, residual <= 1.0e-3,
"accepted-free-residual-exact-zero-scale"});
return;
}
const double scale = (std::max)(internal, external);
const double metric = residual / scale;
AppendInternal(report, {identity, FreeRowCount(projection, rotations), residual,
scale, metric, metric <= 1.0e-2,
"accepted-free-residual-relative"});
}
void AppendGlobalBalance(NonlinearReferenceComparisonReport* report,
const std::string& identity,
const fesa::Vector3& balance,
const double applied_scale,
const double reaction_scale) {
AppendInternal(report, EvaluateGlobalBalance(identity, balance.Norm(),
applied_scale, reaction_scale));
}
void AppendInternalEquilibrium(NonlinearReferenceComparisonReport* report,
const NonlinearReferenceCase& reference_case,
const Hdf5Projection& projection) {
fesa::Vector3 reaction_force;
fesa::Vector3 reaction_moment;
double applied_force_scale = 0.0;
double reaction_force_scale = 0.0;
double applied_moment_scale = 0.0;
double reaction_moment_scale = 0.0;
for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
const std::size_t offset = node * 6U;
const fesa::Vector3 position =
fesa::Vector3{projection.nodes[node].coordinates} +
VectorAt(projection.displacement, offset);
const fesa::Vector3 phi = VectorAt(projection.displacement, offset + 3U);
const fesa::Vector3 applied_force = VectorAt(projection.external_force, offset);
const fesa::Vector3 constrained_force = MaskedResidual(projection, offset);
const fesa::Vector3 applied_moment =
SpatialMoment(phi, VectorAt(projection.external_force, offset + 3U));
const fesa::Vector3 constrained_moment =
SpatialMoment(phi, MaskedResidual(projection, offset + 3U));
const fesa::Vector3 applied_resultant =
position.Cross(applied_force) + applied_moment;
const fesa::Vector3 constrained_resultant =
position.Cross(constrained_force) + constrained_moment;
reaction_force = reaction_force + constrained_force;
reaction_moment = reaction_moment + constrained_resultant;
applied_force_scale += applied_force.Norm();
reaction_force_scale += constrained_force.Norm();
applied_moment_scale += applied_resultant.Norm();
reaction_moment_scale += constrained_resultant.Norm();
}
report->constrained_reaction_resultant = {
reaction_force.Components(), reaction_moment.Components()};
const std::string prefix = ContractFor(reference_case).family_prefix;
AppendGlobalBalance(report, prefix + "global_force_balance",
{projection.global_equilibrium[0U],
projection.global_equilibrium[1U],
projection.global_equilibrium[2U]},
applied_force_scale, reaction_force_scale);
AppendGlobalBalance(report, prefix + "global_moment_balance",
{projection.global_equilibrium[3U],
projection.global_equilibrium[4U],
projection.global_equilibrium[5U]},
applied_moment_scale, reaction_moment_scale);
AppendFreeResidual(report, prefix + "accepted_free_force_residual",
projection, false);
AppendFreeResidual(report, prefix + "accepted_free_moment_residual",
projection, true);
}
std::filesystem::path DerivedLogPath(const std::filesystem::path& copied_input) {
return copied_input.parent_path() /
(copied_input.stem().string() + "_log.txt");
}
std::filesystem::path DerivedResultsPath(
const std::filesystem::path& copied_input) {
return copied_input.parent_path() /
(copied_input.stem().string() + "_results.h5");
}
NonlinearReferenceComparisonReport MakeReport(
const NonlinearReferenceCase& reference_case,
const std::filesystem::path& hdf5_results_path,
const std::vector<ArtifactSnapshot>& before,
const std::vector<ArtifactSnapshot>& after) {
const auto copied_input = CopiedInputPathFor(hdf5_results_path);
return {reference_case.case_id,
reference_case.original_input_path,
copied_input,
hdf5_results_path,
DerivedLogPath(copied_input),
before,
after,
{},
true,
{},
{},
{},
true,
true};
}
void CompareLinearCase(NonlinearReferenceComparisonReport& report,
const NonlinearReferenceCase& reference_case) {
const Mitc4ReferenceCase legacy_case{
reference_case.case_id, "S4", report.copied_input_path,
reference_case.displacement_csv_path, report.generated_results_hdf5};
const auto legacy = Mitc4ReferenceComparison::Compare(legacy_case);
if (!legacy.HasValue()) {
const auto& diagnostics = legacy.GetStatus().Diagnostics();
const std::string message = diagnostics.empty()
? "The linear control comparison failed."
: diagnostics.front().message;
Fail("schema-mismatch", message);
}
for (const auto& metric : legacy.Value().metrics) {
ReferenceToleranceFamilyDecision tolerance{};
tolerance.identity = metric.blocking ? "linear.translation" : "linear.rotation";
tolerance.components = metric.components;
tolerance.blocking = metric.blocking;
tolerance.reference_scale = metric.reference_scale;
tolerance.near_zero_band = metric.near_zero_band;
tolerance.near_zero_count = metric.near_zero_count;
tolerance.relative_rms = metric.relative_rms;
tolerance.maximum_absolute_error = metric.maximum_absolute_error;
tolerance.worst_row = metric.worst_row;
tolerance.rms_passed = metric.rms_passed;
tolerance.passed = metric.passed;
tolerance.diagnostic_code = metric.diagnostic_code;
tolerance.row_count = metric.row_count;
const auto disposition = metric.blocking
? NonlinearReferenceDisposition::kBlocking
: NonlinearReferenceDisposition::kWarningOnly;
report.families.push_back({tolerance.identity, tolerance.components,
disposition, tolerance.row_count, tolerance});
if (metric.blocking && !metric.passed) {
report.external_blocking_passed = false;
}
}
}
void CompareNonlinearCase(NonlinearReferenceComparisonReport& report,
const NonlinearReferenceCase& reference_case) {
const CaseContract& contract = ContractFor(reference_case);
const auto hdf5 = ReadNonlinearHdf5(reference_case,
report.generated_results_hdf5);
report.copied_input_path = hdf5.copied_input_path;
report.generated_log_path = DerivedLogPath(hdf5.copied_input_path);
RequireExactNodeLabels(hdf5.nodes, contract.node_count);
const auto displacement_rows = ReadCsvWideRows(
reference_case.displacement_csv_path, true, kDisplacementComponents,
contract.node_count);
const auto reaction_rows = ReadCsvWideRows(reference_case.reaction_csv_path,
true, kReactionComponents,
contract.node_count);
RequireExactNodeLabels(displacement_rows, contract.node_count);
RequireExactNodeLabels(reaction_rows, contract.node_count);
const std::vector<std::size_t> translations{0U, 1U, 2U};
const std::vector<std::size_t> rotations{3U, 4U, 5U};
const std::vector<std::size_t> inplane_forces{0U, 1U};
const std::vector<std::size_t> transverse_force{2U};
const std::vector<std::size_t> physical_moments{3U, 4U};
const std::vector<std::int64_t> all_nodes;
const auto constrained = RequiredConstrainedLabels(reference_case, hdf5);
AppendFamily(report, contract.family_prefix + "translation", {"U1", "U2", "U3"},
NonlinearReferenceDisposition::kBlocking,
ProjectHdf5(hdf5, hdf5.displacement, kDisplacementComponents,
translations, all_nodes),
ProjectCsv(displacement_rows, kDisplacementComponents,
translations, all_nodes));
AppendFamily(report, contract.family_prefix + "rotation", {"UR1", "UR2", "UR3"},
NonlinearReferenceDisposition::kWarningOnly,
ProjectHdf5(hdf5, hdf5.displacement, kDisplacementComponents,
rotations, all_nodes),
ProjectCsv(displacement_rows, kDisplacementComponents, rotations,
all_nodes));
AppendFamily(report, contract.family_prefix + "local_reaction_inplane_force",
{"RF1", "RF2"}, NonlinearReferenceDisposition::kReviewOnly,
ProjectHdf5(hdf5, hdf5.reaction, kReactionComponents,
inplane_forces, constrained),
ProjectCsv(reaction_rows, kReactionComponents, inplane_forces,
constrained));
AppendFamily(report, contract.family_prefix + "local_reaction_transverse_force",
{"RF3"}, NonlinearReferenceDisposition::kReviewOnly,
ProjectHdf5(hdf5, hdf5.reaction, kReactionComponents,
transverse_force, constrained),
ProjectCsv(reaction_rows, kReactionComponents, transverse_force,
constrained));
AppendFamily(report, contract.family_prefix + "local_reaction_physical_moment",
{"RM1", "RM2"}, NonlinearReferenceDisposition::kReviewOnly,
ProjectHdf5(hdf5, hdf5.reaction, kReactionComponents,
physical_moments, constrained),
ProjectCsv(reaction_rows, kReactionComponents, physical_moments,
constrained));
AppendInternalEquilibrium(&report, reference_case, hdf5);
}
void ComparePressureCase(NonlinearReferenceComparisonReport& report,
const NonlinearReferenceCase& reference_case) {
const CaseContract& contract = ContractFor(reference_case);
const auto hdf5 =
ReadNonlinearHdf5(reference_case, report.generated_results_hdf5);
report.copied_input_path = hdf5.copied_input_path;
report.generated_log_path = DerivedLogPath(hdf5.copied_input_path);
RequireExactNodeLabels(hdf5.nodes, contract.node_count);
RequirePressureInstance(hdf5.nodes);
const auto displacement_rows =
ReadCsvWideRows(reference_case.displacement_csv_path, true,
kDisplacementComponents, contract.node_count);
RequireExactNodeLabels(displacement_rows, contract.node_count);
RequirePressureInstance(displacement_rows);
const std::vector<std::size_t> translations{0U, 1U, 2U};
const std::vector<std::size_t> rotations{3U, 4U, 5U};
const std::vector<std::int64_t> all_nodes;
AppendFamily(report, "pressure_nodal_translation_blocking",
{"U1", "U2", "U3"}, NonlinearReferenceDisposition::kBlocking,
ProjectHdf5(hdf5, hdf5.displacement, kDisplacementComponents,
translations, all_nodes),
ProjectCsv(displacement_rows, kDisplacementComponents,
translations, all_nodes));
AppendFamily(report, "pressure_nodal_rotation_warning", {"UR1", "UR2", "UR3"},
NonlinearReferenceDisposition::kWarningOnly,
ProjectHdf5(hdf5, hdf5.displacement, kDisplacementComponents,
rotations, all_nodes),
ProjectCsv(displacement_rows, kDisplacementComponents, rotations,
all_nodes));
}
} // namespace
NonlinearInternalEquilibriumReport EvaluateGlobalBalance(
const std::string& family_identity, const double numerator,
const double applied_scale, const double reaction_scale) {
const double scale = (std::max)(applied_scale, reaction_scale);
if (scale == 0.0) {
return {family_identity, 3U, numerator, scale, std::nullopt,
numerator == 0.0, "zero-reference-scale-nonzero-error"};
}
const double metric = numerator / scale;
return {family_identity, 3U, numerator, scale, metric,
metric <= kInternalGlobalBalanceRelativeTolerance,
"global-balance-relative"};
}
Result<std::vector<ArtifactSnapshot>> ArtifactCheck(
const NonlinearReferenceCase& reference_case) {
try {
const auto artifacts = DeclaredArtifacts(reference_case);
std::vector<ArtifactSnapshot> snapshots;
snapshots.reserve(artifacts.size());
for (const auto& artifact : artifacts) {
if (artifact.disposition ==
NonlinearReferenceArtifactDisposition::kRequired) {
snapshots.push_back(Snapshot(artifact));
}
}
return Result<std::vector<ArtifactSnapshot>>::Success(std::move(snapshots));
} catch (const ComparisonFailure& exception) {
return Result<std::vector<ArtifactSnapshot>>::Failure(
FailureStatus(reference_case, exception.Code(), exception.what()));
} catch (const std::exception& exception) {
return Result<std::vector<ArtifactSnapshot>>::Failure(
FailureStatus(reference_case, "reference-artifact-check-failed",
exception.what()));
}
}
Result<std::filesystem::path> CopyInputToExecutionDirectory(
const NonlinearReferenceCase& reference_case) {
try {
const auto check = ArtifactCheck(reference_case);
if (!check.HasValue()) {
return Result<std::filesystem::path>::Failure(check.GetStatus());
}
const std::filesystem::path execution_root =
Normalized(reference_case.execution_directory);
if (!IsDescendantOf(execution_root,
std::filesystem::path{FESA_TEST_BINARY_DIR}) ||
IsDescendantOf(execution_root,
reference_case.original_input_path.parent_path())) {
Fail("invalid-execution-directory",
"Reference execution must use a build-local directory outside reference.");
}
std::error_code error;
std::filesystem::create_directories(execution_root, error);
if (error) {
Fail("copy-input-failed", "The build-local execution directory cannot be created.");
}
static std::atomic<unsigned long long> sequence{0U};
const auto tick = std::chrono::steady_clock::now().time_since_epoch().count();
const std::filesystem::path run_directory =
execution_root / (reference_case.case_id + "-" + std::to_string(tick) +
"-" + std::to_string(sequence.fetch_add(1U)));
if (!std::filesystem::create_directory(run_directory, error) || error) {
Fail("copy-input-failed", "A unique execution directory cannot be created.");
}
const std::filesystem::path copied_input =
run_directory / reference_case.original_input_path.filename();
if (!std::filesystem::copy_file(reference_case.original_input_path,
copied_input,
std::filesystem::copy_options::none, error) ||
error || ReadBytes(copied_input) != ReadBytes(reference_case.original_input_path)) {
Fail("copy-input-failed", "The copied input does not preserve source bytes.");
}
return Result<std::filesystem::path>::Success(copied_input);
} catch (const ComparisonFailure& exception) {
return Result<std::filesystem::path>::Failure(
FailureStatus(reference_case, exception.Code(), exception.what()));
} catch (const std::exception& exception) {
return Result<std::filesystem::path>::Failure(
FailureStatus(reference_case, "copy-input-failed", exception.what()));
}
}
Result<GeneratedRunArtifacts> RunFesaApplicationOnCopiedInput(
const NonlinearReferenceCase& reference_case,
const std::filesystem::path& copied_input_path) {
try {
if (Normalized(copied_input_path) ==
Normalized(reference_case.original_input_path)) {
Fail("original-reference-input-forbidden",
"FesaApplication must never receive the original reference input.");
}
if (!IsDescendantOf(copied_input_path,
Normalized(reference_case.execution_directory)) ||
!std::filesystem::is_regular_file(copied_input_path)) {
Fail("invalid-execution-input", "The execution input is not a copied artifact.");
}
const int exit_code = FesaApplication{}.Run({copied_input_path.string()});
if (exit_code != 0) {
Fail("fesa-application-failed",
"FesaApplication did not produce an accepted reference run.");
}
const auto results = DerivedResultsPath(copied_input_path);
const auto log = DerivedLogPath(copied_input_path);
if (!std::filesystem::is_regular_file(results) ||
!std::filesystem::is_regular_file(log)) {
Fail("generated-artifact-missing",
"The copied execution did not produce sibling result and log artifacts.");
}
return Result<GeneratedRunArtifacts>::Success({copied_input_path, results, log});
} catch (const ComparisonFailure& exception) {
return Result<GeneratedRunArtifacts>::Failure(
FailureStatus(reference_case, exception.Code(), exception.what()));
} catch (const std::exception& exception) {
return Result<GeneratedRunArtifacts>::Failure(
FailureStatus(reference_case, "fesa-application-failed", exception.what()));
}
}
Result<NonlinearReferenceComparisonReport> CompareHdf5ToReferenceCsv(
const NonlinearReferenceCase& reference_case,
const std::filesystem::path& hdf5_results_path) {
try {
const auto before = ArtifactCheck(reference_case);
if (!before.HasValue()) {
return Result<NonlinearReferenceComparisonReport>::Failure(
before.GetStatus());
}
if (!std::filesystem::is_regular_file(hdf5_results_path)) {
Fail("schema-mismatch", "The generated HDF5 result is missing.");
}
auto report =
MakeReport(reference_case, hdf5_results_path, before.Value(), {});
if (reference_case.case_id == "shellnl-meshfine-pressure-s4-follower") {
ComparePressureCase(report, reference_case);
} else if (ContractFor(reference_case).nonlinear) {
CompareNonlinearCase(report, reference_case);
} else if (reference_case.case_id == "shell-s4-linear-limit") {
CompareLinearCase(report, reference_case);
} else {
Fail("schema-mismatch", "The requested comparison case is unknown.");
}
const auto after = ArtifactCheck(reference_case);
if (!after.HasValue()) {
return Result<NonlinearReferenceComparisonReport>::Failure(
after.GetStatus());
}
report.post_run_artifacts = after.Value();
report.blocking_passed =
report.external_blocking_passed && report.internal_blocking_passed;
return Result<NonlinearReferenceComparisonReport>::Success(
std::move(report));
} catch (const ComparisonFailure& exception) {
return Result<NonlinearReferenceComparisonReport>::Failure(
FailureStatus(reference_case, exception.Code(), exception.what()));
} catch (const std::exception& exception) {
return Result<NonlinearReferenceComparisonReport>::Failure(
FailureStatus(reference_case, "reference-comparison-failed", exception.what()));
}
}
NonlinearReferenceClassification ClassifyReferenceComparison(
const NonlinearReferenceComparisonReport& report) {
const bool blocking_precheck =
std::any_of(report.precheck_findings.begin(), report.precheck_findings.end(),
[](const NonlinearReferencePrecheckFinding& finding) {
return finding.disposition ==
NonlinearReferenceDisposition::kBlocking;
});
if (!report.blocking_passed || blocking_precheck) {
return NonlinearReferenceClassification::kFailed;
}
const bool warning_only_difference =
std::any_of(report.families.begin(), report.families.end(),
[](const NonlinearReferenceFamilyReport& family) {
return family.disposition ==
NonlinearReferenceDisposition::kWarningOnly &&
!family.tolerance.passed;
});
const bool warning_precheck =
std::any_of(report.precheck_findings.begin(), report.precheck_findings.end(),
[](const NonlinearReferencePrecheckFinding& finding) {
return finding.disposition ==
NonlinearReferenceDisposition::kWarningOnly;
});
return warning_only_difference || warning_precheck
? NonlinearReferenceClassification::kWarningOnlyDifferences
: NonlinearReferenceClassification::kPassed;
}
NonlinearReferenceComparisonReport BuildReferenceComparisonReport(
const NonlinearReferenceCase& reference_case,
const std::filesystem::path& hdf5_results_path) {
NonlinearReferenceComparisonReport report{};
report.case_id = reference_case.case_id;
report.original_input_path = reference_case.original_input_path;
report.generated_results_hdf5 = hdf5_results_path;
report.blocking_passed = false;
const auto artifacts = ArtifactCheck(reference_case);
if (artifacts.HasValue()) {
report.pre_run_artifacts = artifacts.Value();
report.post_run_artifacts = artifacts.Value();
}
if (hdf5_results_path.has_filename()) {
try {
report.copied_input_path = CopiedInputPathFor(hdf5_results_path);
report.generated_log_path = DerivedLogPath(report.copied_input_path);
} catch (const ComparisonFailure&) {
// This builder is non-throwing by contract; comparison reports surface
// the malformed result path through CompareHdf5ToReferenceCsv instead.
}
}
return report;
}
} // namespace fesa::test
Reference in New Issue View Git Blame Copy Permalink