modify tolerance policy

This commit is contained in:
KOKO\Mimi
2026-08-18 02:29:26 +09:00
parent 123a1b6ccc
commit 624ea4872c
36 changed files with 1909 additions and 672 deletions
+255 -165
View File
@@ -29,7 +29,7 @@
#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"
#include "reference_tolerance_policy.h"
namespace fesa::test {
namespace {
@@ -48,9 +48,6 @@ 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:
@@ -166,6 +163,17 @@ struct ReferenceTable {
std::vector<WideReferenceRow> rows;
};
struct SectionReferenceRow {
std::string instance_name;
std::int64_t source_element_label;
std::int64_t source_node_label;
std::vector<double> values;
};
struct SectionReferenceTable {
std::vector<SectionReferenceRow> rows;
};
ReferenceTable ReadReferenceCsv(
const std::filesystem::path& path,
const std::vector<std::string>& expected_header) {
@@ -223,6 +231,70 @@ ReferenceTable ReadReferenceCsv(
return table;
}
SectionReferenceTable ReadSectionReferenceCsv(
const std::filesystem::path& path) {
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();
}
const std::vector<std::string> expected_header = {
"Frame", "Part Instance Name",
"Element Label", "Node Label",
"SF-SF1", "SM-SM1",
"SM-SM2", "SM-SM3"};
if (SplitCsvLine(line) != expected_header) {
Fail("schema-mismatch", "An approved section CSV header is not exact.");
}
SectionReferenceTable 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 section row has invalid schema or frame identity.");
}
SectionReferenceRow row{};
row.instance_name = fields[1U];
row.source_element_label = ParsePositiveLabel(fields[2U]);
row.source_node_label = ParsePositiveLabel(fields[3U]);
row.values.reserve(fields.size() - 4U);
for (std::size_t field = 4U; field < fields.size(); ++field) {
row.values.push_back(ParseFiniteDouble(fields[field]));
}
const auto duplicate = std::find_if(
table.rows.begin(), table.rows.end(),
[&](const SectionReferenceRow& existing) {
return AsciiLower(existing.instance_name) ==
AsciiLower(row.instance_name) &&
existing.source_element_label == row.source_element_label &&
existing.source_node_label == row.source_node_label;
});
if (duplicate != table.rows.end()) {
Fail("schema-mismatch",
"A reference section 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;
@@ -894,76 +966,39 @@ std::vector<const WideReferenceRow*> OrderedRows(
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]));
std::vector<const SectionReferenceRow*> OrderedSectionRows(
const SectionReferenceTable& table, const HdfProjection& hdf) {
const std::size_t endpoint_count = hdf.elements.size() * 2U;
if (table.rows.size() != endpoint_count) {
Fail("schema-mismatch", "The FESA/reference section row sets differ.");
}
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) {
std::vector<const SectionReferenceRow*> ordered;
ordered.reserve(endpoint_count);
for (const auto& element : hdf.elements) {
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];
const std::size_t node_index =
static_cast<std::size_t>(element.node_internal_ids[endpoint]);
if (node_index >= hdf.nodes.size()) {
Fail("schema-mismatch", "An HDF5 section endpoint references no node.");
}
endpoints.push_back({static_cast<EntityIndex>(element),
static_cast<int>(endpoint),
domain.Nodes()[node].source_id,
{},
values});
const auto& node = hdf.nodes[node_index];
const auto found = std::find_if(
table.rows.begin(), table.rows.end(),
[&](const SectionReferenceRow& row) {
return AsciiLower(row.instance_name) ==
AsciiLower(element.instance_name) &&
row.source_element_label == element.source_element_label &&
row.source_node_label == node.source_node_label;
});
if (found == table.rows.end() ||
found->instance_name != element.instance_name) {
Fail("schema-mismatch",
"A reference section endpoint identity does not match HDF5.");
}
ordered.push_back(&*found);
}
}
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;
return ordered;
}
CanonicalComparisonRow CanonicalRow(const HdfNode& node,
@@ -977,6 +1012,8 @@ CanonicalComparisonRow CanonicalRow(const HdfNode& node,
kFrameIndex,
node.instance_name,
node.source_node_label,
0,
-1,
quantity,
std::move(component),
value,
@@ -1003,93 +1040,122 @@ void AppendNodalRows(ComparisonReport& report, const HdfProjection& hdf,
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});
report.rows.push_back({std::move(fesa), std::move(abaqus)});
}
}
}
void AppendSectionRows(ComparisonReport& report, const HdfProjection& hdf,
const std::vector<const WideReferenceRow*>& reference,
const std::vector<NodeStationResultRow>& stations) {
void AppendSectionRows(
ComparisonReport& report, const HdfProjection& hdf,
const std::vector<const SectionReferenceRow*>& reference) {
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});
std::size_t row = 0U;
for (std::size_t element = 0U; element < hdf.elements.size(); ++element) {
for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint, ++row) {
const auto& element_row = hdf.elements[element];
const std::size_t node_index =
static_cast<std::size_t>(element_row.node_internal_ids[endpoint]);
if (node_index >= hdf.nodes.size()) {
Fail("schema-mismatch", "An HDF5 section endpoint references no node.");
}
const auto& node = hdf.nodes[node_index];
const auto canonical = [&](const double value, std::string component,
std::string unit) {
return CanonicalComparisonRow{kModelId,
kStepName,
kFrameIndex,
element_row.instance_name,
node.source_node_label,
element_row.source_element_label,
static_cast<int>(endpoint),
ComparisonQuantity::kSectionResultant,
std::move(component),
value,
std::move(unit),
"beam-local",
kSectionPath};
};
for (std::size_t component = 0U; component < components.size();
++component) {
auto fesa = canonical(
hdf.section_resultants[(element * 2U + endpoint) * 4U + component],
components[component], units[component]);
auto abaqus =
canonical(reference[row]->values[reference_columns[component]],
components[component], units[component]);
report.rows.push_back({std::move(fesa), std::move(abaqus)});
}
}
}
}
double AbsoluteFloor(const ComparisonQuantity quantity, const std::string&) {
return quantity == ComparisonQuantity::kDisplacement ? kKinematicFloor
: kForceMomentFloor;
const char* BranchName(const ReferenceToleranceBranch branch) {
switch (branch) {
case ReferenceToleranceBranch::kNearZero:
return "near-zero";
case ReferenceToleranceBranch::kRelative:
return "relative";
case ReferenceToleranceBranch::kZeroScaleExact:
return "zero-scale-exact";
}
return "unknown";
}
void EvaluateGroup(ComparisonReport& report, const ComparisonQuantity quantity,
const std::string& component) {
void EvaluateFamily(ComparisonReport& report, const ComparisonQuantity quantity,
const std::string& identity,
const std::vector<std::string>& components,
const bool blocking) {
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) {
std::find(components.begin(), components.end(),
report.rows[index].reference.component) != components.end()) {
row_indices.push_back(index);
}
}
if (row_indices.empty()) {
Fail("schema-mismatch", "A canonical comparison component has no rows.");
Fail("schema-mismatch", "A canonical comparison family has no rows.");
}
double reference_scale = 0.0;
ReferenceToleranceFamily family{};
family.identity = identity;
family.components = components;
family.blocking = blocking;
family.rows.reserve(row_indices.size());
for (const std::size_t index : row_indices) {
reference_scale = (std::max)(reference_scale,
std::abs(report.rows[index].reference.value));
family.rows.push_back(
{report.rows[index].fesa.value, 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 auto evaluation =
ReferenceTolerancePolicy::Evaluate({std::move(family)});
if (!evaluation.valid || evaluation.families.size() != 1U) {
Fail("schema-mismatch", "A canonical comparison family is nonfinite.");
}
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.");
const auto& policy = evaluation.families.front();
for (std::size_t local = 0U; local < row_indices.size(); ++local) {
auto& row = report.rows[row_indices[local]];
const auto& decision = policy.rows[local];
row.absolute_error = decision.absolute_error;
row.tolerance = decision.threshold;
row.reference_scale = policy.reference_scale;
row.near_zero_band = policy.near_zero_band;
row.relative_error = decision.relative_error;
row.relative_error_applicable = decision.relative_error_applicable;
row.tolerance_branch = BranchName(decision.branch);
row.passed = decision.passed;
}
const std::size_t worst_row = row_indices[policy.worst_row];
report.metrics.push_back({quantity, policy.identity, policy.components,
policy.reference_scale, policy.near_zero_band,
policy.near_zero_count,
policy.maximum_absolute_error, policy.relative_rms,
worst_row, policy.rms_passed, policy.passed,
policy.diagnostic_code, policy.row_count});
if (blocking && !policy.passed) {
report.passed = false;
}
report.metrics.push_back({quantity, component, reference_scale,
maximum_absolute, maximum_normalized, rms_error,
norm_error, worst_row});
}
PhysicsEvidence MakePhysicsEvidence(const Domain& domain,
@@ -1214,7 +1280,8 @@ void WriteCanonicalRow(std::ostream& stream,
stream << ",\"frame_index\":" << row.frame_index << ",\"instance_name\":";
WriteJsonString(stream, row.instance_name);
stream << ",\"source_node_label\":" << row.source_node_label
<< ",\"quantity\":";
<< ",\"source_element_label\":" << row.source_element_label
<< ",\"endpoint_index\":" << row.endpoint_index << ",\"quantity\":";
WriteJsonString(stream, QuantityName(row.quantity));
stream << ",\"component\":";
WriteJsonString(stream, row.component);
@@ -1231,6 +1298,18 @@ void WriteArray(std::ostream& stream, const std::array<double, 3>& values) {
stream << '[' << values[0U] << ',' << values[1U] << ',' << values[2U] << ']';
}
void WriteStrings(std::ostream& stream,
const std::vector<std::string>& values) {
stream << '[';
for (std::size_t index = 0U; index < values.size(); ++index) {
if (index != 0U) {
stream << ',';
}
WriteJsonString(stream, values[index]);
}
stream << ']';
}
bool FiniteReport(const ComparisonReport& report) {
const auto finite_array = [](const std::array<double, 3>& values) {
return std::all_of(values.begin(), values.end(),
@@ -1245,17 +1324,19 @@ bool FiniteReport(const ComparisonReport& report) {
}
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)) {
!std::isfinite(row.absolute_error) || !std::isfinite(row.tolerance) ||
!std::isfinite(row.reference_scale) ||
!std::isfinite(row.near_zero_band) ||
!std::isfinite(row.relative_error)) {
return false;
}
}
return std::all_of(report.metrics.begin(), report.metrics.end(),
[](const ComponentMetrics& metric) {
return std::isfinite(metric.reference_scale) &&
std::isfinite(metric.near_zero_band) &&
std::isfinite(metric.maximum_absolute_error) &&
std::isfinite(metric.maximum_normalized_error) &&
std::isfinite(metric.rms_error) &&
std::isfinite(metric.norm_error);
std::isfinite(metric.relative_rms);
});
}
@@ -1276,18 +1357,12 @@ Result<ComparisonReport> ReferenceComparison::Compare(
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"});
const SectionReferenceTable section =
ReadSectionReferenceCsv(legacy_reference_directory / kSectionName);
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.");
}
const auto section_rows = OrderedSectionRows(section, hdf);
ComparisonReport report{};
report.passed = true;
@@ -1301,19 +1376,20 @@ Result<ComparisonReport> ReferenceComparison::Compare(
{"force", "force", "force", "force*length", "force*length",
"force*length"},
hdf.reaction, kReactionPath);
AppendSectionRows(report, hdf, section_rows, stations);
AppendSectionRows(report, hdf, section_rows);
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);
}
EvaluateFamily(report, ComparisonQuantity::kDisplacement,
"displacement-translation", {"UX", "UY", "UZ"}, true);
EvaluateFamily(report, ComparisonQuantity::kDisplacement,
"displacement-rotation", {"URX", "URY", "URZ"}, true);
EvaluateFamily(report, ComparisonQuantity::kReaction, "reaction-force",
{"RF1", "RF2", "RF3"}, true);
EvaluateFamily(report, ComparisonQuantity::kReaction, "reaction-moment",
{"RM1", "RM2", "RM3"}, true);
EvaluateFamily(report, ComparisonQuantity::kSectionResultant,
"section-force", {"N"}, true);
EvaluateFamily(report, ComparisonQuantity::kSectionResultant,
"section-moment", {"T", "My", "Mz"}, true);
report.physics_evidence = MakePhysicsEvidence(domain, hdf);
report.stress_comparison_applicable = false;
report.stress_comparison_reason =
@@ -1366,8 +1442,15 @@ Status ReferenceComparison::WriteDeterministicJson(
stream << ",\"reference\":";
WriteCanonicalRow(stream, row.reference);
stream << ",\"absolute_error\":" << row.absolute_error
<< ",\"tolerance\":" << row.tolerance
<< ",\"passed\":" << (row.passed ? "true" : "false") << '}';
<< ",\"threshold\":" << row.tolerance
<< ",\"reference_scale\":" << row.reference_scale
<< ",\"near_zero_band\":" << row.near_zero_band
<< ",\"relative_error\":" << row.relative_error
<< ",\"relative_error_applicable\":"
<< (row.relative_error_applicable ? "true" : "false")
<< ",\"tolerance_branch\":";
WriteJsonString(stream, row.tolerance_branch);
stream << ",\"passed\":" << (row.passed ? "true" : "false") << '}';
}
stream << "],\"metrics\":[";
for (std::size_t index = 0U; index < report.metrics.size(); ++index) {
@@ -1377,15 +1460,22 @@ Status ReferenceComparison::WriteDeterministicJson(
const auto& metric = report.metrics[index];
stream << "{\"quantity\":";
WriteJsonString(stream, QuantityName(metric.quantity));
stream << ",\"component\":";
WriteJsonString(stream, metric.component);
stream << ",\"family_identity\":";
WriteJsonString(stream, metric.family_identity);
stream << ",\"components\":";
WriteStrings(stream, metric.components);
stream << ",\"reference_scale\":" << metric.reference_scale
<< ",\"near_zero_band\":" << metric.near_zero_band
<< ",\"near_zero_count\":" << metric.near_zero_count
<< ",\"row_count\":" << metric.row_count
<< ",\"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 << '}';
<< ",\"relative_rms\":" << metric.relative_rms
<< ",\"worst_row\":" << metric.worst_row
<< ",\"rms_passed\":" << (metric.rms_passed ? "true" : "false")
<< ",\"passed\":" << (metric.passed ? "true" : "false")
<< ",\"diagnostic_code\":";
WriteJsonString(stream, metric.diagnostic_code);
stream << '}';
}
stream << "],\"physics_evidence\":{\"free_residual_norm\":"
<< report.physics_evidence.free_residual_norm << ",\"applied_force\":";