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

1421 lines
56 KiB
C++
Raw Blame History

#include "reference_comparison.h"
#include <hdf5.h>
#include <algorithm>
#include <array>
#include <cctype>
#include <cerrno>
#include <charconv>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iterator>
#include <limits>
#include <locale>
#include <sstream>
#include <stdexcept>
#include <string>
#include <system_error>
#include <utility>
#include <vector>
#include "fesa/analysis/analysis_model.h"
#include "fesa/assembly/load_assembler.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/io/abaqus/domain_mapper.h"
#include "fesa/io/abaqus/input_reader.h"
#include "fesa/results/result_recovery.h"
namespace fesa::test {
namespace {
constexpr const char* kModelId = "cantilever-beam-b33";
constexpr const char* kStepName = "Step-1";
constexpr std::size_t kFrameIndex = 0U;
constexpr const char* kFrameText = "Increment 1: Step Time = 1.000";
constexpr const char* kInputName = "cantilever beam.inp";
constexpr const char* kDisplacementName = "cantilever beam displacements.csv";
constexpr const char* kReactionName = "cantilever beam reactions.csv";
constexpr const char* kSectionName = "cantilever beam elemental forces.csv";
constexpr const char* kDisplacementPath =
"/steps/Step-1/frames/0/nodal/displacement";
constexpr const char* kReactionPath = "/steps/Step-1/frames/0/nodal/reaction";
constexpr const char* kSectionPath =
"/steps/Step-1/frames/0/element/section_resultant";
constexpr const char* kStressPath = "/steps/Step-1/frames/0/element/stress_s11";
constexpr double kKinematicFloor = 1.0e-9;
constexpr double kForceMomentFloor = 1.0e-3;
constexpr double kRelativeCoefficient = 1.0e-6;
class ComparisonFailure final : public std::runtime_error {
public:
ComparisonFailure(std::string code, std::string message)
: std::runtime_error{std::move(message)}, code_{std::move(code)} {}
const std::string& Code() const noexcept { return code_; }
private:
std::string code_;
};
[[noreturn]] void Fail(const std::string& code, const std::string& message) {
throw ComparisonFailure{code, message};
}
Status ComparisonFailureStatus(const std::string& code,
const std::string& message) {
return Status::Failure(FailureCategory::kModel,
{{Severity::kError, code, {}, "", kModelId, message}});
}
std::string Trim(const std::string& value) {
const auto is_space = [](const unsigned char character) {
return std::isspace(character) != 0;
};
const auto begin =
std::find_if_not(value.begin(), value.end(), [&](const char character) {
return is_space(static_cast<unsigned char>(character));
});
const auto end =
std::find_if_not(value.rbegin(), value.rend(), [&](const char character) {
return is_space(static_cast<unsigned char>(character));
}).base();
return begin < end ? std::string{begin, end} : std::string{};
}
std::string CollapseWhitespace(const std::string& value) {
std::string result;
bool pending_space = false;
for (const char character : Trim(value)) {
if (std::isspace(static_cast<unsigned char>(character)) != 0) {
pending_space = !result.empty();
} else {
if (pending_space) {
result.push_back(' ');
}
result.push_back(character);
pending_space = false;
}
}
return result;
}
std::string AsciiLower(std::string value) {
std::transform(value.begin(), value.end(), value.begin(),
[](const char character) {
if (character >= 'A' && character <= 'Z') {
return static_cast<char>(character - 'A' + 'a');
}
return character;
});
return value;
}
std::vector<std::string> SplitCsvLine(const std::string& line) {
std::vector<std::string> fields;
std::size_t start = 0U;
while (true) {
const std::size_t comma = line.find(',', start);
fields.push_back(Trim(line.substr(start, comma - start)));
if (comma == std::string::npos) {
break;
}
start = comma + 1U;
}
return fields;
}
std::int64_t ParsePositiveLabel(const std::string& field) {
std::int64_t value = 0;
const char* const begin = field.data();
const char* const end = begin + field.size();
const auto parsed = std::from_chars(begin, end, value);
if (parsed.ec != std::errc{} || parsed.ptr != end || value <= 0) {
Fail("schema-mismatch", "A CSV or HDF5 source-node label is invalid.");
}
return value;
}
double ParseFiniteDouble(const std::string& field) {
if (field.empty()) {
Fail("schema-mismatch", "A reference numeric field is empty.");
}
errno = 0;
char* end = nullptr;
const double value = std::strtod(field.c_str(), &end);
if (errno == ERANGE || end == field.c_str() || end == nullptr ||
*end != '\0' || !std::isfinite(value)) {
Fail("schema-mismatch",
"A reference numeric field is nonfinite or invalid.");
}
return value;
}
struct WideReferenceRow {
std::string instance_name;
std::int64_t source_node_label;
std::vector<double> values;
};
struct ReferenceTable {
std::vector<WideReferenceRow> rows;
};
ReferenceTable ReadReferenceCsv(
const std::filesystem::path& path,
const std::vector<std::string>& expected_header) {
std::ifstream stream{path};
if (!stream) {
Fail("needs-reference-artifacts", "An approved reference CSV is missing.");
}
std::string line;
if (!std::getline(stream, line)) {
Fail("schema-mismatch", "An approved reference CSV is empty.");
}
if (!line.empty() && line.back() == '\r') {
line.pop_back();
}
if (SplitCsvLine(line) != expected_header) {
Fail("schema-mismatch", "An approved reference CSV header is not exact.");
}
ReferenceTable table;
while (std::getline(stream, line)) {
if (!line.empty() && line.back() == '\r') {
line.pop_back();
}
if (line.empty()) {
Fail("schema-mismatch", "Blank reference CSV rows are not allowed.");
}
const auto fields = SplitCsvLine(line);
if (fields.size() != expected_header.size() ||
CollapseWhitespace(fields[0U]) != kFrameText || fields[1U].empty()) {
Fail("schema-mismatch",
"A reference CSV row has invalid schema or frame identity.");
}
WideReferenceRow row{};
row.instance_name = fields[1U];
row.source_node_label = ParsePositiveLabel(fields[2U]);
row.values.reserve(fields.size() - 3U);
for (std::size_t field = 3U; field < fields.size(); ++field) {
row.values.push_back(ParseFiniteDouble(fields[field]));
}
const auto duplicate = std::find_if(
table.rows.begin(), table.rows.end(),
[&](const WideReferenceRow& existing) {
return AsciiLower(existing.instance_name) ==
AsciiLower(row.instance_name) &&
existing.source_node_label == row.source_node_label;
});
if (duplicate != table.rows.end()) {
Fail("schema-mismatch", "A reference CSV row identity is duplicated.");
}
table.rows.push_back(std::move(row));
}
if (table.rows.empty()) {
Fail("schema-mismatch", "An approved reference CSV has no data rows.");
}
return table;
}
void RequireExactArtifactInventory(
const std::filesystem::path& legacy_directory) {
std::error_code error;
if (!std::filesystem::is_directory(legacy_directory, error) || error) {
Fail("needs-reference-artifacts",
"The approved legacy directory is missing.");
}
std::vector<std::string> names;
for (std::filesystem::directory_iterator iterator{legacy_directory, error},
end;
iterator != end && !error; iterator.increment(error)) {
if (!iterator->is_regular_file(error) || error) {
Fail("needs-reference-artifacts",
"The legacy bundle contains a non-file entry.");
}
names.push_back(iterator->path().filename().string());
}
if (error) {
Fail("needs-reference-artifacts", "The legacy bundle cannot be inspected.");
}
std::sort(names.begin(), names.end());
std::vector<std::string> expected = {kDisplacementName, kInputName,
kReactionName, kSectionName};
std::sort(expected.begin(), expected.end());
if (names != expected) {
Fail("needs-reference-artifacts",
"The legacy bundle inventory is not exact.");
}
}
Domain ReadApprovedDomain(const std::filesystem::path& input_path) {
AbaqusInputReader reader;
auto parsed = reader.Read(input_path);
if (!parsed.HasValue()) {
Fail("needs-reference-artifacts",
"The approved reference input cannot be parsed.");
}
AbaqusDomainMapper mapper;
auto domain = mapper.Map(parsed.Value());
if (!domain.HasValue()) {
Fail("needs-reference-artifacts",
"The approved reference input is not the required B33 model.");
}
return std::move(domain.Value());
}
class Hdf5Handle {
public:
using Closer = herr_t (*)(hid_t);
Hdf5Handle() = default;
Hdf5Handle(const hid_t value, Closer closer)
: value_{value}, closer_{closer} {}
Hdf5Handle(const Hdf5Handle&) = delete;
Hdf5Handle& operator=(const Hdf5Handle&) = delete;
Hdf5Handle(Hdf5Handle&& other) noexcept
: value_{other.value_}, closer_{other.closer_} {
other.value_ = -1;
other.closer_ = nullptr;
}
Hdf5Handle& operator=(Hdf5Handle&& other) noexcept {
if (this != &other) {
Reset();
value_ = other.value_;
closer_ = other.closer_;
other.value_ = -1;
other.closer_ = nullptr;
}
return *this;
}
~Hdf5Handle() { Reset(); }
hid_t Get() const noexcept { return value_; }
private:
void Reset() noexcept {
if (value_ >= 0 && closer_ != nullptr) {
(void)closer_(value_);
}
value_ = -1;
closer_ = nullptr;
}
hid_t value_{-1};
Closer closer_{nullptr};
};
class Hdf5ErrorSilencer {
public:
Hdf5ErrorSilencer() {
if (H5Eget_auto2(H5E_DEFAULT, &callback_, &client_data_) >= 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_, 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) noexcept
: memory_type_{memory_type}, data_space_{data_space}, data_{data} {}
Hdf5VlenReclaimer(const Hdf5VlenReclaimer&) = delete;
Hdf5VlenReclaimer& operator=(const Hdf5VlenReclaimer&) = delete;
~Hdf5VlenReclaimer() {
if (active_) {
(void)H5Dvlen_reclaim(memory_type_, data_space_, H5P_DEFAULT, data_);
}
}
herr_t Reclaim() noexcept {
active_ = false;
return H5Dvlen_reclaim(memory_type_, data_space_, H5P_DEFAULT, data_);
}
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 char* path) {
return {RequireId(H5Dopen2(file, path, H5P_DEFAULT),
"A required HDF5 dataset is missing."),
H5Dclose};
}
std::vector<hsize_t> DatasetDimensions(const hid_t dataset) {
Hdf5Handle space{
RequireId(H5Dget_space(dataset), "Unable to inspect an HDF5 dataspace."),
H5Sclose};
const int rank = H5Sget_simple_extent_ndims(space.Get());
if (rank < 0) {
Fail("schema-mismatch", "Unable to inspect an HDF5 dataset rank.");
}
std::vector<hsize_t> dimensions(static_cast<std::size_t>(rank));
if (rank > 0) {
RequireHdf(
H5Sget_simple_extent_dims(space.Get(), dimensions.data(), nullptr),
"Unable to inspect HDF5 dataset dimensions.");
}
return dimensions;
}
std::string ReadStringAttribute(const hid_t object, const char* name) {
Hdf5Handle attribute{
RequireId(H5Aopen(object, name, H5P_DEFAULT),
"A required HDF5 string attribute is missing."),
H5Aclose};
Hdf5Handle type{RequireId(H5Aget_type(attribute.Get()),
"Unable to inspect an HDF5 string attribute."),
H5Tclose};
if (H5Tget_class(type.Get()) != H5T_STRING ||
H5Tis_variable_str(type.Get()) <= 0 ||
H5Tget_cset(type.Get()) != H5T_CSET_UTF8) {
Fail("schema-mismatch", "An HDF5 string attribute has the wrong type.");
}
char* raw = nullptr;
RequireHdf(H5Aread(attribute.Get(), type.Get(), &raw),
"Unable to read an HDF5 string attribute.");
if (raw == nullptr) {
Fail("schema-mismatch", "An HDF5 string attribute is null.");
}
const std::string value{raw};
RequireHdf(H5free_memory(raw), "Unable to release HDF5 string memory.");
return value;
}
std::uint64_t ReadUint64Attribute(const hid_t object, const char* name) {
Hdf5Handle attribute{
RequireId(H5Aopen(object, name, H5P_DEFAULT),
"A required HDF5 integer attribute is missing."),
H5Aclose};
Hdf5Handle type{RequireId(H5Aget_type(attribute.Get()),
"Unable to inspect an HDF5 integer attribute."),
H5Tclose};
if (H5Tget_class(type.Get()) != H5T_INTEGER ||
H5Tget_size(type.Get()) != sizeof(std::uint64_t) ||
H5Tget_sign(type.Get()) != H5T_SGN_NONE) {
Fail("schema-mismatch", "An HDF5 integer attribute has the wrong type.");
}
std::uint64_t value = 0U;
RequireHdf(H5Aread(attribute.Get(), H5T_NATIVE_UINT64, &value),
"Unable to read an HDF5 integer attribute.");
return value;
}
void RequireStringAttribute(const hid_t object, const char* name,
const char* expected) {
if (ReadStringAttribute(object, name) != expected) {
Fail("schema-mismatch", "An HDF5 string attribute has the wrong value.");
}
}
void RequireResultAttributes(const hid_t dataset, const char* components,
const char* units, const char* coordinate_system,
const char* location) {
RequireStringAttribute(dataset, "component_names", components);
RequireStringAttribute(dataset, "component_unit_dimensions", units);
RequireStringAttribute(dataset, "coordinate_system", coordinate_system);
RequireStringAttribute(dataset, "location", location);
RequireStringAttribute(dataset, "step_name", kStepName);
if (ReadUint64Attribute(dataset, "frame_index") != kFrameIndex) {
Fail("schema-mismatch", "An HDF5 result has the wrong frame identity.");
}
}
std::vector<double> ReadDoubleDataset(
const hid_t file, const char* path,
const std::vector<hsize_t>& expected_dimensions, const char* components,
const char* units, const char* coordinate_system, const char* location) {
auto dataset = OpenDataset(file, path);
if (DatasetDimensions(dataset.Get()) != expected_dimensions) {
Fail("schema-mismatch", "An HDF5 result dataset has the wrong shape.");
}
Hdf5Handle type{RequireId(H5Dget_type(dataset.Get()),
"Unable to inspect an HDF5 result type."),
H5Tclose};
if (H5Tget_class(type.Get()) != H5T_FLOAT ||
H5Tget_size(type.Get()) != sizeof(double) ||
H5Tequal(type.Get(), H5T_IEEE_F64LE) <= 0) {
Fail("schema-mismatch", "An HDF5 result dataset is not float64.");
}
RequireResultAttributes(dataset.Get(), components, units, coordinate_system,
location);
std::size_t count = 1U;
for (const hsize_t dimension : expected_dimensions) {
if (dimension > (std::numeric_limits<std::size_t>::max)() / count) {
Fail("schema-mismatch", "An HDF5 result shape overflows size_t.");
}
count *= static_cast<std::size_t>(dimension);
}
std::vector<double> values(count);
if (!values.empty()) {
RequireHdf(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
H5P_DEFAULT, values.data()),
"Unable to read an HDF5 result dataset.");
}
if (!std::all_of(values.begin(), values.end(),
[](const double value) { return std::isfinite(value); })) {
Fail("schema-mismatch", "An HDF5 comparison value is nonfinite.");
}
return values;
}
void RequireCompoundMembers(const hid_t dataset,
const std::vector<std::string>& expected) {
Hdf5Handle type{RequireId(H5Dget_type(dataset),
"Unable to inspect an HDF5 compound type."),
H5Tclose};
if (H5Tget_class(type.Get()) != H5T_COMPOUND ||
H5Tget_nmembers(type.Get()) != static_cast<int>(expected.size())) {
Fail("schema-mismatch", "An HDF5 compound dataset has the wrong schema.");
}
for (std::size_t index = 0U; index < expected.size(); ++index) {
char* raw = H5Tget_member_name(type.Get(), static_cast<unsigned>(index));
if (raw == nullptr) {
Fail("schema-mismatch", "Unable to inspect an HDF5 member name.");
}
const std::string actual{raw};
RequireHdf(H5free_memory(raw), "Unable to release an HDF5 member name.");
if (actual != expected[index]) {
Fail("schema-mismatch", "An HDF5 compound member has the wrong name.");
}
}
}
Hdf5Handle MakeVariableStringType() {
Hdf5Handle type{
RequireId(H5Tcopy(H5T_C_S1), "Unable to create an HDF5 string type."),
H5Tclose};
RequireHdf(H5Tset_size(type.Get(), H5T_VARIABLE),
"Unable to size an HDF5 string type.");
RequireHdf(H5Tset_cset(type.Get(), H5T_CSET_UTF8),
"Unable to configure an HDF5 string type.");
return type;
}
struct NodeMemoryRow {
std::uint64_t internal_node_id;
char* instance_name;
char* source_label;
double coordinates[3];
};
struct ElementMemoryRow {
std::uint64_t internal_element_id;
char* instance_name;
char* source_label;
std::uint64_t node_internal_ids[2];
double local_axes[9];
};
struct HdfNode {
std::uint64_t internal_node_id;
std::string instance_name;
std::int64_t source_node_label;
std::string source_node_label_text;
std::array<double, 3> coordinates;
};
struct HdfElement {
std::uint64_t internal_element_id;
std::string instance_name;
std::int64_t source_element_label;
std::string source_element_label_text;
std::array<std::uint64_t, 2> node_internal_ids;
std::array<double, 9> local_axes;
};
std::vector<HdfNode> ReadNodeRows(const hid_t file) {
auto dataset = OpenDataset(file, "/model/nodes");
const auto dimensions = DatasetDimensions(dataset.Get());
if (dimensions.size() != 1U || dimensions[0U] == 0U) {
Fail("schema-mismatch", "The HDF5 node table has the wrong shape.");
}
RequireCompoundMembers(dataset.Get(), {"internal_node_id", "instance_name",
"source_label", "coordinates"});
RequireStringAttribute(dataset.Get(), "coordinate_system",
"global-cartesian");
RequireStringAttribute(dataset.Get(), "units_label", "length");
auto string_type = MakeVariableStringType();
const hsize_t coordinate_dimensions[1] = {3U};
Hdf5Handle coordinate_type{
RequireId(H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinate_dimensions),
"Unable to create the node coordinate memory type."),
H5Tclose};
Hdf5Handle memory_type{
RequireId(H5Tcreate(H5T_COMPOUND, sizeof(NodeMemoryRow)),
"Unable to create the node memory type."),
H5Tclose};
RequireHdf(
H5Tinsert(memory_type.Get(), "internal_node_id",
HOFFSET(NodeMemoryRow, internal_node_id), H5T_NATIVE_UINT64),
"Unable to define the node ID memory field.");
RequireHdf(
H5Tinsert(memory_type.Get(), "instance_name",
HOFFSET(NodeMemoryRow, instance_name), string_type.Get()),
"Unable to define the node instance memory field.");
RequireHdf(H5Tinsert(memory_type.Get(), "source_label",
HOFFSET(NodeMemoryRow, source_label), string_type.Get()),
"Unable to define the node label memory field.");
RequireHdf(
H5Tinsert(memory_type.Get(), "coordinates",
HOFFSET(NodeMemoryRow, coordinates), coordinate_type.Get()),
"Unable to define the node coordinate memory field.");
std::vector<NodeMemoryRow> raw(static_cast<std::size_t>(dimensions[0U]));
Hdf5Handle space{RequireId(H5Dget_space(dataset.Get()),
"Unable to reopen the node dataspace."),
H5Sclose};
RequireHdf(H5Dread(dataset.Get(), memory_type.Get(), H5S_ALL, H5S_ALL,
H5P_DEFAULT, raw.data()),
"Unable to read the HDF5 node table.");
Hdf5VlenReclaimer strings{memory_type.Get(), space.Get(), raw.data()};
std::vector<HdfNode> rows;
rows.reserve(raw.size());
for (const auto& row : raw) {
if (row.instance_name == nullptr || row.source_label == nullptr) {
Fail("schema-mismatch", "An HDF5 node identity is null.");
}
const std::array<double, 3> coordinates = {
row.coordinates[0U], row.coordinates[1U], row.coordinates[2U]};
if (!std::all_of(coordinates.begin(), coordinates.end(),
[](const double value) { return std::isfinite(value); })) {
Fail("schema-mismatch", "An HDF5 node coordinate is nonfinite.");
}
rows.push_back({row.internal_node_id, row.instance_name,
ParsePositiveLabel(row.source_label), row.source_label,
coordinates});
}
RequireHdf(strings.Reclaim(), "Unable to reclaim HDF5 node strings.");
return rows;
}
std::vector<HdfElement> ReadElementRows(const hid_t file) {
auto dataset = OpenDataset(file, "/model/elements");
const auto dimensions = DatasetDimensions(dataset.Get());
if (dimensions.size() != 1U || dimensions[0U] == 0U) {
Fail("schema-mismatch", "The HDF5 element table has the wrong shape.");
}
RequireCompoundMembers(
dataset.Get(), {"internal_element_id", "instance_name", "source_label",
"node_internal_ids", "local_axes"});
RequireStringAttribute(dataset.Get(), "formulation",
"B33-3D-Euler-Bernoulli");
auto string_type = MakeVariableStringType();
const hsize_t connectivity_dimensions[1] = {2U};
const hsize_t axes_dimensions[2] = {3U, 3U};
Hdf5Handle connectivity_type{
RequireId(H5Tarray_create2(H5T_NATIVE_UINT64, 1, connectivity_dimensions),
"Unable to create the connectivity memory type."),
H5Tclose};
Hdf5Handle axes_type{
RequireId(H5Tarray_create2(H5T_NATIVE_DOUBLE, 2, axes_dimensions),
"Unable to create the local-axis memory type."),
H5Tclose};
Hdf5Handle memory_type{
RequireId(H5Tcreate(H5T_COMPOUND, sizeof(ElementMemoryRow)),
"Unable to create the element memory type."),
H5Tclose};
RequireHdf(H5Tinsert(memory_type.Get(), "internal_element_id",
HOFFSET(ElementMemoryRow, internal_element_id),
H5T_NATIVE_UINT64),
"Unable to define the element ID memory field.");
RequireHdf(
H5Tinsert(memory_type.Get(), "instance_name",
HOFFSET(ElementMemoryRow, instance_name), string_type.Get()),
"Unable to define the element instance memory field.");
RequireHdf(
H5Tinsert(memory_type.Get(), "source_label",
HOFFSET(ElementMemoryRow, source_label), string_type.Get()),
"Unable to define the element label memory field.");
RequireHdf(H5Tinsert(memory_type.Get(), "node_internal_ids",
HOFFSET(ElementMemoryRow, node_internal_ids),
connectivity_type.Get()),
"Unable to define the connectivity memory field.");
RequireHdf(H5Tinsert(memory_type.Get(), "local_axes",
HOFFSET(ElementMemoryRow, local_axes), axes_type.Get()),
"Unable to define the local-axis memory field.");
std::vector<ElementMemoryRow> raw(static_cast<std::size_t>(dimensions[0U]));
Hdf5Handle space{RequireId(H5Dget_space(dataset.Get()),
"Unable to reopen the element dataspace."),
H5Sclose};
RequireHdf(H5Dread(dataset.Get(), memory_type.Get(), H5S_ALL, H5S_ALL,
H5P_DEFAULT, raw.data()),
"Unable to read the HDF5 element table.");
Hdf5VlenReclaimer strings{memory_type.Get(), space.Get(), raw.data()};
std::vector<HdfElement> rows;
rows.reserve(raw.size());
for (const auto& row : raw) {
if (row.instance_name == nullptr || row.source_label == nullptr) {
Fail("schema-mismatch", "An HDF5 element identity is null.");
}
HdfElement converted{};
converted.internal_element_id = row.internal_element_id;
converted.instance_name = row.instance_name;
converted.source_element_label = ParsePositiveLabel(row.source_label);
converted.source_element_label_text = row.source_label;
converted.node_internal_ids = {row.node_internal_ids[0U],
row.node_internal_ids[1U]};
std::copy(std::begin(row.local_axes), std::end(row.local_axes),
converted.local_axes.begin());
if (!std::all_of(converted.local_axes.begin(), converted.local_axes.end(),
[](const double value) { return std::isfinite(value); })) {
Fail("schema-mismatch", "An HDF5 local axis is nonfinite.");
}
rows.push_back(std::move(converted));
}
RequireHdf(strings.Reclaim(), "Unable to reclaim HDF5 element strings.");
return rows;
}
void RequireFiniteStress(const hid_t file) {
auto dataset = OpenDataset(file, kStressPath);
const auto dimensions = DatasetDimensions(dataset.Get());
if (dimensions.size() != 1U || dimensions[0U] == 0U) {
Fail("schema-mismatch", "The mandatory stress dataset has no rows.");
}
RequireCompoundMembers(dataset.Get(),
{"internal_element_id", "gauss_point_index",
"section_point_index", "x1", "x2", "source", "S11"});
RequireResultAttributes(dataset.Get(), "S11", "force/length^2", "beam-local",
"section-point");
struct StressValue {
double s11;
};
Hdf5Handle memory_type{RequireId(H5Tcreate(H5T_COMPOUND, sizeof(StressValue)),
"Unable to create a stress memory type."),
H5Tclose};
RequireHdf(H5Tinsert(memory_type.Get(), "S11", HOFFSET(StressValue, s11),
H5T_NATIVE_DOUBLE),
"Unable to define the stress memory field.");
std::vector<StressValue> values(static_cast<std::size_t>(dimensions[0U]));
RequireHdf(H5Dread(dataset.Get(), memory_type.Get(), H5S_ALL, H5S_ALL,
H5P_DEFAULT, values.data()),
"Unable to read the stress dataset.");
if (!std::all_of(values.begin(), values.end(), [](const StressValue& value) {
return std::isfinite(value.s11);
})) {
Fail("schema-mismatch", "The mandatory stress dataset is nonfinite.");
}
}
std::array<double, 9> ExpectedLocalAxes(const Domain& domain,
const EulerBeam3DDefinition& element) {
const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates;
const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates;
const auto& guide = domain.Sections()[element.section_index].first_axis;
const std::array<double, 3> delta = {
second[0U] - first[0U], second[1U] - first[1U], second[2U] - first[2U]};
const double length = std::hypot(delta[0U], delta[1U], delta[2U]);
const std::array<double, 3> x = {delta[0U] / length, delta[1U] / length,
delta[2U] / length};
const double projection =
guide[0U] * x[0U] + guide[1U] * x[1U] + guide[2U] * x[2U];
const std::array<double, 3> y_trial = {guide[0U] - projection * x[0U],
guide[1U] - projection * x[1U],
guide[2U] - projection * x[2U]};
const double y_norm = std::hypot(y_trial[0U], y_trial[1U], y_trial[2U]);
const std::array<double, 3> y = {y_trial[0U] / y_norm, y_trial[1U] / y_norm,
y_trial[2U] / y_norm};
const std::array<double, 3> z = {x[1U] * y[2U] - x[2U] * y[1U],
x[2U] * y[0U] - x[0U] * y[2U],
x[0U] * y[1U] - x[1U] * y[0U]};
return {x[0U], x[1U], x[2U], y[0U], y[1U], y[2U], z[0U], z[1U], z[2U]};
}
struct HdfProjection {
std::vector<HdfNode> nodes;
std::vector<HdfElement> elements;
std::vector<double> displacement;
std::vector<double> reaction;
std::vector<double> section_resultants;
};
HdfProjection ReadHdfProjection(const std::filesystem::path& results,
const std::filesystem::path& input,
const Domain& domain) {
std::error_code error;
if (!std::filesystem::is_regular_file(results, error) || error) {
Fail("needs-solver-results",
"The authoritative FESA results.h5 is missing.");
}
Hdf5ErrorSilencer silence;
if (H5Fis_hdf5(results.string().c_str()) <= 0) {
Fail("schema-mismatch", "The solver result is not an HDF5 file.");
}
Hdf5Handle file{
RequireId(H5Fopen(results.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT),
"Unable to open the solver HDF5 file read-only."),
H5Fclose};
Hdf5Handle metadata{RequireId(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT),
"The HDF5 metadata group is missing."),
H5Gclose};
if (ReadUint64Attribute(metadata.Get(), "schema_version") != 0U ||
ReadUint64Attribute(metadata.Get(), "frame_index") != kFrameIndex) {
Fail("schema-mismatch", "The HDF5 schema or frame version is wrong.");
}
RequireStringAttribute(metadata.Get(), "feature_id",
"linear-static-3d-euler-beam");
RequireStringAttribute(metadata.Get(), "unit_system_label",
"user-consistent-unspecified");
RequireStringAttribute(metadata.Get(), "coordinate_convention",
"global-cartesian; beam-local=(t,n1,t-cross-n1)");
RequireStringAttribute(metadata.Get(), "element_formulation",
"B33-3D-Euler-Bernoulli");
RequireStringAttribute(metadata.Get(), "step_name", kStepName);
const std::string source_identity =
ReadStringAttribute(metadata.Get(), "source_input_identity");
const std::string normalized_input =
std::filesystem::absolute(input).lexically_normal().generic_u8string();
const std::string expected_identity =
"path=" + normalized_input +
";content_identity=" + domain.SourceContentIdentity();
if (source_identity != expected_identity) {
Fail("schema-mismatch", "The HDF5 source-input identity is inconsistent.");
}
HdfProjection projection{};
projection.nodes = ReadNodeRows(file.Get());
projection.elements = ReadElementRows(file.Get());
if (projection.nodes.size() != domain.Nodes().size() ||
projection.elements.size() != domain.Elements().size()) {
Fail("schema-mismatch",
"HDF5 model identity counts do not match the input.");
}
for (std::size_t node = 0U; node < projection.nodes.size(); ++node) {
const auto& actual = projection.nodes[node];
const auto& expected = domain.Nodes()[node];
if (actual.internal_node_id != node ||
actual.instance_name != expected.source_id.instance_name ||
actual.source_node_label != expected.source_id.source_label ||
actual.source_node_label_text != expected.source_id.source_label_text ||
actual.coordinates != expected.coordinates) {
Fail("schema-mismatch",
"An HDF5 node identity does not match the input.");
}
}
for (std::size_t element = 0U; element < projection.elements.size();
++element) {
const auto& actual = projection.elements[element];
const auto& expected = domain.Elements()[element];
if (actual.internal_element_id != element ||
actual.instance_name != expected.source_id.instance_name ||
actual.source_element_label != expected.source_id.source_label ||
actual.source_element_label_text !=
expected.source_id.source_label_text ||
actual.node_internal_ids[0U] != expected.node_indices[0U] ||
actual.node_internal_ids[1U] != expected.node_indices[1U]) {
Fail("schema-mismatch",
"An HDF5 element identity does not match the input.");
}
const auto axes = ExpectedLocalAxes(domain, expected);
for (std::size_t component = 0U; component < axes.size(); ++component) {
if (std::abs(actual.local_axes[component] - axes[component]) > 1.0e-12) {
Fail("schema-mismatch", "An HDF5 local axis does not match the input.");
}
}
}
const hsize_t node_count = static_cast<hsize_t>(projection.nodes.size());
const hsize_t element_count =
static_cast<hsize_t>(projection.elements.size());
projection.displacement = ReadDoubleDataset(
file.Get(), kDisplacementPath, {node_count, 6U}, "UX,UY,UZ,URX,URY,URZ",
"length,length,length,radian,radian,radian", "global-cartesian", "nodal");
projection.reaction = ReadDoubleDataset(
file.Get(), kReactionPath, {node_count, 6U}, "RF1,RF2,RF3,RM1,RM2,RM3",
"force,force,force,force*length,force*length,force*length",
"global-cartesian", "nodal");
projection.section_resultants = ReadDoubleDataset(
file.Get(), kSectionPath, {element_count, 2U, 4U}, "N,T,My,Mz",
"force,force*length,force*length,force*length", "beam-local",
"endpoint-positive-local-x-section-cut");
RequireFiniteStress(file.Get());
return projection;
}
std::vector<const WideReferenceRow*> OrderedRows(
const ReferenceTable& table, const std::vector<HdfNode>& nodes) {
if (table.rows.size() != nodes.size()) {
Fail("schema-mismatch", "The FESA/reference projected row sets differ.");
}
std::vector<const WideReferenceRow*> ordered;
ordered.reserve(nodes.size());
for (const auto& node : nodes) {
const auto found = std::find_if(
table.rows.begin(), table.rows.end(), [&](const WideReferenceRow& row) {
return AsciiLower(row.instance_name) ==
AsciiLower(node.instance_name) &&
row.source_node_label == node.source_node_label;
});
if (found == table.rows.end() ||
found->instance_name != node.instance_name) {
Fail("schema-mismatch", "A reference row identity does not match HDF5.");
}
ordered.push_back(&*found);
}
return ordered;
}
double TableScale(const ReferenceTable& table, const std::size_t value_index) {
double scale = 0.0;
for (const auto& row : table.rows) {
if (value_index >= row.values.size()) {
Fail("schema-mismatch", "A reference row has the wrong component arity.");
}
scale = (std::max)(scale, std::abs(row.values[value_index]));
}
return scale;
}
std::vector<NodeStationResultRow> NormalizeStations(
const Domain& domain, const HdfProjection& hdf,
const ReferenceTable& section_table) {
auto model_result = AnalysisModel::Create(domain);
if (!model_result.HasValue()) {
Fail("schema-mismatch",
"The approved input cannot create an analysis view.");
}
const AnalysisModel model = std::move(model_result.Value());
const std::array<double, 4> tolerances = {
kForceMomentFloor + kRelativeCoefficient * TableScale(section_table, 0U),
kForceMomentFloor + kRelativeCoefficient * TableScale(section_table, 3U),
kForceMomentFloor + kRelativeCoefficient * TableScale(section_table, 1U),
kForceMomentFloor + kRelativeCoefficient * TableScale(section_table, 2U)};
std::vector<EndpointResultRow> endpoints;
endpoints.reserve(hdf.elements.size() * 2U);
for (std::size_t element = 0U; element < hdf.elements.size(); ++element) {
for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
const auto node = static_cast<std::size_t>(
hdf.elements[element].node_internal_ids[endpoint]);
std::array<double, 4> values{};
for (std::size_t component = 0U; component < values.size(); ++component) {
values[component] =
hdf.section_resultants[(element * 2U + endpoint) * 4U + component];
}
endpoints.push_back({static_cast<EntityIndex>(element),
static_cast<int>(endpoint),
domain.Nodes()[node].source_id,
{},
values});
}
}
auto normalized = ResultRecovery::NormalizeSectionResultantsToNodeStations(
model, endpoints, tolerances);
if (!normalized.HasValue()) {
const auto& diagnostics = normalized.GetStatus().Diagnostics();
const std::string code =
diagnostics.empty() ? std::string{} : diagnostics[0U].code;
if (code == "node-station-tolerance-failure") {
Fail("tolerance-failure",
"Interior endpoint section resultants disagree.");
}
Fail("schema-mismatch",
"A node station is not eligible for legacy projection.");
}
return std::move(normalized.Value());
}
const NodeStationResultRow& FindStation(
const std::vector<NodeStationResultRow>& stations, const HdfNode& node) {
const auto found = std::find_if(
stations.begin(), stations.end(), [&](const NodeStationResultRow& row) {
return row.node.instance_name == node.instance_name &&
row.node.source_label == node.source_node_label;
});
if (found == stations.end()) {
Fail("schema-mismatch", "A projected HDF5 node station is missing.");
}
return *found;
}
CanonicalComparisonRow CanonicalRow(const HdfNode& node,
const ComparisonQuantity quantity,
std::string component, const double value,
std::string unit,
std::string coordinate_system,
std::string dataset_path) {
return {kModelId,
kStepName,
kFrameIndex,
node.instance_name,
node.source_node_label,
quantity,
std::move(component),
value,
std::move(unit),
std::move(coordinate_system),
std::move(dataset_path)};
}
void AppendNodalRows(ComparisonReport& report, const HdfProjection& hdf,
const std::vector<const WideReferenceRow*>& reference,
const ComparisonQuantity quantity,
const std::array<std::string, 6>& components,
const std::array<std::string, 6>& units,
const std::vector<double>& fesa_values,
const char* dataset_path) {
for (std::size_t node = 0U; node < hdf.nodes.size(); ++node) {
for (std::size_t component = 0U; component < components.size();
++component) {
auto fesa =
CanonicalRow(hdf.nodes[node], quantity, components[component],
fesa_values[node * 6U + component], units[component],
"global-cartesian", dataset_path);
auto abaqus =
CanonicalRow(hdf.nodes[node], quantity, components[component],
reference[node]->values[component], units[component],
"global-cartesian", dataset_path);
report.rows.push_back(
{std::move(fesa), std::move(abaqus), 0.0, 0.0, false});
}
}
}
void AppendSectionRows(ComparisonReport& report, const HdfProjection& hdf,
const std::vector<const WideReferenceRow*>& reference,
const std::vector<NodeStationResultRow>& stations) {
const std::array<std::string, 4> components = {"N", "T", "My", "Mz"};
const std::array<std::string, 4> units = {"force", "force*length",
"force*length", "force*length"};
const std::array<std::size_t, 4> reference_columns = {0U, 3U, 1U, 2U};
for (std::size_t node = 0U; node < hdf.nodes.size(); ++node) {
const auto& station = FindStation(stations, hdf.nodes[node]);
for (std::size_t component = 0U; component < components.size();
++component) {
auto fesa = CanonicalRow(
hdf.nodes[node], ComparisonQuantity::kSectionResultant,
components[component], station.section_resultant[component],
units[component], "beam-local", kSectionPath);
auto abaqus =
CanonicalRow(hdf.nodes[node], ComparisonQuantity::kSectionResultant,
components[component],
reference[node]->values[reference_columns[component]],
units[component], "beam-local", kSectionPath);
report.rows.push_back(
{std::move(fesa), std::move(abaqus), 0.0, 0.0, false});
}
}
}
double AbsoluteFloor(const ComparisonQuantity quantity, const std::string&) {
return quantity == ComparisonQuantity::kDisplacement ? kKinematicFloor
: kForceMomentFloor;
}
void EvaluateGroup(ComparisonReport& report, const ComparisonQuantity quantity,
const std::string& component) {
std::vector<std::size_t> row_indices;
for (std::size_t index = 0U; index < report.rows.size(); ++index) {
if (report.rows[index].reference.quantity == quantity &&
report.rows[index].reference.component == component) {
row_indices.push_back(index);
}
}
if (row_indices.empty()) {
Fail("schema-mismatch", "A canonical comparison component has no rows.");
}
double reference_scale = 0.0;
for (const std::size_t index : row_indices) {
reference_scale = (std::max)(reference_scale,
std::abs(report.rows[index].reference.value));
}
const double tolerance = AbsoluteFloor(quantity, component) +
kRelativeCoefficient * reference_scale;
double maximum_absolute = -1.0;
double maximum_normalized = 0.0;
std::size_t worst_row = row_indices.front();
long double squared_error = 0.0L;
for (const std::size_t index : row_indices) {
auto& row = report.rows[index];
row.absolute_error = std::abs(row.fesa.value - row.reference.value);
if (!std::isfinite(row.absolute_error)) {
Fail("schema-mismatch", "A canonical row error is nonfinite.");
}
row.tolerance = tolerance;
row.passed = row.absolute_error <= tolerance;
report.passed = report.passed && row.passed;
const double normalized = row.absolute_error / tolerance;
if (row.absolute_error > maximum_absolute) {
maximum_absolute = row.absolute_error;
worst_row = index;
}
maximum_normalized = (std::max)(maximum_normalized, normalized);
const long double error = static_cast<long double>(row.absolute_error);
squared_error += error * error;
}
const double norm_error = std::sqrt(static_cast<double>(squared_error));
const double rms_error = std::sqrt(static_cast<double>(
squared_error / static_cast<long double>(row_indices.size())));
if (!std::isfinite(norm_error) || !std::isfinite(rms_error)) {
Fail("schema-mismatch", "A component aggregate error is nonfinite.");
}
report.metrics.push_back({quantity, component, reference_scale,
maximum_absolute, maximum_normalized, rms_error,
norm_error, worst_row});
}
PhysicsEvidence MakePhysicsEvidence(const Domain& domain,
const HdfProjection& hdf) {
auto model_result = AnalysisModel::Create(domain);
if (!model_result.HasValue()) {
Fail("schema-mismatch",
"The approved input cannot create physics evidence.");
}
const AnalysisModel model = std::move(model_result.Value());
auto dofs_result = DofManager::Create(model);
if (!dofs_result.HasValue()) {
Fail("schema-mismatch", "The approved input cannot create a DOF map.");
}
const DofManager dofs = std::move(dofs_result.Value());
auto load_result = LoadAssembler::AssembleFullNodalLoad(model, dofs);
if (!load_result.HasValue()) {
Fail("schema-mismatch", "The approved input load cannot be assembled.");
}
const Vector load = std::move(load_result.Value());
if (load.Size() != hdf.reaction.size()) {
Fail("schema-mismatch", "The load and reaction spaces are inconsistent.");
}
PhysicsEvidence evidence{};
long double residual_squared = 0.0L;
for (const std::size_t free_dof : dofs.FreeDofs()) {
const long double value = static_cast<long double>(hdf.reaction[free_dof]);
residual_squared += value * value;
}
evidence.free_residual_norm =
std::sqrt(static_cast<double>(residual_squared));
for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) {
const auto& coordinates = domain.Nodes()[node].coordinates;
const std::array<double, 3> applied = {
load[node * 6U + 0U], load[node * 6U + 1U], load[node * 6U + 2U]};
const std::array<double, 3> reaction = {hdf.reaction[node * 6U + 0U],
hdf.reaction[node * 6U + 1U],
hdf.reaction[node * 6U + 2U]};
for (std::size_t component = 0U; component < 3U; ++component) {
evidence.applied_force[component] += applied[component];
evidence.reaction_force[component] += reaction[component];
evidence.applied_moment_about_origin[component] +=
load[node * 6U + 3U + component];
evidence.reaction_moment_about_origin[component] +=
hdf.reaction[node * 6U + 3U + component];
}
evidence.applied_moment_about_origin[0U] +=
coordinates[1U] * applied[2U] - coordinates[2U] * applied[1U];
evidence.applied_moment_about_origin[1U] +=
coordinates[2U] * applied[0U] - coordinates[0U] * applied[2U];
evidence.applied_moment_about_origin[2U] +=
coordinates[0U] * applied[1U] - coordinates[1U] * applied[0U];
evidence.reaction_moment_about_origin[0U] +=
coordinates[1U] * reaction[2U] - coordinates[2U] * reaction[1U];
evidence.reaction_moment_about_origin[1U] +=
coordinates[2U] * reaction[0U] - coordinates[0U] * reaction[2U];
evidence.reaction_moment_about_origin[2U] +=
coordinates[0U] * reaction[1U] - coordinates[1U] * reaction[0U];
}
evidence.endpoint_consistency_passed = true;
return evidence;
}
const char* QuantityName(const ComparisonQuantity quantity) {
switch (quantity) {
case ComparisonQuantity::kDisplacement:
return "displacement";
case ComparisonQuantity::kReaction:
return "reaction";
case ComparisonQuantity::kSectionResultant:
return "section_resultant";
}
return "unknown";
}
void WriteJsonString(std::ostream& stream, const std::string& value) {
static constexpr char kDigits[] = "0123456789abcdef";
stream.put('"');
for (const unsigned char character : value) {
switch (character) {
case '"':
stream << "\\\"";
break;
case '\\':
stream << "\\\\";
break;
case '\b':
stream << "\\b";
break;
case '\f':
stream << "\\f";
break;
case '\n':
stream << "\\n";
break;
case '\r':
stream << "\\r";
break;
case '\t':
stream << "\\t";
break;
default:
if (character < 0x20U) {
stream << "\\u00" << kDigits[character >> 4U]
<< kDigits[character & 0x0fU];
} else {
stream.put(static_cast<char>(character));
}
break;
}
}
stream.put('"');
}
void WriteCanonicalRow(std::ostream& stream,
const CanonicalComparisonRow& row) {
stream << "{\"model_id\":";
WriteJsonString(stream, row.model_id);
stream << ",\"step_name\":";
WriteJsonString(stream, row.step_name);
stream << ",\"frame_index\":" << row.frame_index << ",\"instance_name\":";
WriteJsonString(stream, row.instance_name);
stream << ",\"source_node_label\":" << row.source_node_label
<< ",\"quantity\":";
WriteJsonString(stream, QuantityName(row.quantity));
stream << ",\"component\":";
WriteJsonString(stream, row.component);
stream << ",\"value\":" << row.value << ",\"unit_dimension\":";
WriteJsonString(stream, row.unit_dimension);
stream << ",\"coordinate_system\":";
WriteJsonString(stream, row.coordinate_system);
stream << ",\"hdf5_dataset_path\":";
WriteJsonString(stream, row.hdf5_dataset_path);
stream << '}';
}
void WriteArray(std::ostream& stream, const std::array<double, 3>& values) {
stream << '[' << values[0U] << ',' << values[1U] << ',' << values[2U] << ']';
}
bool FiniteReport(const ComparisonReport& report) {
const auto finite_array = [](const std::array<double, 3>& values) {
return std::all_of(values.begin(), values.end(),
[](const double value) { return std::isfinite(value); });
};
if (!std::isfinite(report.physics_evidence.free_residual_norm) ||
!finite_array(report.physics_evidence.applied_force) ||
!finite_array(report.physics_evidence.reaction_force) ||
!finite_array(report.physics_evidence.applied_moment_about_origin) ||
!finite_array(report.physics_evidence.reaction_moment_about_origin)) {
return false;
}
for (const auto& row : report.rows) {
if (!std::isfinite(row.fesa.value) || !std::isfinite(row.reference.value) ||
!std::isfinite(row.absolute_error) || !std::isfinite(row.tolerance)) {
return false;
}
}
return std::all_of(report.metrics.begin(), report.metrics.end(),
[](const ComponentMetrics& metric) {
return std::isfinite(metric.reference_scale) &&
std::isfinite(metric.maximum_absolute_error) &&
std::isfinite(metric.maximum_normalized_error) &&
std::isfinite(metric.rms_error) &&
std::isfinite(metric.norm_error);
});
}
} // namespace
Result<ComparisonReport> ReferenceComparison::Compare(
const std::filesystem::path& results_hdf5,
const std::filesystem::path& legacy_reference_directory) {
try {
RequireExactArtifactInventory(legacy_reference_directory);
const auto input = legacy_reference_directory / kInputName;
Domain domain = ReadApprovedDomain(input);
const ReferenceTable displacement =
ReadReferenceCsv(legacy_reference_directory / kDisplacementName,
{"Frame", "Part Instance Name", "Node Label", "U-U1",
"U-U2", "U-U3", "UR-UR1", "UR-UR2", "UR-UR3"});
const ReferenceTable reaction =
ReadReferenceCsv(legacy_reference_directory / kReactionName,
{"Frame", "Part Instance Name", "Node Label", "RF-RF1",
"RF-RF2", "RF-RF3", "RM-RM1", "RM-RM2", "RM-RM3"});
const ReferenceTable section =
ReadReferenceCsv(legacy_reference_directory / kSectionName,
{"Frame", "Part Instance Name", "Node Label", "SF-SF1",
"SM-SM1", "SM-SM2", "SM-SM3"});
HdfProjection hdf = ReadHdfProjection(results_hdf5, input, domain);
const auto displacement_rows = OrderedRows(displacement, hdf.nodes);
const auto reaction_rows = OrderedRows(reaction, hdf.nodes);
const auto section_rows = OrderedRows(section, hdf.nodes);
const auto stations = NormalizeStations(domain, hdf, section);
if (stations.size() != hdf.nodes.size()) {
Fail("schema-mismatch", "The HDF5 node-station row set is incomplete.");
}
ComparisonReport report{};
report.passed = true;
AppendNodalRows(
report, hdf, displacement_rows, ComparisonQuantity::kDisplacement,
{"UX", "UY", "UZ", "URX", "URY", "URZ"},
{"length", "length", "length", "radian", "radian", "radian"},
hdf.displacement, kDisplacementPath);
AppendNodalRows(report, hdf, reaction_rows, ComparisonQuantity::kReaction,
{"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"},
{"force", "force", "force", "force*length", "force*length",
"force*length"},
hdf.reaction, kReactionPath);
AppendSectionRows(report, hdf, section_rows, stations);
for (const std::string& component :
{"UX", "UY", "UZ", "URX", "URY", "URZ"}) {
EvaluateGroup(report, ComparisonQuantity::kDisplacement, component);
}
for (const std::string& component :
{"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"}) {
EvaluateGroup(report, ComparisonQuantity::kReaction, component);
}
for (const std::string& component : {"N", "T", "My", "Mz"}) {
EvaluateGroup(report, ComparisonQuantity::kSectionResultant, component);
}
report.physics_evidence = MakePhysicsEvidence(domain, hdf);
report.stress_comparison_applicable = false;
report.stress_comparison_reason =
"Abaqus beam stress comparison is N/A; analytical/unit and HDF5 "
"schema tests provide stress evidence.";
return Result<ComparisonReport>::Success(std::move(report));
} catch (const ComparisonFailure& failure) {
return Result<ComparisonReport>::Failure(
ComparisonFailureStatus(failure.Code(), failure.what()));
} catch (const std::exception& failure) {
return Result<ComparisonReport>::Failure(
ComparisonFailureStatus("schema-mismatch", failure.what()));
}
}
Status ReferenceComparison::WriteDeterministicJson(
const ComparisonReport& report, const std::filesystem::path& output_json) {
if (output_json.empty() || output_json.filename().empty() ||
!FiniteReport(report)) {
return Status::Failure(
FailureCategory::kOutput,
{{Severity::kError,
"comparison-json-write-failure",
{},
"",
kModelId,
"The deterministic comparison report or output path is invalid."}});
}
std::ofstream stream{output_json, std::ios::binary | std::ios::trunc};
if (!stream) {
return Status::Failure(
FailureCategory::kOutput,
{{Severity::kError,
"comparison-json-write-failure",
{},
"",
kModelId,
"The deterministic comparison JSON cannot be opened."}});
}
stream.imbue(std::locale::classic());
stream << std::setprecision(std::numeric_limits<double>::max_digits10);
stream << "{\"rows\":[";
for (std::size_t index = 0U; index < report.rows.size(); ++index) {
if (index != 0U) {
stream << ',';
}
const auto& row = report.rows[index];
stream << "{\"fesa\":";
WriteCanonicalRow(stream, row.fesa);
stream << ",\"reference\":";
WriteCanonicalRow(stream, row.reference);
stream << ",\"absolute_error\":" << row.absolute_error
<< ",\"tolerance\":" << row.tolerance
<< ",\"passed\":" << (row.passed ? "true" : "false") << '}';
}
stream << "],\"metrics\":[";
for (std::size_t index = 0U; index < report.metrics.size(); ++index) {
if (index != 0U) {
stream << ',';
}
const auto& metric = report.metrics[index];
stream << "{\"quantity\":";
WriteJsonString(stream, QuantityName(metric.quantity));
stream << ",\"component\":";
WriteJsonString(stream, metric.component);
stream << ",\"reference_scale\":" << metric.reference_scale
<< ",\"maximum_absolute_error\":" << metric.maximum_absolute_error
<< ",\"maximum_normalized_error\":"
<< metric.maximum_normalized_error
<< ",\"rms_error\":" << metric.rms_error
<< ",\"norm_error\":" << metric.norm_error
<< ",\"worst_row\":" << metric.worst_row << '}';
}
stream << "],\"physics_evidence\":{\"free_residual_norm\":"
<< report.physics_evidence.free_residual_norm << ",\"applied_force\":";
WriteArray(stream, report.physics_evidence.applied_force);
stream << ",\"reaction_force\":";
WriteArray(stream, report.physics_evidence.reaction_force);
stream << ",\"applied_moment_about_origin\":";
WriteArray(stream, report.physics_evidence.applied_moment_about_origin);
stream << ",\"reaction_moment_about_origin\":";
WriteArray(stream, report.physics_evidence.reaction_moment_about_origin);
stream << ",\"endpoint_consistency_passed\":"
<< (report.physics_evidence.endpoint_consistency_passed ? "true"
: "false")
<< "},\"stress_comparison_applicable\":"
<< (report.stress_comparison_applicable ? "true" : "false")
<< ",\"stress_comparison_reason\":";
WriteJsonString(stream, report.stress_comparison_reason);
stream << ",\"passed\":" << (report.passed ? "true" : "false") << "}\n";
if (!stream) {
return Status::Failure(
FailureCategory::kOutput,
{{Severity::kError,
"comparison-json-write-failure",
{},
"",
kModelId,
"The deterministic comparison JSON write failed."}});
}
return Status::Ok();
}
} // namespace fesa::test
Reference in New Issue View Git Blame Copy Permalink