diff --git a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md index 766e442..e0bba4c 100644 --- a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md +++ b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md @@ -1209,3 +1209,83 @@ - concerns: no remaining critical implementation or upstream-contract blocker. The HDF5 static-link and shared-runtime environment blockers were resolved and are recorded separately from the valid missing-API RED. + +## Step 24 — linear-static-cli + +- task_id: `TASK-24` +- status: `completed` +- changed_files: `include/fesa/analysis/linear_static_analysis.hpp`, + `src/fesa/analysis/linear_static_analysis.cpp`, + `include/fesa/app/fesa_application.hpp`, + `src/fesa/app/fesa_application.cpp`, `src/fesa/app/main.cpp`, + `tests/integration/analysis/linear_static_analysis_test.cpp`, + `tests/integration/app/fesa_application_test.cpp`, + `tests/reference/reference_comparison.hpp`, + `tests/reference/reference_comparison.cpp`, + `tests/reference/reference_comparison_test.cpp`, + `tests/reference/b33_reference_comparison_test.cpp`, + `src/fesa/CMakeLists.txt`, `tests/CMakeLists.txt`, + `docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md`, + `phases/linear-static-3d-euler-beam/index.json` +- requirement_ids: `FESA-REQ-LS3DEB-001`, `FESA-REQ-LS3DEB-002`, + `FESA-REQ-LS3DEB-005`, `FESA-REQ-LS3DEB-020`, + `FESA-REQ-LS3DEB-026`, `FESA-REQ-LS3DEB-028`, + `FESA-REQ-LS3DEB-030`, `FESA-REQ-LS3DEB-032`, + `FESA-REQ-LS3DEB-033`, `FESA-REQ-LS3DEB-034`, + `FESA-REQ-LS3DEB-036` through `FESA-REQ-LS3DEB-044` +- test_ids: `T24-ORCH-001`, `T24-ORCH-002`, `T24-CLI-001`, + `T24-CLI-002`, `T24-CLI-003`, `T24-REF-001`, `T24-REF-002`, + `T24-REF-003`, `T24-REF-004`, `T24-B33-001` + +| stage | exact command | exit_code | expected_or_observed_result | evidence_tail | +| --- | --- | ---: | --- | --- | +| RED-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 1 | Exact ten tests were registered while orchestration, CLI, and comparator APIs were absent | MSVC C1083 reported only missing `fesa/analysis/linear_static_analysis.hpp`, `fesa/app/fesa_application.hpp`, and `reference_comparison.hpp`; the pre-existing unit target still built | +| RED-targeted | `ctest --test-dir .harness/build -C Debug -R "LinearStaticCli|ReferenceComparisonContract|B33ReferenceComparison" --output-on-failure` | 0 | Missing executables leave no runnable Step 24 tests after the valid compile RED | CTest reported `No tests were found`; this is secondary evidence, not the RED acceptance result | +| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | Minimal production lifecycle, application boundary, test-only comparator, and all three test executables compile/link | `fesa_solver.lib`, `fesa_unit_tests.exe`, `fesa_integration_tests.exe`, and `fesa_reference_tests.exe` built without a FESA warning under `/W4 /WX` | +| GREEN-targeted | `ctest --test-dir .harness/build -C Debug -R "LinearStaticCli|ReferenceComparisonContract|B33ReferenceComparison" --output-on-failure` | 0 | Exact lifecycle/CLI/reference contract is green | Exact Step 24 suite passed 10/10: integration 5 and reference 5 | +| VERIFY-configure | `cmake -S . -B .harness/build -A x64 -DFESA_GTEST_SOURCE_DIR=C:/git/googletest "-DMKL_DIR=C:/Program Files (x86)/Intel/oneAPI/mkl/2026.1/lib/cmake/mkl" "-DTBB_DIR=C:/Program Files (x86)/Intel/oneAPI/tbb/2023.1/lib/cmake/tbb" "-DHDF5_DIR=C:/Program Files/HDF_Group/HDF5/2.1.1/cmake"` | 0 | Approved explicit-dependency MSVC x64 build tree generates | Visual Studio 18/MSVC x64, Windows SDK 10.0.26100.0, oneMKL 2026.1, oneTBB, and shared HDF5 resolved | +| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build includes the actual CLI and all test executables | `.harness/build/src/fesa/Debug/fesa.exe` built under `/W4 /WX`; selected MKL/TBB/HDF5/Intel runtime DLLs were staged beside it | +| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R "LinearStaticCli|ReferenceComparisonContract|B33ReferenceComparison" --output-on-failure` | 0 | Focused Step 24 suite remains green after the full build | Exact tests passed 10/10 | +| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` plus JSON count audit | 0 | All accumulated and exact Step 24 tests are discovered | 83 total tests, exactly 10 Step 24 names | +| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 83/83 tests passed | +| VERIFY-cli | `.harness/build/src/fesa/Debug/fesa.exe "reference/cantilever beam/cantilever beam.inp" --output .harness/build/cli-smoke-step24/results.h5` | 0 | The actual executable strips argv[0], runs the approved pipeline, and writes HDF5 | Process exit 0; `results.h5` exists and is 25,336 bytes; an empty-argument process smoke test returned exact usage exit 2 | +| VERIFY-determinism | Run exact `B33ReferenceComparison` twice and hash `comparison.json` after each run | 0 | Build-local comparison evidence is byte deterministic | Both runs passed 1/1 and produced SHA-256 `258347AEA791D981AEA9B2BCAD85DE5344D4859ECA3692DC5E7AA01A848F8E0D` | +| VERIFY-contract-scans | Exact-test count, public-backend, forbidden-scope, reference-write, whitespace, and reference diff/status scans | 0 | Approved module and artifact boundaries remain isolated | exact tests 10; public backend leaks 0; B31/DLOAD/fallback/averaging scope 0; comparator has only its JSON output writer; diff clean; reference unchanged | + +- contract_checks: `Analysis::run` invokes the exact eight hook stages and + returns immediately on failure. `LinearStaticAnalysis` keeps the immutable + Domain alive under the non-owning AnalysisModel, owns the DOF/state/matrix + handoff, and calls `factorize(Kff)` before any nodal load assembly. Solve is + substitution-only, followed by full reconstruction, recovery, and the + injected ResultsWriter. +- contract_checks: `FesaApplication::run` receives argv[0]-excluded operands, + accepts only `` or ` --output `, defaults to the current + directory `results.h5`, constructs TBB/MKL/HDF5 adapters behind their public + boundaries, emits seven deterministic diagnostic fields, and returns exact + process codes 0/2/3/4/5/6. Domain warnings remain in the HDF5 diagnostics and + output requests do not filter mandatory result datasets. +- contract_checks: the test-only comparator opens the exact legacy four-file + bundle and authoritative HDF5 read-only, requires B33/schema/source/row + identity/finite values before comparison, and produces exactly 176 stable + decisions plus 16 component metrics. It uses Abaqus-only component scales, + exact SI floors, no zero clamp/drop, and Step 22 endpoint consistency before + selecting the smaller element identity without averaging. +- generated_evidence: `.harness/build/src/fesa/Debug/fesa.exe`, + `.harness/build/reference/cantilever-beam-b33/results.h5` (25,336 bytes), + `.harness/build/reference/cantilever-beam-b33/comparison.json` (128,118 + bytes). The report has `passed=true`, 0 failed rows, 176 matched rows, 16 + metrics, stress comparison N/A, and worst normalized error + `0.026602795021994977` (`UZ`). +- physics_evidence: free residual norm `9.356339321107032e-7`; applied force + `[0,0,-1000000]`; reaction force `[0,0,999999.9999998808]`; applied moment + about origin `[0,10000000,0]`; reaction moment about origin + `[0,-9999999.999997258,0]`; endpoint consistency passed. +- reference_diff: unchanged; `git diff --exit-code -- reference/` and + `git status --short -- reference/` both report no change. Contract-fixture + mutations occur only under `.harness/build/reference/contract-fixtures/`. +- handoff: Step 25 can independently rebuild and execute all 83 tests. Step 26 + can rerun the already-built exact B33 CTest and consume the deterministic + `comparison.json` without changing source, tests, tolerances, or references; + Step 27 can consume only its `physics_evidence` after Step 26 passes. +- concerns: none; no critical implementation, environment, numerical, or + upstream-contract blocker remains. diff --git a/include/fesa/analysis/linear_static_analysis.hpp b/include/fesa/analysis/linear_static_analysis.hpp new file mode 100644 index 0000000..fcbeef8 --- /dev/null +++ b/include/fesa/analysis/linear_static_analysis.hpp @@ -0,0 +1,78 @@ +#pragma once + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/analysis/analysis_state.hpp" +#include "fesa/constraints/essential_constraints.hpp" +#include "fesa/core/status.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/math/sparse_matrix.hpp" +#include "fesa/math/vector.hpp" +#include "fesa/model/domain.hpp" + +#include +#include +#include + +namespace fesa { + +class LinearSolver; +class ParallelFor; +class ResultsWriter; + +struct AnalysisRequest { + std::filesystem::path inputPath; + std::filesystem::path outputPath; +}; + +// Fixes the public V0 lifecycle while leaving each analysis procedure to +// implement its approved stages. +class Analysis { +public: + virtual ~Analysis() = default; + Status run(const AnalysisRequest& request); + +protected: + virtual Status initialize(const AnalysisRequest& request) = 0; + virtual Status buildAnalysisModel() = 0; + virtual Status buildDofMapAndSparsePattern() = 0; + virtual Status assembleAndPartitionStiffness() = 0; + virtual Status factorize() = 0; + virtual Status assembleLoadsAndEffectiveRhs() = 0; + virtual Status substituteAndReconstruct() = 0; + virtual Status recoverAndWriteResults() = 0; +}; + +// Orchestrates the single-step B33 procedure through injected backend +// boundaries so factorization and substitution remain independently visible. +class LinearStaticAnalysis final : public Analysis { +public: + LinearStaticAnalysis(const ParallelFor& parallelFor, + LinearSolver& linearSolver, + ResultsWriter& resultsWriter); + +protected: + Status initialize(const AnalysisRequest& request) override; + Status buildAnalysisModel() override; + Status buildDofMapAndSparsePattern() override; + Status assembleAndPartitionStiffness() override; + Status factorize() override; + Status assembleLoadsAndEffectiveRhs() override; + Status substituteAndReconstruct() override; + Status recoverAndWriteResults() override; + +private: + const ParallelFor& parallelFor_; + LinearSolver& linearSolver_; + ResultsWriter& resultsWriter_; + AnalysisRequest request_; + std::unique_ptr domain_; + std::unique_ptr model_; + std::unique_ptr dofs_; + std::unique_ptr state_; + std::unique_ptr fullStiffness_; + std::unique_ptr partitionedStiffness_; + std::unique_ptr effectiveRhs_; + std::vector diagnostics_; +}; + +} // namespace fesa diff --git a/include/fesa/app/fesa_application.hpp b/include/fesa/app/fesa_application.hpp new file mode 100644 index 0000000..9111a05 --- /dev/null +++ b/include/fesa/app/fesa_application.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include +#include + +namespace fesa { + +// Owns the argv-independent command-line contract and stable process codes. +class FesaApplication { +public: + int run(const std::vector& arguments); +}; + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index c966e6a..cbaa81b 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -3,6 +3,8 @@ add_library( STATIC analysis/analysis_model.cpp analysis/analysis_state.cpp + analysis/linear_static_analysis.cpp + app/fesa_application.cpp assembly/load_assembler.cpp assembly/parallel_for.cpp assembly/sparse_assembler.cpp @@ -43,3 +45,60 @@ target_compile_options( PRIVATE $<$:/W4 /WX> ) + +add_executable(fesa_cli app/main.cpp) +set_target_properties(fesa_cli PROPERTIES OUTPUT_NAME fesa) +target_link_libraries(fesa_cli PRIVATE fesa_solver) +target_compile_options( + fesa_cli + PRIVATE + $<$:/W4 /WX> +) + +if(WIN32) + if(NOT TARGET MKL::mkl_intel_thread OR NOT TARGET MKL::mkl_core OR + NOT OMP_DLL_DIR OR NOT OMP_DLLNAME) + message(FATAL_ERROR "The oneMKL runtime files required by fesa.exe were not resolved") + endif() + if(NOT EXISTS "${OMP_DLL_DIR}/libmmd.dll") + message(FATAL_ERROR "The Intel runtime required by the selected HDF5 DLL was not resolved") + endif() + + file(GLOB _fesa_cli_mkl_dispatch_dlls "${MKL_ROOT}/bin/mkl_def.*.dll") + if(NOT _fesa_cli_mkl_dispatch_dlls) + message(FATAL_ERROR "The oneMKL default dispatch runtime was not resolved") + endif() + list(SORT _fesa_cli_mkl_dispatch_dlls COMPARE NATURAL ORDER DESCENDING) + list(GET _fesa_cli_mkl_dispatch_dlls 0 _fesa_cli_mkl_dispatch_dll) + + # Keep the installed CLI runnable without assuming oneAPI or HDF5 on PATH. + add_custom_command( + TARGET fesa_cli + POST_BUILD + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${OMP_DLL_DIR}/${OMP_DLLNAME}" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${OMP_DLL_DIR}/libmmd.dll" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${_fesa_cli_mkl_dispatch_dll}" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND_EXPAND_LISTS + ) + + unset(_fesa_cli_mkl_dispatch_dll) + unset(_fesa_cli_mkl_dispatch_dlls) +endif() diff --git a/src/fesa/analysis/linear_static_analysis.cpp b/src/fesa/analysis/linear_static_analysis.cpp new file mode 100644 index 0000000..4c877e9 --- /dev/null +++ b/src/fesa/analysis/linear_static_analysis.cpp @@ -0,0 +1,170 @@ +#include "fesa/analysis/linear_static_analysis.hpp" + +#include "fesa/assembly/load_assembler.hpp" +#include "fesa/assembly/parallel_for.hpp" +#include "fesa/assembly/sparse_assembler.hpp" +#include "fesa/io/abaqus/domain_mapper.hpp" +#include "fesa/io/abaqus/input_reader.hpp" +#include "fesa/results/result_recovery.hpp" +#include "fesa/results/results_writer.hpp" +#include "fesa/solvers/linear/linear_solver.hpp" + +#include + +namespace fesa { + +Status Analysis::run(const AnalysisRequest& request) { + Status status = initialize(request); + if (!status.isOk()) { + return status; + } + status = buildAnalysisModel(); + if (!status.isOk()) { + return status; + } + status = buildDofMapAndSparsePattern(); + if (!status.isOk()) { + return status; + } + status = assembleAndPartitionStiffness(); + if (!status.isOk()) { + return status; + } + status = factorize(); + if (!status.isOk()) { + return status; + } + status = assembleLoadsAndEffectiveRhs(); + if (!status.isOk()) { + return status; + } + status = substituteAndReconstruct(); + if (!status.isOk()) { + return status; + } + return recoverAndWriteResults(); +} + +LinearStaticAnalysis::LinearStaticAnalysis( + const ParallelFor& parallelFor, + LinearSolver& linearSolver, + ResultsWriter& resultsWriter) + : parallelFor_{parallelFor}, + linearSolver_{linearSolver}, + resultsWriter_{resultsWriter} {} + +Status LinearStaticAnalysis::initialize(const AnalysisRequest& request) { + // Clear dependent objects in reverse ownership order so a reused analysis + // never exposes a view into a Domain from an earlier run. + effectiveRhs_.reset(); + partitionedStiffness_.reset(); + fullStiffness_.reset(); + state_.reset(); + dofs_.reset(); + model_.reset(); + domain_.reset(); + diagnostics_.clear(); + request_ = request; + + const auto parsed = AbaqusInputReader{}.read(request_.inputPath); + if (!parsed.hasValue()) { + return parsed.status(); + } + auto domain = AbaqusDomainMapper{}.map(parsed.value()); + if (!domain.hasValue()) { + return domain.status(); + } + + domain_ = std::make_unique(std::move(domain.value())); + diagnostics_ = domain_->warnings(); + sortDiagnostics(diagnostics_); + return Status::ok(); +} + +Status LinearStaticAnalysis::buildAnalysisModel() { + auto model = AnalysisModel::create(*domain_); + if (!model.hasValue()) { + return model.status(); + } + model_ = std::make_unique(std::move(model.value())); + return Status::ok(); +} + +Status LinearStaticAnalysis::buildDofMapAndSparsePattern() { + auto dofs = DofManager::create(*model_); + if (!dofs.hasValue()) { + return dofs.status(); + } + dofs_ = std::make_unique(std::move(dofs.value())); + state_ = std::make_unique( + AnalysisState::create(*dofs_, {"Step-1", 0U})); + return Status::ok(); +} + +Status LinearStaticAnalysis::assembleAndPartitionStiffness() { + auto stiffness = SparseAssembler::assembleStiffness( + *model_, *dofs_, parallelFor_); + if (!stiffness.hasValue()) { + return stiffness.status(); + } + fullStiffness_ = + std::make_unique(std::move(stiffness.value())); + + auto partitioned = EssentialConstraints::partition( + *fullStiffness_, *dofs_); + if (!partitioned.hasValue()) { + return partitioned.status(); + } + partitionedStiffness_ = std::make_unique( + std::move(partitioned.value())); + return Status::ok(); +} + +Status LinearStaticAnalysis::factorize() { + // This call intentionally precedes all load assembly in Analysis::run. + return linearSolver_.factorize(partitionedStiffness_->kff); +} + +Status LinearStaticAnalysis::assembleLoadsAndEffectiveRhs() { + auto fullLoad = LoadAssembler::assembleFullNodalLoad(*model_, *dofs_); + if (!fullLoad.hasValue()) { + return fullLoad.status(); + } + state_->externalForce() = std::move(fullLoad.value()); + + auto rhs = LoadAssembler::effectiveFreeRhs( + state_->externalForce(), + partitionedStiffness_->kfc, + dofs_->prescribedValues(), + *dofs_); + if (!rhs.hasValue()) { + return rhs.status(); + } + effectiveRhs_ = std::make_unique(std::move(rhs.value())); + return Status::ok(); +} + +Status LinearStaticAnalysis::substituteAndReconstruct() { + Vector freeDisplacement{dofs_->freeDofCount()}; + const Status solveStatus = + linearSolver_.solve(*effectiveRhs_, freeDisplacement); + if (!solveStatus.isOk()) { + return solveStatus; + } + + state_->displacement() = EssentialConstraints::reconstructFull( + freeDisplacement, dofs_->prescribedValues(), *dofs_); + return Status::ok(); +} + +Status LinearStaticAnalysis::recoverAndWriteResults() { + const Status recoveryStatus = ResultRecovery::recover( + *model_, *dofs_, *fullStiffness_, *state_); + if (!recoveryStatus.isOk()) { + return recoveryStatus; + } + return resultsWriter_.write( + request_.outputPath, *domain_, *state_, diagnostics_); +} + +} // namespace fesa diff --git a/src/fesa/app/fesa_application.cpp b/src/fesa/app/fesa_application.cpp new file mode 100644 index 0000000..b7dfad5 --- /dev/null +++ b/src/fesa/app/fesa_application.cpp @@ -0,0 +1,112 @@ +#include "fesa/app/fesa_application.hpp" + +#include "fesa/analysis/linear_static_analysis.hpp" +#include "fesa/assembly/parallel_for.hpp" +#include "fesa/core/diagnostic.hpp" +#include "fesa/io/hdf5/hdf5_results_writer.hpp" +#include "fesa/solvers/linear/mkl_pardiso_solver.hpp" + +#include +#include +#include +#include + +namespace fesa { +namespace { + +constexpr int kSuccessExitCode = 0; +constexpr int kUsageExitCode = 2; +constexpr int kInputExitCode = 3; +constexpr int kModelExitCode = 4; +constexpr int kSolverExitCode = 5; +constexpr int kOutputExitCode = 6; + +bool startsWithOption(const std::string& argument) { + return !argument.empty() && argument.front() == '-'; +} + +Diagnostic usageDiagnostic() { + return { + Severity::error, + "cli-usage", + {{}, 0U}, + "", + "", + "Usage: fesa.exe [--output ]."}; +} + +const char* severityName(const Severity severity) { + return severity == Severity::warning ? "warning" : "error"; +} + +void writeDiagnostics(std::vector diagnostics) { + sortDiagnostics(diagnostics); + for (const auto& diagnostic : diagnostics) { + // Stable field labels and tab separators keep empty source fields + // explicit without depending on locale-specific formatting. + std::cerr + << "severity=" << severityName(diagnostic.severity) + << '\t' << "code=" << diagnostic.code + << '\t' << "file=" + << diagnostic.location.file.generic_u8string() + << '\t' << "line=" << diagnostic.location.line + << '\t' << "keyword=" << diagnostic.keyword + << '\t' << "entity_identity=" << diagnostic.entityIdentity + << '\t' << "message=" << diagnostic.message + << '\n'; + } +} + +int exitCodeFor(const Status& status) { + switch (status.failureCategory().value_or(FailureCategory::input)) { + case FailureCategory::input: + return kInputExitCode; + case FailureCategory::model: + return kModelExitCode; + case FailureCategory::solver: + return kSolverExitCode; + case FailureCategory::output: + return kOutputExitCode; + } + return kInputExitCode; +} + +} // namespace + +int FesaApplication::run(const std::vector& arguments) { + const bool defaultOutputForm = + arguments.size() == 1U && + !arguments[0U].empty() && + !startsWithOption(arguments[0U]); + const bool explicitOutputForm = + arguments.size() == 3U && + !arguments[0U].empty() && + !startsWithOption(arguments[0U]) && + arguments[1U] == "--output" && + !arguments[2U].empty(); + if (!defaultOutputForm && !explicitOutputForm) { + writeDiagnostics({usageDiagnostic()}); + return kUsageExitCode; + } + + AnalysisRequest request; + request.inputPath = arguments[0U]; + request.outputPath = explicitOutputForm + ? std::filesystem::path{arguments[2U]} + : std::filesystem::current_path() / "results.h5"; + + TbbParallelFor parallelFor; + MklPardisoSolver linearSolver; + Hdf5ResultsWriter resultsWriter; + LinearStaticAnalysis analysis{ + parallelFor, linearSolver, resultsWriter}; + const Status status = analysis.run(request); + if (status.isOk()) { + return kSuccessExitCode; + } + + writeDiagnostics(status.diagnostics()); + return exitCodeFor(status); +} + +} // namespace fesa diff --git a/src/fesa/app/main.cpp b/src/fesa/app/main.cpp new file mode 100644 index 0000000..cf77fc3 --- /dev/null +++ b/src/fesa/app/main.cpp @@ -0,0 +1,17 @@ +#include "fesa/app/fesa_application.hpp" + +#include +#include +#include + +int main(const int argc, char* argv[]) { + std::vector arguments; + if (argc > 1) { + arguments.reserve(static_cast(argc - 1)); + } + // The application boundary receives only operands and options, not argv[0]. + for (int index = 1; index < argc; ++index) { + arguments.emplace_back(argv[index]); + } + return fesa::FesaApplication{}.run(arguments); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7bfd53e..dfced5c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -38,6 +38,42 @@ target_link_libraries( GTest::gtest_main ) +add_executable( + fesa_integration_tests + integration/analysis/linear_static_analysis_test.cpp + integration/app/fesa_application_test.cpp +) + +target_link_libraries( + fesa_integration_tests + PRIVATE + fesa_solver + Fesa::HDF5 + GTest::gtest_main +) + +add_executable( + fesa_reference_tests + reference/reference_comparison.cpp + reference/reference_comparison_test.cpp + reference/b33_reference_comparison_test.cpp +) + +target_link_libraries( + fesa_reference_tests + PRIVATE + fesa_solver + Fesa::HDF5 + GTest::gtest_main +) + +target_compile_definitions( + fesa_reference_tests + PRIVATE + FESA_TEST_SOURCE_DIR="${PROJECT_SOURCE_DIR}" + FESA_TEST_BINARY_DIR="${PROJECT_BINARY_DIR}" +) + if(WIN32) if(NOT TARGET MKL::mkl_intel_thread OR NOT TARGET MKL::mkl_core OR NOT OMP_DLL_DIR OR NOT OMP_DLLNAME) @@ -54,36 +90,41 @@ if(WIN32) list(SORT _fesa_mkl_dispatch_dlls COMPARE NATURAL ORDER DESCENDING) list(GET _fesa_mkl_dispatch_dlls 0 _fesa_mkl_dispatch_dll) - # Stage dynamic backend dependencies before post-build GoogleTest - # discovery; the exact validation commands do not assume oneAPI in PATH. - add_custom_command( - TARGET fesa_unit_tests - POST_BUILD - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "$" - "$" - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "$" - "$" - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "$" - "$" - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "${OMP_DLL_DIR}/${OMP_DLLNAME}" - "$" - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "${OMP_DLL_DIR}/libmmd.dll" - "$" - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "${_fesa_mkl_dispatch_dll}" - "$" - # Stage runtime DLLs selected through imported targets, including the - # normalized HDF5 backend, without naming a package-specific target. - COMMAND "${CMAKE_COMMAND}" -E copy_if_different - "$" - "$" - COMMAND_EXPAND_LISTS - ) + function(fesa_stage_test_runtime target_name) + # Stage dynamic backend dependencies before post-build GoogleTest + # discovery; the exact validation commands do not assume oneAPI in PATH. + add_custom_command( + TARGET "${target_name}" + POST_BUILD + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${OMP_DLL_DIR}/${OMP_DLLNAME}" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${OMP_DLL_DIR}/libmmd.dll" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${_fesa_mkl_dispatch_dll}" + "$" + # Imported-target runtime closure includes the normalized HDF5 DLL. + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND_EXPAND_LISTS + ) + endfunction() + + fesa_stage_test_runtime(fesa_unit_tests) + fesa_stage_test_runtime(fesa_integration_tests) + fesa_stage_test_runtime(fesa_reference_tests) unset(_fesa_mkl_dispatch_dll) unset(_fesa_mkl_dispatch_dlls) @@ -95,4 +136,22 @@ gtest_discover_tests( LABELS "linear-static-3d-euler-beam;unit" ) -add_custom_target(fesa_tests DEPENDS fesa_unit_tests) +gtest_discover_tests( + fesa_integration_tests + PROPERTIES + LABELS "linear-static-3d-euler-beam;integration" +) + +gtest_discover_tests( + fesa_reference_tests + PROPERTIES + LABELS "linear-static-3d-euler-beam;reference" +) + +add_custom_target( + fesa_tests + DEPENDS + fesa_unit_tests + fesa_integration_tests + fesa_reference_tests +) diff --git a/tests/integration/analysis/linear_static_analysis_test.cpp b/tests/integration/analysis/linear_static_analysis_test.cpp new file mode 100644 index 0000000..867c053 --- /dev/null +++ b/tests/integration/analysis/linear_static_analysis_test.cpp @@ -0,0 +1,300 @@ +#include "fesa/analysis/linear_static_analysis.hpp" + +#include "fesa/assembly/parallel_for.hpp" +#include "fesa/results/results_writer.hpp" +#include "fesa/solvers/linear/mkl_pardiso_solver.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +class TempDirectory { +public: + explicit TempDirectory(const std::string& label) { + static std::atomic sequence{0U}; + const auto tick = std::chrono::steady_clock::now() + .time_since_epoch() + .count(); + path_ = std::filesystem::temp_directory_path() / + ("fesa-step24-analysis-" + label + "-" + + std::to_string(tick) + "-" + + std::to_string(sequence.fetch_add(1U))); + std::error_code error; + if (!std::filesystem::create_directory(path_, error) || error) { + throw std::runtime_error{"Unable to create the Step 24 analysis fixture."}; + } + } + + TempDirectory(const TempDirectory&) = delete; + TempDirectory& operator=(const TempDirectory&) = delete; + + ~TempDirectory() { + std::error_code ignored; + std::filesystem::remove_all(path_, ignored); + } + + const std::filesystem::path& path() const noexcept { return path_; } + +private: + std::filesystem::path path_; +}; + +void writeText(const std::filesystem::path& path, const std::string& text) { + std::ofstream stream{path, std::ios::binary | std::ios::trunc}; + stream.write(text.data(), static_cast(text.size())); + if (!stream) { + throw std::runtime_error{"Unable to write the Step 24 analysis input."}; + } +} + +std::string axialDeck(const double rootUx, const double tipForce) { + return R"inp(*Part, name=BeamPart +*Node +1, 0., 0., 0. +2, 2., 0., 0. +*Element, type=B33 +1, 1, 2 +*Elset, elset=BeamSet +1 +*Beam General Section, elset=BeamSet, material=Steel, section=GENERAL +2., 0.5, 0., 0.75, 0.25 +0., 1., 0. +*End Part +*Assembly, name=Assembly +*Instance, name=Beam-1, part=BeamPart +*End Instance +*Nset, nset=Root, instance=Beam-1 +1 +*Nset, nset=Tip, instance=Beam-1 +2 +*End Assembly +*Material, name=Steel +*Elastic +100., 0.25 +*Boundary +Root, 1, 1, )inp" + std::to_string(rootUx) + R"inp( +Root, 2, 6 +Tip, 2, 6 +*Step, name=Load, nlgeom=NO +*Static +0.1, 1., 0.01, 1. +*Cload +Tip, 1, )inp" + std::to_string(tipForce) + R"inp( +*End Step +)inp"; +} + +// The pure Template Method spy makes the eight public lifecycle hooks observable +// without coupling the ordering assertion to any solver backend. +class SpyAnalysis final : public fesa::Analysis { +public: + const std::vector& events() const noexcept { return events_; } + +protected: + fesa::Status initialize(const fesa::AnalysisRequest&) override { + return record("initialize"); + } + fesa::Status buildAnalysisModel() override { + return record("build-analysis-model"); + } + fesa::Status buildDofMapAndSparsePattern() override { + return record("build-dof-map-and-sparse-pattern"); + } + fesa::Status assembleAndPartitionStiffness() override { + return record("assemble-and-partition-stiffness"); + } + fesa::Status factorize() override { return record("factorize"); } + fesa::Status assembleLoadsAndEffectiveRhs() override { + return record("assemble-loads-and-effective-rhs"); + } + fesa::Status substituteAndReconstruct() override { + return record("substitute-and-reconstruct"); + } + fesa::Status recoverAndWriteResults() override { + return record("recover-and-write-results"); + } + +private: + fesa::Status record(const char* event) { + events_.emplace_back(event); + return fesa::Status::ok(); + } + + std::vector events_; +}; + +// The solver spy records only the adapter-boundary operations. In particular, +// solve() cannot conceal a second factorization call. +class SpyLinearSolver final : public fesa::LinearSolver { +public: + explicit SpyLinearSolver(std::vector& events) + : events_{events} {} + + fesa::Status factorize(const fesa::SparseMatrix&) override { + ++factorizeCalls_; + events_.emplace_back("solver-factorize"); + return fesa::Status::ok(); + } + + fesa::Status solve( + const fesa::Vector& rhs, fesa::Vector& solution) const override { + ++solveCalls_; + events_.emplace_back("solver-solve"); + for (std::size_t index = 0U; + index < rhs.size() && index < solution.size(); ++index) { + solution[index] = 0.0; + } + return fesa::Status::ok(); + } + + int factorizeCalls() const noexcept { return factorizeCalls_; } + int solveCalls() const noexcept { return solveCalls_; } + +private: + std::vector& events_; + int factorizeCalls_{0}; + mutable int solveCalls_{0}; +}; + +class SpyResultsWriter final : public fesa::ResultsWriter { +public: + explicit SpyResultsWriter(std::vector& events) + : events_{events} {} + + fesa::Status write( + const std::filesystem::path&, + const fesa::Domain&, + const fesa::AnalysisState&, + const std::vector&) override { + ++writeCalls_; + events_.emplace_back("writer-write"); + return fesa::Status::ok(); + } + + int writeCalls() const noexcept { return writeCalls_; } + +private: + std::vector& events_; + int writeCalls_{0}; +}; + +class CapturingResultsWriter final : public fesa::ResultsWriter { +public: + fesa::Status write( + const std::filesystem::path& outputPath, + const fesa::Domain& domain, + const fesa::AnalysisState& state, + const std::vector& diagnostics) override { + outputPath_ = outputPath; + nodeCount_ = domain.nodes().size(); + state_ = std::make_unique(state); + diagnostics_ = diagnostics; + return fesa::Status::ok(); + } + + const fesa::AnalysisState& state() const { + if (!state_) { + throw std::logic_error{"No AnalysisState was captured."}; + } + return *state_; + } + + const std::filesystem::path& outputPath() const noexcept { + return outputPath_; + } + std::size_t nodeCount() const noexcept { return nodeCount_; } + const std::vector& diagnostics() const noexcept { + return diagnostics_; + } + +private: + std::filesystem::path outputPath_; + std::size_t nodeCount_{0U}; + std::unique_ptr state_; + std::vector diagnostics_; +}; + +} // namespace + +TEST(LinearStaticCli, FactorizesBeforeLoadAndSolvesWithoutRefactorization) { + SpyAnalysis lifecycle; + const fesa::AnalysisRequest emptyRequest{}; + ASSERT_TRUE(lifecycle.run(emptyRequest).isOk()); + EXPECT_EQ( + lifecycle.events(), + (std::vector{ + "initialize", + "build-analysis-model", + "build-dof-map-and-sparse-pattern", + "assemble-and-partition-stiffness", + "factorize", + "assemble-loads-and-effective-rhs", + "substitute-and-reconstruct", + "recover-and-write-results"})); + + TempDirectory directory{"order"}; + const auto input = directory.path() / "order.inp"; + const auto output = directory.path() / "results.h5"; + writeText(input, axialDeck(0.0, 0.0)); + + std::vector adapterEvents; + fesa::SerialParallelFor serial; + SpyLinearSolver solver{adapterEvents}; + SpyResultsWriter writer{adapterEvents}; + fesa::LinearStaticAnalysis analysis{serial, solver, writer}; + + ASSERT_TRUE(analysis.run({input, output}).isOk()); + EXPECT_EQ(solver.factorizeCalls(), 1); + EXPECT_EQ(solver.solveCalls(), 1); + EXPECT_EQ(writer.writeCalls(), 1); + EXPECT_EQ( + adapterEvents, + (std::vector{ + "solver-factorize", "solver-solve", "writer-write"})); +} + +TEST(LinearStaticCli, RealPipelineHandlesAnalyticalAndNonzeroPrescription) { + TempDirectory directory{"analytical"}; + const auto input = directory.path() / "prescribed-axial.inp"; + const auto output = directory.path() / "captured-results.h5"; + writeText(input, axialDeck(0.1, 10.0)); + + fesa::SerialParallelFor serial; + fesa::MklPardisoSolver solver; + CapturingResultsWriter writer; + fesa::LinearStaticAnalysis analysis{serial, solver, writer}; + + const auto status = analysis.run({input, output}); + ASSERT_TRUE(status.isOk()); + EXPECT_EQ(writer.outputPath(), output); + EXPECT_EQ(writer.nodeCount(), 2U); + EXPECT_TRUE(writer.diagnostics().empty()); + + const auto& state = writer.state(); + ASSERT_EQ(state.displacement().size(), 12U); + EXPECT_EQ(state.identity().stepName, "Step-1"); + EXPECT_EQ(state.identity().frameIndex, 0U); + + // EA/L = 100 for this fixture, so u_tip = 0.1 + 10/100 = 0.2. + EXPECT_NEAR(state.displacement()[0U], 0.1, 1.0e-12); + EXPECT_NEAR(state.displacement()[6U], 0.2, 2.0e-10); + EXPECT_NEAR(state.externalForce()[6U], 10.0, 1.0e-12); + EXPECT_NEAR(state.internalForce()[0U], -10.0, 1.0e-9); + EXPECT_NEAR(state.internalForce()[6U], 10.0, 1.0e-9); + EXPECT_NEAR(state.reaction()[0U], -10.0, 1.0e-9); + EXPECT_NEAR(state.residual()[6U], 0.0, 1.0e-9); + EXPECT_EQ(state.endpointResults().size(), 2U); + EXPECT_EQ(state.gaussResults().size(), 2U); + EXPECT_EQ(state.stressResults().size(), 2U); +} diff --git a/tests/integration/app/fesa_application_test.cpp b/tests/integration/app/fesa_application_test.cpp new file mode 100644 index 0000000..57ce559 --- /dev/null +++ b/tests/integration/app/fesa_application_test.cpp @@ -0,0 +1,400 @@ +#include "fesa/app/fesa_application.hpp" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr const char* kStepRoot = "/steps/Step-1/frames/0"; + +class TempDirectory { +public: + explicit TempDirectory(const std::string& label) { + static std::atomic sequence{0U}; + const auto tick = std::chrono::steady_clock::now() + .time_since_epoch() + .count(); + path_ = std::filesystem::temp_directory_path() / + ("fesa-step24-app-" + label + "-" + + std::to_string(tick) + "-" + + std::to_string(sequence.fetch_add(1U))); + std::error_code error; + if (!std::filesystem::create_directory(path_, error) || error) { + throw std::runtime_error{"Unable to create the Step 24 app fixture."}; + } + } + + TempDirectory(const TempDirectory&) = delete; + TempDirectory& operator=(const TempDirectory&) = delete; + + ~TempDirectory() { + std::error_code ignored; + std::filesystem::remove_all(path_, ignored); + } + + const std::filesystem::path& path() const noexcept { return path_; } + +private: + std::filesystem::path path_; +}; + +class CurrentDirectoryGuard { +public: + explicit CurrentDirectoryGuard(const std::filesystem::path& replacement) + : original_{std::filesystem::current_path()} { + std::filesystem::current_path(replacement); + } + + CurrentDirectoryGuard(const CurrentDirectoryGuard&) = delete; + CurrentDirectoryGuard& operator=(const CurrentDirectoryGuard&) = delete; + + ~CurrentDirectoryGuard() { + std::error_code ignored; + std::filesystem::current_path(original_, ignored); + } + +private: + std::filesystem::path original_; +}; + +class Hdf5Handle { +public: + using Closer = herr_t (*)(hid_t); + + Hdf5Handle(const hid_t value, Closer closer) + : value_{value}, closer_{closer} {} + Hdf5Handle(const Hdf5Handle&) = delete; + Hdf5Handle& operator=(const Hdf5Handle&) = delete; + Hdf5Handle(Hdf5Handle&& other) noexcept + : value_{other.value_}, closer_{other.closer_} { + other.value_ = -1; + other.closer_ = nullptr; + } + ~Hdf5Handle() { + if (value_ >= 0 && closer_ != nullptr) { + (void)closer_(value_); + } + } + + hid_t get() const noexcept { return value_; } + +private: + hid_t value_{-1}; + Closer closer_{nullptr}; +}; + +void writeText(const std::filesystem::path& path, const std::string& text) { + std::ofstream stream{path, std::ios::binary | std::ios::trunc}; + stream.write(text.data(), static_cast(text.size())); + if (!stream) { + throw std::runtime_error{"Unable to write the Step 24 app input."}; + } +} + +std::string axialDeck( + const bool constrained, + const bool zeroLength, + const bool outputRequests) { + const std::string secondNode = zeroLength + ? "2, 0., 0., 0.\n" + : "2, 2., 0., 0.\n"; + const std::string boundaries = constrained + ? "*Boundary\nRoot, 1, 6\nTip, 2, 6\n" + : ""; + const std::string outputs = outputRequests + ? R"inp(*Output, field +*Node Output +U, RF +*Element Output, directions=YES +S, SF +*Output, history +*Contact Output +)inp" + : ""; + + return std::string{R"inp(*Part, name=BeamPart +*Node +1, 0., 0., 0. +)inp"} + secondNode + R"inp(*Element, type=B33 +1, 1, 2 +*Elset, elset=BeamSet +1 +*Beam General Section, elset=BeamSet, material=Steel, section=GENERAL +2., 0.5, 0., 0.75, 0.25 +0., 1., 0. +*End Part +*Assembly, name=Assembly +*Instance, name=Beam-1, part=BeamPart +*End Instance +*Nset, nset=Root, instance=Beam-1 +1 +*Nset, nset=Tip, instance=Beam-1 +2 +*End Assembly +*Material, name=Steel +*Elastic +100., 0.25 +)inp" + boundaries + R"inp(*Step, name=Load, nlgeom=NO +*Static +0.1, 1., 0.01, 1. +*Cload +Tip, 1, 10. +)inp" + outputs + R"inp(*End Step +)inp"; +} + +Hdf5Handle openFile(const std::filesystem::path& path) { + const hid_t file = H5Fopen( + path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + if (file < 0) { + throw std::runtime_error{"Unable to open the CLI HDF5 artifact."}; + } + return Hdf5Handle{file, H5Fclose}; +} + +Hdf5Handle openDataset(const hid_t file, const std::string& path) { + const hid_t dataset = H5Dopen2(file, path.c_str(), H5P_DEFAULT); + if (dataset < 0) { + throw std::runtime_error{"Unable to open mandatory dataset: " + path}; + } + return Hdf5Handle{dataset, H5Dclose}; +} + +std::vector datasetDimensions( + const hid_t file, const std::string& path) { + const auto dataset = openDataset(file, path); + Hdf5Handle space{H5Dget_space(dataset.get()), H5Sclose}; + const int rank = H5Sget_simple_extent_ndims(space.get()); + if (space.get() < 0 || rank < 0) { + throw std::runtime_error{"Unable to inspect mandatory dataset dimensions."}; + } + std::vector dimensions(static_cast(rank)); + if (rank > 0 && + H5Sget_simple_extent_dims(space.get(), dimensions.data(), nullptr) < 0) { + throw std::runtime_error{"Unable to read mandatory dataset dimensions."}; + } + return dimensions; +} + +std::vector readDoubleDataset( + const hid_t file, const std::string& path) { + const auto dimensions = datasetDimensions(file, path); + std::size_t count = 1U; + for (const hsize_t dimension : dimensions) { + count *= static_cast(dimension); + } + const auto dataset = openDataset(file, path); + std::vector values(count); + if (!values.empty() && + H5Dread(dataset.get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()) < 0) { + throw std::runtime_error{"Unable to read mandatory numeric results."}; + } + return values; +} + +std::string readStringAttribute(const hid_t object, const char* name) { + Hdf5Handle attribute{H5Aopen(object, name, H5P_DEFAULT), H5Aclose}; + Hdf5Handle type{H5Aget_type(attribute.get()), H5Tclose}; + if (attribute.get() < 0 || type.get() < 0 || + H5Tget_class(type.get()) != H5T_STRING || + H5Tis_variable_str(type.get()) <= 0) { + throw std::runtime_error{"Expected a variable-length string attribute."}; + } + char* raw = nullptr; + if (H5Aread(attribute.get(), type.get(), &raw) < 0 || raw == nullptr) { + throw std::runtime_error{"Unable to read the HDF5 identity attribute."}; + } + const std::string value{raw}; + (void)H5free_memory(raw); + return value; +} + +void expectMandatoryInventory(const hid_t file) { + for (const char* path : { + "/metadata", + "/model/nodes", + "/model/elements", + "/steps/Step-1/frames/0/nodal/displacement", + "/steps/Step-1/frames/0/nodal/reaction", + "/steps/Step-1/frames/0/element/end_force_local", + "/steps/Step-1/frames/0/element/section_resultant", + "/steps/Step-1/frames/0/element/generalized_strain", + "/steps/Step-1/frames/0/element/generalized_resultant", + "/steps/Step-1/frames/0/element/stress_s11", + "/diagnostics"}) { + EXPECT_GT(H5Lexists(file, path, H5P_DEFAULT), 0) << path; + } +} + +void expectFesaHdf5Identity( + const std::filesystem::path& output, + const std::filesystem::path& input) { + ASSERT_TRUE(std::filesystem::exists(output)); + ASSERT_GT(H5Fis_hdf5(output.string().c_str()), 0); + const auto file = openFile(output); + expectMandatoryInventory(file.get()); + Hdf5Handle metadata{ + H5Gopen2(file.get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(metadata.get(), 0); + EXPECT_EQ( + readStringAttribute(metadata.get(), "feature_id"), + "linear-static-3d-euler-beam"); + const std::string normalizedInput = + std::filesystem::absolute(input).lexically_normal().generic_u8string(); + EXPECT_EQ( + readStringAttribute(metadata.get(), "source_input_identity").find( + "path=" + normalizedInput + ";content_identity="), + 0U); +} + +struct AppRun { + int exitCode; + std::string standardError; +}; + +AppRun runApplication(const std::vector& arguments) { + testing::internal::CaptureStderr(); + try { + const int exitCode = fesa::FesaApplication{}.run(arguments); + return {exitCode, testing::internal::GetCapturedStderr()}; + } catch (...) { + (void)testing::internal::GetCapturedStderr(); + throw; + } +} + +void expectDiagnosticFieldOrder(const std::string& text) { + ASSERT_FALSE(text.empty()); + std::size_t cursor = 0U; + for (const char* field : { + "severity", "code", "file", "line", "keyword", + "entity_identity", "message"}) { + const auto position = text.find(field, cursor); + ASSERT_NE(position, std::string::npos) + << "Missing or out-of-order diagnostic field: " << field + << "\nstderr:\n" << text; + cursor = position + std::string{field}.size(); + } +} + +std::vector explicitOutputArguments( + const std::filesystem::path& input, + const std::filesystem::path& output) { + // FesaApplication receives argv[0]-excluded operands and options. + return {input.string(), "--output", output.string()}; +} + +} // namespace + +TEST(LinearStaticCli, DefaultAndExplicitOutputProduceFesaHdf5) { + TempDirectory directory{"paths"}; + const auto input = directory.path() / "model.inp"; + writeText(input, axialDeck(true, false, false)); + + const auto defaultOutput = directory.path() / "results.h5"; + { + CurrentDirectoryGuard currentDirectory{directory.path()}; + const auto result = runApplication({input.string()}); + ASSERT_EQ(result.exitCode, 0) << result.standardError; + } + expectFesaHdf5Identity(defaultOutput, input); + + const auto explicitOutput = directory.path() / "named-output.h5"; + const auto result = runApplication( + explicitOutputArguments(input, explicitOutput)); + ASSERT_EQ(result.exitCode, 0) << result.standardError; + expectFesaHdf5Identity(explicitOutput, input); +} + +TEST(LinearStaticCli, ReturnsEveryExactExitCodeAndOrderedDiagnostic) { + TempDirectory directory{"exit-codes"}; + const auto validInput = directory.path() / "valid.inp"; + const auto modelInput = directory.path() / "invalid-model.inp"; + const auto solverInput = directory.path() / "singular.inp"; + writeText(validInput, axialDeck(true, false, false)); + writeText(modelInput, axialDeck(true, true, false)); + writeText(solverInput, axialDeck(false, false, false)); + + const auto success = runApplication(explicitOutputArguments( + validInput, directory.path() / "success.h5")); + const auto usage = runApplication({}); + const auto missingInput = directory.path() / "missing.inp"; + const auto input = runApplication({missingInput.string()}); + const auto repeatedInput = runApplication({missingInput.string()}); + const auto model = runApplication(explicitOutputArguments( + modelInput, directory.path() / "model-failure.h5")); + const auto solver = runApplication(explicitOutputArguments( + solverInput, directory.path() / "solver-failure.h5")); + const auto output = runApplication(explicitOutputArguments( + validInput, + directory.path() / "nonexistent-parent" / "results.h5")); + + EXPECT_EQ(success.exitCode, 0) << success.standardError; + EXPECT_EQ(usage.exitCode, 2); + EXPECT_EQ(input.exitCode, 3); + EXPECT_EQ(model.exitCode, 4); + EXPECT_EQ(solver.exitCode, 5); + EXPECT_EQ(output.exitCode, 6); + + expectDiagnosticFieldOrder(usage.standardError); + expectDiagnosticFieldOrder(input.standardError); + expectDiagnosticFieldOrder(model.standardError); + expectDiagnosticFieldOrder(solver.standardError); + expectDiagnosticFieldOrder(output.standardError); + EXPECT_EQ(repeatedInput.exitCode, 3); + EXPECT_EQ(repeatedInput.standardError, input.standardError); +} + +TEST(LinearStaticCli, OutputRequestsDoNotFilterMandatoryResults) { + TempDirectory directory{"output-requests"}; + const auto plainInput = directory.path() / "plain.inp"; + const auto requestedInput = directory.path() / "requested.inp"; + const auto plainOutput = directory.path() / "plain.h5"; + const auto requestedOutput = directory.path() / "requested.h5"; + writeText(plainInput, axialDeck(true, false, false)); + writeText(requestedInput, axialDeck(true, false, true)); + + const auto plain = runApplication( + explicitOutputArguments(plainInput, plainOutput)); + const auto requested = runApplication( + explicitOutputArguments(requestedInput, requestedOutput)); + ASSERT_EQ(plain.exitCode, 0) << plain.standardError; + ASSERT_EQ(requested.exitCode, 0) << requested.standardError; + + const auto plainFile = openFile(plainOutput); + const auto requestedFile = openFile(requestedOutput); + expectMandatoryInventory(plainFile.get()); + expectMandatoryInventory(requestedFile.get()); + + for (const char* suffix : { + "/nodal/displacement", + "/nodal/reaction", + "/element/end_force_local", + "/element/section_resultant", + "/element/generalized_strain", + "/element/generalized_resultant"}) { + const std::string path = std::string{kStepRoot} + suffix; + EXPECT_EQ( + readDoubleDataset(requestedFile.get(), path), + readDoubleDataset(plainFile.get(), path)) + << path; + } + EXPECT_EQ(datasetDimensions(plainFile.get(), "/diagnostics"), + std::vector({0U})); + const auto requestedDiagnostics = + datasetDimensions(requestedFile.get(), "/diagnostics"); + ASSERT_EQ(requestedDiagnostics.size(), 1U); + EXPECT_GT(requestedDiagnostics[0U], 0U); +} diff --git a/tests/reference/b33_reference_comparison_test.cpp b/tests/reference/b33_reference_comparison_test.cpp new file mode 100644 index 0000000..e34fdc6 --- /dev/null +++ b/tests/reference/b33_reference_comparison_test.cpp @@ -0,0 +1,276 @@ +#include "reference_comparison.hpp" + +#include "fesa/app/fesa_application.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef FESA_TEST_SOURCE_DIR +#error FESA_TEST_SOURCE_DIR must identify the repository root. +#endif + +#ifndef FESA_TEST_BINARY_DIR +#error FESA_TEST_BINARY_DIR must identify the CMake binary root. +#endif + +namespace { + +constexpr const char* kStressPath = + "/steps/Step-1/frames/0/element/stress_s11"; +constexpr std::size_t kExpectedRowCount = 176U; +constexpr std::size_t kExpectedMetricCount = 16U; + +class Hdf5Handle { +public: + using Closer = herr_t (*)(hid_t); + + Hdf5Handle(const hid_t value, Closer closer) + : value_{value}, closer_{closer} {} + Hdf5Handle(const Hdf5Handle&) = delete; + Hdf5Handle& operator=(const Hdf5Handle&) = delete; + ~Hdf5Handle() { + if (value_ >= 0 && closer_ != nullptr) { + (void)closer_(value_); + } + } + + hid_t get() const noexcept { return value_; } + +private: + hid_t value_; + Closer closer_; +}; + +struct ReferenceSnapshotEntry { + std::filesystem::path relativePath; + bool isDirectory; + std::string bytes; + std::filesystem::file_time_type lastWriteTime; +}; + +std::string readBytes(const std::filesystem::path& path) { + std::ifstream stream{path, std::ios::binary}; + if (!stream) { + throw std::runtime_error{"Unable to read reference evidence: " + + path.string()}; + } + return {std::istreambuf_iterator{stream}, + std::istreambuf_iterator{}}; +} + +std::vector snapshotTree( + const std::filesystem::path& root) { + std::vector entries; + for (const auto& entry : std::filesystem::recursive_directory_iterator{root}) { + const bool isDirectory = entry.is_directory(); + if (!isDirectory && !entry.is_regular_file()) { + throw std::runtime_error{"Unexpected reference-tree entry type."}; + } + entries.push_back({ + std::filesystem::relative(entry.path(), root), + isDirectory, + isDirectory ? std::string{} : readBytes(entry.path()), + entry.last_write_time()}); + } + std::sort( + entries.begin(), + entries.end(), + [](const ReferenceSnapshotEntry& left, + const ReferenceSnapshotEntry& right) { + return left.relativePath.generic_string() < + right.relativePath.generic_string(); + }); + return entries; +} + +void expectTreeUnchanged( + const std::vector& before, + const std::vector& after) { + ASSERT_EQ(after.size(), before.size()); + for (std::size_t index = 0U; index < before.size(); ++index) { + EXPECT_EQ(after[index].relativePath, before[index].relativePath); + EXPECT_EQ(after[index].isDirectory, before[index].isDirectory); + EXPECT_EQ(after[index].bytes, before[index].bytes) + << before[index].relativePath.string(); + EXPECT_EQ(after[index].lastWriteTime, before[index].lastWriteTime) + << before[index].relativePath.string(); + } +} + +double norm(const std::array& value) { + return std::sqrt( + value[0U] * value[0U] + value[1U] * value[1U] + + value[2U] * value[2U]); +} + +std::array sum( + const std::array& left, + const std::array& right) { + return { + left[0U] + right[0U], + left[1U] + right[1U], + left[2U] + right[2U]}; +} + +std::size_t stressRowCount(const std::filesystem::path& results) { + const hid_t fileId = + H5Fopen(results.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + if (fileId < 0) { + throw std::runtime_error{"Unable to open authoritative HDF5 output."}; + } + const Hdf5Handle file{fileId, H5Fclose}; + if (H5Lexists(file.get(), kStressPath, H5P_DEFAULT) <= 0) { + throw std::runtime_error{"Mandatory stress_s11 dataset is missing."}; + } + const hid_t datasetId = H5Dopen2(file.get(), kStressPath, H5P_DEFAULT); + if (datasetId < 0) { + throw std::runtime_error{"Unable to open mandatory stress_s11 dataset."}; + } + const Hdf5Handle dataset{datasetId, H5Dclose}; + const hid_t spaceId = H5Dget_space(dataset.get()); + if (spaceId < 0) { + throw std::runtime_error{"Unable to inspect stress_s11 dataspace."}; + } + const Hdf5Handle space{spaceId, H5Sclose}; + if (H5Sget_simple_extent_ndims(space.get()) != 1) { + throw std::runtime_error{"stress_s11 must be a flat row dataset."}; + } + hsize_t count = 0U; + if (H5Sget_simple_extent_dims(space.get(), &count, nullptr) < 0) { + throw std::runtime_error{"Unable to read stress_s11 extent."}; + } + return static_cast(count); +} + +} // namespace + +TEST(B33ReferenceComparison, + GeneratesAuthoritativeHdf5AndComparisonEvidence) { + const std::filesystem::path sourceRoot{FESA_TEST_SOURCE_DIR}; + const std::filesystem::path binaryRoot{FESA_TEST_BINARY_DIR}; + const auto referenceDirectory = + sourceRoot / "reference" / "cantilever beam"; + const auto input = referenceDirectory / "cantilever beam.inp"; + const auto outputDirectory = + binaryRoot / "reference" / "cantilever-beam-b33"; + const auto results = outputDirectory / "results.h5"; + const auto comparison = outputDirectory / "comparison.json"; + + // Only the exact build-local evidence directory is reset; the approved + // reference tree is snapshotted and subsequently opened read-only. + std::error_code error; + std::filesystem::remove_all(outputDirectory, error); + error.clear(); + ASSERT_TRUE(std::filesystem::create_directories(outputDirectory, error)); + ASSERT_FALSE(error); + const auto referenceBefore = snapshotTree(referenceDirectory); + ASSERT_EQ(referenceBefore.size(), 4U); + EXPECT_EQ( + referenceBefore[0U].relativePath, + std::filesystem::path{"cantilever beam displacements.csv"}); + EXPECT_EQ( + referenceBefore[1U].relativePath, + std::filesystem::path{"cantilever beam elemental forces.csv"}); + EXPECT_EQ( + referenceBefore[2U].relativePath, + std::filesystem::path{"cantilever beam reactions.csv"}); + EXPECT_EQ( + referenceBefore[3U].relativePath, + std::filesystem::path{"cantilever beam.inp"}); + + fesa::FesaApplication application; + // FesaApplication receives application operands/options; main strips argv[0]. + ASSERT_EQ( + application.run( + {input.string(), "--output", results.string()}), + 0); + ASSERT_TRUE(std::filesystem::is_regular_file(results)); + ASSERT_GT(H5Fis_hdf5(results.string().c_str()), 0); + EXPECT_GT(stressRowCount(results), 0U); + + auto comparisonResult = fesa::test::ReferenceComparison::compare( + results, referenceDirectory); + ASSERT_TRUE(comparisonResult.hasValue()); + const auto& report = comparisonResult.value(); + ASSERT_TRUE(report.passed); + ASSERT_EQ(report.rows.size(), kExpectedRowCount); + ASSERT_EQ(report.metrics.size(), kExpectedMetricCount); + EXPECT_TRUE(std::all_of( + report.rows.begin(), + report.rows.end(), + [](const fesa::test::RowDecision& row) { + return row.passed && std::isfinite(row.absoluteError) && + std::isfinite(row.tolerance) && row.tolerance > 0.0 && + row.fesa.modelId == "cantilever-beam-b33" && + row.reference.modelId == "cantilever-beam-b33" && + row.fesa.stepName == "Step-1" && + row.reference.stepName == "Step-1" && + row.fesa.frameIndex == 0U && + row.reference.frameIndex == 0U && + row.fesa.instanceName == "PART-1_1-1" && + row.reference.instanceName == "PART-1_1-1" && + !row.fesa.hdf5DatasetPath.empty(); + })); + EXPECT_TRUE(std::all_of( + report.metrics.begin(), + report.metrics.end(), + [](const fesa::test::ComponentMetrics& metric) { + return std::isfinite(metric.referenceScale) && + std::isfinite(metric.maximumAbsoluteError) && + std::isfinite(metric.maximumNormalizedError) && + std::isfinite(metric.rmsError) && + std::isfinite(metric.normError) && + metric.maximumNormalizedError <= 1.0; + })); + + EXPECT_FALSE(report.stressComparisonApplicable); + EXPECT_NE(report.stressComparisonReason.find("N/A"), std::string::npos); + EXPECT_NE(report.stressComparisonReason.find("HDF5"), std::string::npos); + EXPECT_TRUE(report.physicsEvidence.endpointConsistencyPassed); + EXPECT_TRUE(std::isfinite(report.physicsEvidence.freeResidualNorm)); + EXPECT_LE(report.physicsEvidence.freeResidualNorm, 1.0e-3); + EXPECT_LE( + norm(sum( + report.physicsEvidence.appliedForce, + report.physicsEvidence.reactionForce)), + 1.0e-3); + EXPECT_LE( + norm(sum( + report.physicsEvidence.appliedMomentAboutOrigin, + report.physicsEvidence.reactionMomentAboutOrigin)), + 1.0e-2); + + ASSERT_TRUE( + fesa::test::ReferenceComparison::writeDeterministicJson( + report, comparison) + .isOk()); + ASSERT_TRUE(std::filesystem::is_regular_file(comparison)); + const std::string json = readBytes(comparison); + EXPECT_NE(json.find("\"stress_comparison_applicable\":false"), + std::string::npos); + EXPECT_NE(json.find("\"physics_evidence\""), std::string::npos); + + std::vector generatedNames; + for (const auto& entry : + std::filesystem::directory_iterator{outputDirectory}) { + generatedNames.push_back(entry.path().filename().string()); + } + std::sort(generatedNames.begin(), generatedNames.end()); + EXPECT_EQ( + generatedNames, + (std::vector{"comparison.json", "results.h5"})); + + expectTreeUnchanged(referenceBefore, snapshotTree(referenceDirectory)); +} diff --git a/tests/reference/reference_comparison.cpp b/tests/reference/reference_comparison.cpp new file mode 100644 index 0000000..bc8c6ec --- /dev/null +++ b/tests/reference/reference_comparison.cpp @@ -0,0 +1,1510 @@ +#include "reference_comparison.hpp" + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/assembly/load_assembler.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/io/abaqus/domain_mapper.hpp" +#include "fesa/io/abaqus/input_reader.hpp" +#include "fesa/results/result_recovery.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa::test { +namespace { + +constexpr const char* 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::model, + {{Severity::error, code, {}, "", kModelId, message}}); +} + +std::string trim(const std::string& value) { + const auto isSpace = [](const unsigned char character) { + return std::isspace(character) != 0; + }; + const auto begin = std::find_if_not( + value.begin(), value.end(), [&](const char character) { + return isSpace(static_cast(character)); + }); + const auto end = std::find_if_not( + value.rbegin(), value.rend(), [&](const char character) { + return isSpace(static_cast(character)); + }).base(); + return begin < end ? std::string{begin, end} : std::string{}; +} + +std::string collapseWhitespace(const std::string& value) { + std::string result; + bool pendingSpace = false; + for (const char character : trim(value)) { + if (std::isspace(static_cast(character)) != 0) { + pendingSpace = !result.empty(); + } else { + if (pendingSpace) { + result.push_back(' '); + } + result.push_back(character); + pendingSpace = 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(character - 'A' + 'a'); + } + return character; + }); + return value; +} + +std::vector splitCsvLine(const std::string& line) { + std::vector fields; + std::size_t start = 0U; + while (true) { + const std::size_t comma = line.find(',', start); + fields.push_back(trim(line.substr(start, comma - start))); + if (comma == std::string::npos) { + break; + } + start = comma + 1U; + } + return fields; +} + +std::int64_t parsePositiveLabel(const std::string& field) { + std::int64_t value = 0; + const char* const begin = field.data(); + const char* const end = begin + field.size(); + const auto parsed = std::from_chars(begin, end, value); + if (parsed.ec != std::errc{} || parsed.ptr != end || value <= 0) { + fail("schema-mismatch", "A 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 instanceName; + std::int64_t sourceNodeLabel; + std::vector values; +}; + +struct ReferenceTable { + std::vector rows; +}; + +ReferenceTable readReferenceCsv( + const std::filesystem::path& path, + const std::vector& expectedHeader) { + 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) != expectedHeader) { + 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() != expectedHeader.size() || + collapseWhitespace(fields[0U]) != kFrameText || + fields[1U].empty()) { + fail("schema-mismatch", "A reference CSV row has invalid schema or frame identity."); + } + WideReferenceRow row{}; + row.instanceName = fields[1U]; + row.sourceNodeLabel = 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.instanceName) == + asciiLower(row.instanceName) && + existing.sourceNodeLabel == row.sourceNodeLabel; + }); + 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& legacyDirectory) { + std::error_code error; + if (!std::filesystem::is_directory(legacyDirectory, error) || error) { + fail("needs-reference-artifacts", "The approved legacy directory is missing."); + } + std::vector names; + for (std::filesystem::directory_iterator iterator{legacyDirectory, 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 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& inputPath) { + AbaqusInputReader reader; + auto parsed = reader.read(inputPath); + 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_, &clientData_) >= 0 && + H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr) >= 0) { + active_ = true; + } + } + Hdf5ErrorSilencer(const Hdf5ErrorSilencer&) = delete; + Hdf5ErrorSilencer& operator=(const Hdf5ErrorSilencer&) = delete; + ~Hdf5ErrorSilencer() { + if (active_) { + (void)H5Eset_auto2(H5E_DEFAULT, callback_, clientData_); + } + } + +private: + H5E_auto2_t callback_{nullptr}; + void* clientData_{nullptr}; + bool active_{false}; +}; + +class Hdf5VlenReclaimer { +public: + Hdf5VlenReclaimer( + const hid_t memoryType, + const hid_t dataSpace, + void* const data) noexcept + : memoryType_{memoryType}, dataSpace_{dataSpace}, data_{data} {} + Hdf5VlenReclaimer(const Hdf5VlenReclaimer&) = delete; + Hdf5VlenReclaimer& operator=(const Hdf5VlenReclaimer&) = delete; + ~Hdf5VlenReclaimer() { + if (active_) { + (void)H5Dvlen_reclaim( + memoryType_, dataSpace_, H5P_DEFAULT, data_); + } + } + + herr_t reclaim() noexcept { + active_ = false; + return H5Dvlen_reclaim( + memoryType_, dataSpace_, H5P_DEFAULT, data_); + } + +private: + hid_t memoryType_; + hid_t dataSpace_; + void* data_; + bool active_{true}; +}; + +hid_t requireId(const hid_t value, const char* message) { + if (value < 0) { + fail("schema-mismatch", message); + } + return value; +} + +void requireHdf(const herr_t value, const char* message) { + if (value < 0) { + fail("schema-mismatch", message); + } +} + +Hdf5Handle openDataset(const hid_t file, const char* path) { + return {requireId(H5Dopen2(file, path, H5P_DEFAULT), + "A required HDF5 dataset is missing."), + H5Dclose}; +} + +std::vector 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 dimensions(static_cast(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* coordinateSystem, + const char* location) { + requireStringAttribute(dataset, "component_names", components); + requireStringAttribute(dataset, "component_unit_dimensions", units); + requireStringAttribute(dataset, "coordinate_system", coordinateSystem); + 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 readDoubleDataset( + const hid_t file, + const char* path, + const std::vector& expectedDimensions, + const char* components, + const char* units, + const char* coordinateSystem, + const char* location) { + auto dataset = openDataset(file, path); + if (datasetDimensions(dataset.get()) != expectedDimensions) { + 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, coordinateSystem, location); + std::size_t count = 1U; + for (const hsize_t dimension : expectedDimensions) { + if (dimension > (std::numeric_limits::max)() / count) { + fail("schema-mismatch", "An HDF5 result shape overflows size_t."); + } + count *= static_cast(dimension); + } + std::vector 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& expected) { + Hdf5Handle type{ + requireId(H5Dget_type(dataset), "Unable to inspect an HDF5 compound type."), + H5Tclose}; + if (H5Tget_class(type.get()) != H5T_COMPOUND || + H5Tget_nmembers(type.get()) != static_cast(expected.size())) { + fail("schema-mismatch", "An HDF5 compound dataset has the wrong schema."); + } + for (std::size_t index = 0U; index < expected.size(); ++index) { + char* raw = H5Tget_member_name(type.get(), static_cast(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 internalNodeId; + char* instanceName; + char* sourceLabel; + double coordinates[3]; +}; + +struct ElementMemoryRow { + std::uint64_t internalElementId; + char* instanceName; + char* sourceLabel; + std::uint64_t nodeInternalIds[2]; + double localAxes[9]; +}; + +struct HdfNode { + std::uint64_t internalNodeId; + std::string instanceName; + std::int64_t sourceNodeLabel; + std::string sourceNodeLabelText; + std::array coordinates; +}; + +struct HdfElement { + std::uint64_t internalElementId; + std::string instanceName; + std::int64_t sourceElementLabel; + std::string sourceElementLabelText; + std::array nodeInternalIds; + std::array localAxes; +}; + +std::vector 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 stringType = makeVariableStringType(); + const hsize_t coordinateDimensions[1] = {3U}; + Hdf5Handle coordinateType{ + requireId( + H5Tarray_create2(H5T_NATIVE_DOUBLE, 1, coordinateDimensions), + "Unable to create the node coordinate memory type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireId(H5Tcreate(H5T_COMPOUND, sizeof(NodeMemoryRow)), + "Unable to create the node memory type."), + H5Tclose}; + requireHdf( + H5Tinsert(memoryType.get(), "internal_node_id", + HOFFSET(NodeMemoryRow, internalNodeId), H5T_NATIVE_UINT64), + "Unable to define the node ID memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "instance_name", + HOFFSET(NodeMemoryRow, instanceName), stringType.get()), + "Unable to define the node instance memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "source_label", + HOFFSET(NodeMemoryRow, sourceLabel), stringType.get()), + "Unable to define the node label memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "coordinates", + HOFFSET(NodeMemoryRow, coordinates), coordinateType.get()), + "Unable to define the node coordinate memory field."); + + std::vector raw(static_cast(dimensions[0U])); + Hdf5Handle space{ + requireId(H5Dget_space(dataset.get()), + "Unable to reopen the node dataspace."), + H5Sclose}; + requireHdf( + H5Dread(dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL, + H5P_DEFAULT, raw.data()), + "Unable to read the HDF5 node table."); + Hdf5VlenReclaimer strings{memoryType.get(), space.get(), raw.data()}; + std::vector rows; + rows.reserve(raw.size()); + for (const auto& row : raw) { + if (row.instanceName == nullptr || row.sourceLabel == nullptr) { + fail("schema-mismatch", "An HDF5 node identity is null."); + } + const std::array 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.internalNodeId, + row.instanceName, + parsePositiveLabel(row.sourceLabel), + row.sourceLabel, + coordinates}); + } + requireHdf( + strings.reclaim(), + "Unable to reclaim HDF5 node strings."); + return rows; +} + +std::vector 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 stringType = makeVariableStringType(); + const hsize_t connectivityDimensions[1] = {2U}; + const hsize_t axesDimensions[2] = {3U, 3U}; + Hdf5Handle connectivityType{ + requireId( + H5Tarray_create2(H5T_NATIVE_UINT64, 1, connectivityDimensions), + "Unable to create the connectivity memory type."), + H5Tclose}; + Hdf5Handle axesType{ + requireId( + H5Tarray_create2(H5T_NATIVE_DOUBLE, 2, axesDimensions), + "Unable to create the local-axis memory type."), + H5Tclose}; + Hdf5Handle memoryType{ + requireId(H5Tcreate(H5T_COMPOUND, sizeof(ElementMemoryRow)), + "Unable to create the element memory type."), + H5Tclose}; + requireHdf( + H5Tinsert(memoryType.get(), "internal_element_id", + HOFFSET(ElementMemoryRow, internalElementId), + H5T_NATIVE_UINT64), + "Unable to define the element ID memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "instance_name", + HOFFSET(ElementMemoryRow, instanceName), stringType.get()), + "Unable to define the element instance memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "source_label", + HOFFSET(ElementMemoryRow, sourceLabel), stringType.get()), + "Unable to define the element label memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "node_internal_ids", + HOFFSET(ElementMemoryRow, nodeInternalIds), + connectivityType.get()), + "Unable to define the connectivity memory field."); + requireHdf( + H5Tinsert(memoryType.get(), "local_axes", + HOFFSET(ElementMemoryRow, localAxes), axesType.get()), + "Unable to define the local-axis memory field."); + + std::vector raw(static_cast(dimensions[0U])); + Hdf5Handle space{ + requireId(H5Dget_space(dataset.get()), + "Unable to reopen the element dataspace."), + H5Sclose}; + requireHdf( + H5Dread(dataset.get(), memoryType.get(), H5S_ALL, H5S_ALL, + H5P_DEFAULT, raw.data()), + "Unable to read the HDF5 element table."); + Hdf5VlenReclaimer strings{memoryType.get(), space.get(), raw.data()}; + std::vector rows; + rows.reserve(raw.size()); + for (const auto& row : raw) { + if (row.instanceName == nullptr || row.sourceLabel == nullptr) { + fail("schema-mismatch", "An HDF5 element identity is null."); + } + HdfElement converted{}; + converted.internalElementId = row.internalElementId; + converted.instanceName = row.instanceName; + converted.sourceElementLabel = parsePositiveLabel(row.sourceLabel); + converted.sourceElementLabelText = row.sourceLabel; + converted.nodeInternalIds = { + row.nodeInternalIds[0U], row.nodeInternalIds[1U]}; + std::copy( + std::begin(row.localAxes), std::end(row.localAxes), + converted.localAxes.begin()); + if (!std::all_of( + converted.localAxes.begin(), converted.localAxes.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 memoryType{ + requireId(H5Tcreate(H5T_COMPOUND, sizeof(StressValue)), + "Unable to create a stress memory type."), + H5Tclose}; + requireHdf( + H5Tinsert(memoryType.get(), "S11", HOFFSET(StressValue, s11), + H5T_NATIVE_DOUBLE), + "Unable to define the stress memory field."); + std::vector values(static_cast(dimensions[0U])); + requireHdf( + H5Dread(dataset.get(), memoryType.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 expectedLocalAxes( + const Domain& domain, const EulerBeam3DDefinition& element) { + const auto& first = domain.nodes()[element.nodeIndices[0U]].coordinates; + const auto& second = domain.nodes()[element.nodeIndices[1U]].coordinates; + const auto& guide = domain.sections()[element.sectionIndex].firstAxis; + const std::array 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 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 yTrial = { + guide[0U] - projection * x[0U], + guide[1U] - projection * x[1U], + guide[2U] - projection * x[2U]}; + const double yNorm = std::hypot(yTrial[0U], yTrial[1U], yTrial[2U]); + const std::array y = { + yTrial[0U] / yNorm, yTrial[1U] / yNorm, yTrial[2U] / yNorm}; + const std::array 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 nodes; + std::vector elements; + std::vector displacement; + std::vector reaction; + std::vector sectionResultants; +}; + +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 sourceIdentity = + readStringAttribute(metadata.get(), "source_input_identity"); + const std::string normalizedInput = + std::filesystem::absolute(input).lexically_normal().generic_u8string(); + const std::string expectedIdentity = + "path=" + normalizedInput + + ";content_identity=" + domain.sourceContentIdentity(); + if (sourceIdentity != expectedIdentity) { + 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.internalNodeId != node || + actual.instanceName != expected.sourceId.instanceName || + actual.sourceNodeLabel != expected.sourceId.sourceLabel || + actual.sourceNodeLabelText != expected.sourceId.sourceLabelText || + 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.internalElementId != element || + actual.instanceName != expected.sourceId.instanceName || + actual.sourceElementLabel != expected.sourceId.sourceLabel || + actual.sourceElementLabelText != expected.sourceId.sourceLabelText || + actual.nodeInternalIds[0U] != expected.nodeIndices[0U] || + actual.nodeInternalIds[1U] != expected.nodeIndices[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.localAxes[component] - axes[component]) > 1.0e-12) { + fail("schema-mismatch", "An HDF5 local axis does not match the input."); + } + } + } + + const hsize_t nodeCount = static_cast(projection.nodes.size()); + const hsize_t elementCount = static_cast(projection.elements.size()); + projection.displacement = readDoubleDataset( + file.get(), kDisplacementPath, {nodeCount, 6U}, + "UX,UY,UZ,URX,URY,URZ", + "length,length,length,radian,radian,radian", + "global-cartesian", "nodal"); + projection.reaction = readDoubleDataset( + file.get(), kReactionPath, {nodeCount, 6U}, + "RF1,RF2,RF3,RM1,RM2,RM3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "nodal"); + projection.sectionResultants = readDoubleDataset( + file.get(), kSectionPath, {elementCount, 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 orderedRows( + const ReferenceTable& table, const std::vector& nodes) { + if (table.rows.size() != nodes.size()) { + fail("schema-mismatch", "The FESA/reference projected row sets differ."); + } + std::vector 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.instanceName) == + asciiLower(node.instanceName) && + row.sourceNodeLabel == node.sourceNodeLabel; + }); + if (found == table.rows.end() || found->instanceName != node.instanceName) { + 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 valueIndex) { + double scale = 0.0; + for (const auto& row : table.rows) { + if (valueIndex >= row.values.size()) { + fail("schema-mismatch", "A reference row has the wrong component arity."); + } + scale = (std::max)(scale, std::abs(row.values[valueIndex])); + } + return scale; +} + +std::vector normalizeStations( + const Domain& domain, + const HdfProjection& hdf, + const ReferenceTable& sectionTable) { + auto modelResult = AnalysisModel::create(domain); + if (!modelResult.hasValue()) { + fail("schema-mismatch", "The approved input cannot create an analysis view."); + } + const AnalysisModel model = std::move(modelResult.value()); + const std::array tolerances = { + kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 0U), + kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 3U), + kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 1U), + kForceMomentFloor + kRelativeCoefficient * tableScale(sectionTable, 2U)}; + std::vector 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( + hdf.elements[element].nodeInternalIds[endpoint]); + std::array values{}; + for (std::size_t component = 0U; component < values.size(); ++component) { + values[component] = + hdf.sectionResultants[(element * 2U + endpoint) * 4U + component]; + } + endpoints.push_back({ + static_cast(element), + static_cast(endpoint), + domain.nodes()[node].sourceId, + {}, + values}); + } + } + auto normalized = ResultRecovery::normalizeSectionResultantsToNodeStations( + model, endpoints, tolerances); + if (!normalized.hasValue()) { + const auto& diagnostics = normalized.status().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& stations, + const HdfNode& node) { + const auto found = std::find_if( + stations.begin(), stations.end(), [&](const NodeStationResultRow& row) { + return row.node.instanceName == node.instanceName && + row.node.sourceLabel == node.sourceNodeLabel; + }); + 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 coordinateSystem, + std::string datasetPath) { + return { + kModelId, + kStepName, + kFrameIndex, + node.instanceName, + node.sourceNodeLabel, + quantity, + std::move(component), + value, + std::move(unit), + std::move(coordinateSystem), + std::move(datasetPath)}; +} + +void appendNodalRows( + ComparisonReport& report, + const HdfProjection& hdf, + const std::vector& reference, + const ComparisonQuantity quantity, + const std::array& components, + const std::array& units, + const std::vector& fesaValues, + const char* datasetPath) { + 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], + fesaValues[node * 6U + component], units[component], + "global-cartesian", datasetPath); + auto abaqus = canonicalRow( + hdf.nodes[node], quantity, components[component], + reference[node]->values[component], units[component], + "global-cartesian", datasetPath); + 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& reference, + const std::vector& stations) { + const std::array components = {"N", "T", "My", "Mz"}; + const std::array units = { + "force", "force*length", "force*length", "force*length"}; + const std::array referenceColumns = {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::sectionResultant, + components[component], station.sectionResultant[component], + units[component], "beam-local", kSectionPath); + auto abaqus = canonicalRow( + hdf.nodes[node], ComparisonQuantity::sectionResultant, + components[component], + reference[node]->values[referenceColumns[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::displacement + ? kKinematicFloor + : kForceMomentFloor; +} + +void evaluateGroup( + ComparisonReport& report, + const ComparisonQuantity quantity, + const std::string& component) { + std::vector rowIndices; + for (std::size_t index = 0U; index < report.rows.size(); ++index) { + if (report.rows[index].reference.quantity == quantity && + report.rows[index].reference.component == component) { + rowIndices.push_back(index); + } + } + if (rowIndices.empty()) { + fail("schema-mismatch", "A canonical comparison component has no rows."); + } + double referenceScale = 0.0; + for (const std::size_t index : rowIndices) { + referenceScale = (std::max)( + referenceScale, std::abs(report.rows[index].reference.value)); + } + const double tolerance = + absoluteFloor(quantity, component) + + kRelativeCoefficient * referenceScale; + double maximumAbsolute = -1.0; + double maximumNormalized = 0.0; + std::size_t worstRow = rowIndices.front(); + long double squaredError = 0.0L; + for (const std::size_t index : rowIndices) { + auto& row = report.rows[index]; + row.absoluteError = std::abs(row.fesa.value - row.reference.value); + if (!std::isfinite(row.absoluteError)) { + fail("schema-mismatch", "A canonical row error is nonfinite."); + } + row.tolerance = tolerance; + row.passed = row.absoluteError <= tolerance; + report.passed = report.passed && row.passed; + const double normalized = row.absoluteError / tolerance; + if (row.absoluteError > maximumAbsolute) { + maximumAbsolute = row.absoluteError; + worstRow = index; + } + maximumNormalized = (std::max)(maximumNormalized, normalized); + const long double error = static_cast(row.absoluteError); + squaredError += error * error; + } + const double normError = std::sqrt(static_cast(squaredError)); + const double rmsError = std::sqrt( + static_cast(squaredError / + static_cast(rowIndices.size()))); + if (!std::isfinite(normError) || !std::isfinite(rmsError)) { + fail("schema-mismatch", "A component aggregate error is nonfinite."); + } + report.metrics.push_back({ + quantity, + component, + referenceScale, + maximumAbsolute, + maximumNormalized, + rmsError, + normError, + worstRow}); +} + +PhysicsEvidence makePhysicsEvidence( + const Domain& domain, + const HdfProjection& hdf) { + auto modelResult = AnalysisModel::create(domain); + if (!modelResult.hasValue()) { + fail("schema-mismatch", "The approved input cannot create physics evidence."); + } + const AnalysisModel model = std::move(modelResult.value()); + auto dofsResult = DofManager::create(model); + if (!dofsResult.hasValue()) { + fail("schema-mismatch", "The approved input cannot create a DOF map."); + } + const DofManager dofs = std::move(dofsResult.value()); + auto loadResult = LoadAssembler::assembleFullNodalLoad(model, dofs); + if (!loadResult.hasValue()) { + fail("schema-mismatch", "The approved input load cannot be assembled."); + } + const Vector load = std::move(loadResult.value()); + if (load.size() != hdf.reaction.size()) { + fail("schema-mismatch", "The load and reaction spaces are inconsistent."); + } + + PhysicsEvidence evidence{}; + long double residualSquared = 0.0L; + for (const std::size_t freeDof : dofs.freeDofs()) { + const long double value = + static_cast(hdf.reaction[freeDof]); + residualSquared += value * value; + } + evidence.freeResidualNorm = + std::sqrt(static_cast(residualSquared)); + for (std::size_t node = 0U; node < domain.nodes().size(); ++node) { + const auto& coordinates = domain.nodes()[node].coordinates; + const std::array applied = { + load[node * 6U + 0U], + load[node * 6U + 1U], + load[node * 6U + 2U]}; + const std::array 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.appliedForce[component] += applied[component]; + evidence.reactionForce[component] += reaction[component]; + evidence.appliedMomentAboutOrigin[component] += + load[node * 6U + 3U + component]; + evidence.reactionMomentAboutOrigin[component] += + hdf.reaction[node * 6U + 3U + component]; + } + evidence.appliedMomentAboutOrigin[0U] += + coordinates[1U] * applied[2U] - coordinates[2U] * applied[1U]; + evidence.appliedMomentAboutOrigin[1U] += + coordinates[2U] * applied[0U] - coordinates[0U] * applied[2U]; + evidence.appliedMomentAboutOrigin[2U] += + coordinates[0U] * applied[1U] - coordinates[1U] * applied[0U]; + evidence.reactionMomentAboutOrigin[0U] += + coordinates[1U] * reaction[2U] - coordinates[2U] * reaction[1U]; + evidence.reactionMomentAboutOrigin[1U] += + coordinates[2U] * reaction[0U] - coordinates[0U] * reaction[2U]; + evidence.reactionMomentAboutOrigin[2U] += + coordinates[0U] * reaction[1U] - coordinates[1U] * reaction[0U]; + } + evidence.endpointConsistencyPassed = true; + return evidence; +} + +const char* quantityName(const ComparisonQuantity quantity) { + switch (quantity) { + case ComparisonQuantity::displacement: + return "displacement"; + case ComparisonQuantity::reaction: + return "reaction"; + case ComparisonQuantity::sectionResultant: + return "section_resultant"; + } + return "unknown"; +} + +void writeJsonString(std::ostream& stream, const std::string& value) { + static constexpr char digits[] = "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" << digits[character >> 4U] + << digits[character & 0x0fU]; + } else { + stream.put(static_cast(character)); + } + break; + } + } + stream.put('"'); +} + +void writeCanonicalRow( + std::ostream& stream, const CanonicalComparisonRow& row) { + stream << "{\"model_id\":"; + writeJsonString(stream, row.modelId); + stream << ",\"step_name\":"; + writeJsonString(stream, row.stepName); + stream << ",\"frame_index\":" << row.frameIndex + << ",\"instance_name\":"; + writeJsonString(stream, row.instanceName); + stream << ",\"source_node_label\":" << row.sourceNodeLabel + << ",\"quantity\":"; + writeJsonString(stream, quantityName(row.quantity)); + stream << ",\"component\":"; + writeJsonString(stream, row.component); + stream << ",\"value\":" << row.value << ",\"unit_dimension\":"; + writeJsonString(stream, row.unitDimension); + stream << ",\"coordinate_system\":"; + writeJsonString(stream, row.coordinateSystem); + stream << ",\"hdf5_dataset_path\":"; + writeJsonString(stream, row.hdf5DatasetPath); + stream << '}'; +} + +void writeArray(std::ostream& stream, const std::array& values) { + stream << '[' << values[0U] << ',' << values[1U] << ',' << values[2U] + << ']'; +} + +bool finiteReport(const ComparisonReport& report) { + const auto finiteArray = [](const std::array& values) { + return std::all_of(values.begin(), values.end(), [](const double value) { + return std::isfinite(value); + }); + }; + if (!std::isfinite(report.physicsEvidence.freeResidualNorm) || + !finiteArray(report.physicsEvidence.appliedForce) || + !finiteArray(report.physicsEvidence.reactionForce) || + !finiteArray(report.physicsEvidence.appliedMomentAboutOrigin) || + !finiteArray(report.physicsEvidence.reactionMomentAboutOrigin)) { + return false; + } + for (const auto& row : report.rows) { + if (!std::isfinite(row.fesa.value) || + !std::isfinite(row.reference.value) || + !std::isfinite(row.absoluteError) || + !std::isfinite(row.tolerance)) { + return false; + } + } + return std::all_of( + report.metrics.begin(), report.metrics.end(), + [](const ComponentMetrics& metric) { + return std::isfinite(metric.referenceScale) && + std::isfinite(metric.maximumAbsoluteError) && + std::isfinite(metric.maximumNormalizedError) && + std::isfinite(metric.rmsError) && + std::isfinite(metric.normError); + }); +} + +} // namespace + +Result ReferenceComparison::compare( + const std::filesystem::path& resultsHdf5, + const std::filesystem::path& legacyReferenceDirectory) { + try { + requireExactArtifactInventory(legacyReferenceDirectory); + const auto input = legacyReferenceDirectory / kInputName; + Domain domain = readApprovedDomain(input); + const ReferenceTable displacement = readReferenceCsv( + legacyReferenceDirectory / kDisplacementName, + {"Frame", "Part Instance Name", "Node Label", "U-U1", "U-U2", + "U-U3", "UR-UR1", "UR-UR2", "UR-UR3"}); + const ReferenceTable reaction = readReferenceCsv( + legacyReferenceDirectory / kReactionName, + {"Frame", "Part Instance Name", "Node Label", "RF-RF1", "RF-RF2", + "RF-RF3", "RM-RM1", "RM-RM2", "RM-RM3"}); + const ReferenceTable section = readReferenceCsv( + legacyReferenceDirectory / kSectionName, + {"Frame", "Part Instance Name", "Node Label", "SF-SF1", "SM-SM1", + "SM-SM2", "SM-SM3"}); + HdfProjection hdf = readHdfProjection(resultsHdf5, input, domain); + const auto displacementRows = orderedRows(displacement, hdf.nodes); + const auto reactionRows = orderedRows(reaction, hdf.nodes); + const auto sectionRows = 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, + displacementRows, + ComparisonQuantity::displacement, + {"UX", "UY", "UZ", "URX", "URY", "URZ"}, + {"length", "length", "length", "radian", "radian", "radian"}, + hdf.displacement, + kDisplacementPath); + appendNodalRows( + report, + hdf, + reactionRows, + ComparisonQuantity::reaction, + {"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"}, + {"force", "force", "force", "force*length", "force*length", + "force*length"}, + hdf.reaction, + kReactionPath); + appendSectionRows(report, hdf, sectionRows, stations); + + for (const std::string& component : + {"UX", "UY", "UZ", "URX", "URY", "URZ"}) { + evaluateGroup(report, ComparisonQuantity::displacement, component); + } + for (const std::string& component : + {"RF1", "RF2", "RF3", "RM1", "RM2", "RM3"}) { + evaluateGroup(report, ComparisonQuantity::reaction, component); + } + for (const std::string& component : {"N", "T", "My", "Mz"}) { + evaluateGroup(report, ComparisonQuantity::sectionResultant, component); + } + report.physicsEvidence = makePhysicsEvidence(domain, hdf); + report.stressComparisonApplicable = false; + report.stressComparisonReason = + "Abaqus beam stress comparison is N/A; analytical/unit and HDF5 " + "schema tests provide stress evidence."; + return Result::success(std::move(report)); + } catch (const ComparisonFailure& failure) { + return Result::failure( + comparisonFailureStatus(failure.code(), failure.what())); + } catch (const std::exception& failure) { + return Result::failure(comparisonFailureStatus( + "schema-mismatch", failure.what())); + } +} + +Status ReferenceComparison::writeDeterministicJson( + const ComparisonReport& report, + const std::filesystem::path& outputJson) { + if (outputJson.empty() || outputJson.filename().empty() || + !finiteReport(report)) { + return Status::failure( + FailureCategory::output, + {{Severity::error, + "comparison-json-write-failure", + {}, + "", + kModelId, + "The deterministic comparison report or output path is invalid."}}); + } + std::ofstream stream{outputJson, std::ios::binary | std::ios::trunc}; + if (!stream) { + return Status::failure( + FailureCategory::output, + {{Severity::error, + "comparison-json-write-failure", + {}, + "", + kModelId, + "The deterministic comparison JSON cannot be opened."}}); + } + stream.imbue(std::locale::classic()); + stream << std::setprecision(std::numeric_limits::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.absoluteError + << ",\"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.referenceScale + << ",\"maximum_absolute_error\":" + << metric.maximumAbsoluteError + << ",\"maximum_normalized_error\":" + << metric.maximumNormalizedError + << ",\"rms_error\":" << metric.rmsError + << ",\"norm_error\":" << metric.normError + << ",\"worst_row\":" << metric.worstRow << '}'; + } + stream << "],\"physics_evidence\":{\"free_residual_norm\":" + << report.physicsEvidence.freeResidualNorm + << ",\"applied_force\":"; + writeArray(stream, report.physicsEvidence.appliedForce); + stream << ",\"reaction_force\":"; + writeArray(stream, report.physicsEvidence.reactionForce); + stream << ",\"applied_moment_about_origin\":"; + writeArray(stream, report.physicsEvidence.appliedMomentAboutOrigin); + stream << ",\"reaction_moment_about_origin\":"; + writeArray(stream, report.physicsEvidence.reactionMomentAboutOrigin); + stream << ",\"endpoint_consistency_passed\":" + << (report.physicsEvidence.endpointConsistencyPassed ? "true" : "false") + << "},\"stress_comparison_applicable\":" + << (report.stressComparisonApplicable ? "true" : "false") + << ",\"stress_comparison_reason\":"; + writeJsonString(stream, report.stressComparisonReason); + stream << ",\"passed\":" << (report.passed ? "true" : "false") << "}\n"; + if (!stream) { + return Status::failure( + FailureCategory::output, + {{Severity::error, + "comparison-json-write-failure", + {}, + "", + kModelId, + "The deterministic comparison JSON write failed."}}); + } + return Status::ok(); +} + +} // namespace fesa::test diff --git a/tests/reference/reference_comparison.hpp b/tests/reference/reference_comparison.hpp new file mode 100644 index 0000000..28114e2 --- /dev/null +++ b/tests/reference/reference_comparison.hpp @@ -0,0 +1,83 @@ +#pragma once + +#include "fesa/core/status.hpp" + +#include +#include +#include +#include +#include +#include + +namespace fesa::test { + +enum class ComparisonQuantity { + displacement, + reaction, + sectionResultant +}; + +struct CanonicalComparisonRow { + std::string modelId; + std::string stepName; + std::size_t frameIndex; + std::string instanceName; + std::int64_t sourceNodeLabel; + ComparisonQuantity quantity; + std::string component; + double value; + std::string unitDimension; + std::string coordinateSystem; + std::string hdf5DatasetPath; +}; + +struct RowDecision { + CanonicalComparisonRow fesa; + CanonicalComparisonRow reference; + double absoluteError; + double tolerance; + bool passed; +}; + +struct ComponentMetrics { + ComparisonQuantity quantity; + std::string component; + double referenceScale; + double maximumAbsoluteError; + double maximumNormalizedError; + double rmsError; + double normError; + std::size_t worstRow; +}; + +struct PhysicsEvidence { + double freeResidualNorm; + std::array appliedForce; + std::array reactionForce; + std::array appliedMomentAboutOrigin; + std::array reactionMomentAboutOrigin; + bool endpointConsistencyPassed; +}; + +struct ComparisonReport { + std::vector rows; + std::vector metrics; + PhysicsEvidence physicsEvidence; + bool stressComparisonApplicable; + std::string stressComparisonReason; + bool passed; +}; + +// Test-only comparison support keeps Abaqus artifacts read-only and exposes no +// backend handles to the implementation or downstream verification Steps. +class ReferenceComparison { +public: + static Result compare( + const std::filesystem::path& resultsHdf5, + const std::filesystem::path& legacyReferenceDirectory); + static Status writeDeterministicJson( + const ComparisonReport& report, + const std::filesystem::path& outputJson); +}; + +} // namespace fesa::test diff --git a/tests/reference/reference_comparison_test.cpp b/tests/reference/reference_comparison_test.cpp new file mode 100644 index 0000000..3d430bf --- /dev/null +++ b/tests/reference/reference_comparison_test.cpp @@ -0,0 +1,730 @@ +#include "reference_comparison.hpp" + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/analysis/analysis_state.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/io/abaqus/input_reader.hpp" +#include "fesa/io/hdf5/hdf5_results_writer.hpp" +#include "fesa/model/domain.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef FESA_TEST_SOURCE_DIR +#error FESA_TEST_SOURCE_DIR must identify the repository root. +#endif + +#ifndef FESA_TEST_BINARY_DIR +#error FESA_TEST_BINARY_DIR must identify the CMake binary root. +#endif + +namespace { + +constexpr const char* 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* kInstanceName = "PART-1_1-1"; +constexpr std::size_t kNodeCount = 11U; +constexpr std::size_t kElementCount = 10U; +constexpr std::size_t kExpectedRowCount = 176U; +constexpr std::size_t kExpectedMetricCount = 16U; + +using NodalValues = std::array, kNodeCount>; +using EndpointValues = + std::array, 2>, kElementCount>; + +struct ComparisonValues { + NodalValues displacement{}; + NodalValues reaction{}; + EndpointValues sectionResultants{}; +}; + +const std::filesystem::path& sourceRoot() { + static const std::filesystem::path root{FESA_TEST_SOURCE_DIR}; + return root; +} + +const std::filesystem::path& binaryRoot() { + static const std::filesystem::path root{FESA_TEST_BINARY_DIR}; + return root; +} + +std::string readBytes(const std::filesystem::path& path) { + std::ifstream stream{path, std::ios::binary}; + if (!stream) { + throw std::runtime_error{"Unable to read fixture: " + path.string()}; + } + return {std::istreambuf_iterator{stream}, + std::istreambuf_iterator{}}; +} + +void writeBytes(const std::filesystem::path& path, const std::string& contents) { + std::ofstream stream{path, std::ios::binary | std::ios::trunc}; + if (!stream) { + throw std::runtime_error{"Unable to write build-local fixture: " + + path.string()}; + } + stream.write(contents.data(), static_cast(contents.size())); + if (!stream) { + throw std::runtime_error{"Unable to finish build-local fixture write."}; + } +} + +std::vector readLines(const std::filesystem::path& path) { + std::ifstream stream{path}; + if (!stream) { + throw std::runtime_error{"Unable to read fixture lines."}; + } + std::vector lines; + for (std::string line; std::getline(stream, line);) { + if (!line.empty() && line.back() == '\r') { + line.pop_back(); + } + lines.push_back(std::move(line)); + } + return lines; +} + +void writeLines( + const std::filesystem::path& path, + const std::vector& lines) { + std::ofstream stream{path, std::ios::trunc}; + if (!stream) { + throw std::runtime_error{"Unable to write build-local fixture lines."}; + } + for (const auto& line : lines) { + stream << line << '\n'; + } + if (!stream) { + throw std::runtime_error{"Unable to finish build-local line write."}; + } +} + +void replaceFirst( + std::string& contents, + const std::string& from, + const std::string& to) { + const std::size_t position = contents.find(from); + if (position == std::string::npos) { + throw std::runtime_error{"Fixture token was not found: " + from}; + } + contents.replace(position, from.size(), to); +} + +ComparisonValues referenceValues() { + ComparisonValues values{}; + const std::array uz = { + -1.0e-30, + -2.761905780e-4, + -1.066667140e-3, + -2.314286540e-3, + -3.961906300e-3, + -5.952383390e-3, + -8.228574880e-3, + -1.073333810e-2, + -1.340952890e-2, + -1.620000600e-2, + -1.904762720e-2}; + const std::array ury = { + 1.0e-29, + 5.428573350e-4, + 1.028571860e-3, + 1.457143460e-3, + 1.828572130e-3, + 2.142857990e-3, + 2.400001050e-3, + 2.600000940e-3, + 2.742858140e-3, + 2.828572640e-3, + 2.857143990e-3}; + std::array, kNodeCount> stations{}; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + values.displacement[node][2U] = uz[node]; + values.displacement[node][4U] = ury[node]; + stations[node][2U] = node < kElementCount + ? 1.0e7 - 1.0e6 * static_cast(node) + : -1.56e-2; + } + values.reaction[0U][2U] = 1.0e6; + values.reaction[0U][4U] = -1.0e7; + for (std::size_t element = 0U; element < kElementCount; ++element) { + values.sectionResultants[element][0U] = stations[element]; + values.sectionResultants[element][1U] = stations[element + 1U]; + } + return values; +} + +fesa::ModelDefinition makeDefinition( + const std::filesystem::path& input, + std::string sourceContentIdentity) { + fesa::ModelDefinition definition{}; + definition.sourcePath = input; + definition.sourceContentIdentity = std::move(sourceContentIdentity); + for (std::size_t node = 0U; node < kNodeCount; ++node) { + const auto label = static_cast(node + 1U); + definition.nodes.push_back({ + {kInstanceName, label, std::to_string(label)}, + {static_cast(node), 0.0, 0.0}, + {input, node + 1U}}); + } + definition.materials.push_back( + {"Material-1", 2.1e11, 0.3, {input, 20U}}); + definition.sections.push_back({ + "Section-1", + 1.0, + 0.0833333, + 0.0, + 0.0833333, + 0.140833, + {0.0, 1.0, 0.0}, + {}, + {input, 30U}}); + for (std::size_t element = 0U; element < kElementCount; ++element) { + const auto label = static_cast(element + 1U); + definition.elements.push_back({ + {kInstanceName, label, std::to_string(label)}, + {static_cast(element), + static_cast(element + 1U)}, + 0U, + 0U, + {input, 40U + element}}); + } + definition.steps.push_back( + {"Step-1", {}, {}, 1.0, 1.0, 1.0e-5, 1.0, {input, 60U}}); + return definition; +} + +void writeResultsFixture( + const std::filesystem::path& output, + const std::filesystem::path& input, + const ComparisonValues& values) { + auto parsedInput = fesa::AbaqusInputReader{}.read(input); + if (!parsedInput.hasValue()) { + throw std::runtime_error{"Reference fixture input identity read failed."}; + } + auto domainResult = fesa::Domain::create( + makeDefinition(input, parsedInput.value().sourceContentIdentity)); + if (!domainResult.hasValue()) { + throw std::runtime_error{"Reference fixture Domain construction failed."}; + } + fesa::Domain domain = std::move(domainResult.value()); + auto modelResult = fesa::AnalysisModel::create(domain); + if (!modelResult.hasValue()) { + throw std::runtime_error{"Reference fixture AnalysisModel construction failed."}; + } + fesa::AnalysisModel model = std::move(modelResult.value()); + auto dofsResult = fesa::DofManager::create(model); + if (!dofsResult.hasValue()) { + throw std::runtime_error{"Reference fixture DofManager construction failed."}; + } + fesa::DofManager dofs = std::move(dofsResult.value()); + fesa::AnalysisState state = + fesa::AnalysisState::create(dofs, {"Step-1", 0U}); + + for (std::size_t node = 0U; node < kNodeCount; ++node) { + for (std::size_t component = 0U; component < 6U; ++component) { + const std::size_t index = node * 6U + component; + state.displacement()[index] = values.displacement[node][component]; + state.reaction()[index] = values.reaction[node][component]; + state.residual()[index] = values.reaction[node][component]; + } + } + for (std::size_t element = 0U; element < kElementCount; ++element) { + for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) { + const std::size_t node = element + endpoint; + state.endpointResults().push_back({ + static_cast(element), + static_cast(endpoint), + domain.nodes()[node].sourceId, + {}, + values.sectionResultants[element][endpoint]}); + } + state.gaussResults().push_back( + {static_cast(element), 1, {}, {}}); + state.gaussResults().push_back( + {static_cast(element), 2, {}, {}}); + state.stressResults().push_back({ + static_cast(element), + 1, + 0U, + 0.0, + 0.0, + 0.0, + "fesa-default"}); + state.stressResults().push_back({ + static_cast(element), + 2, + 0U, + 0.0, + 0.0, + 0.0, + "fesa-default"}); + } + + fesa::Hdf5ResultsWriter writer; + const fesa::Status status = writer.write(output, domain, state, {}); + if (!status.isOk()) { + throw std::runtime_error{"Reference fixture HDF5 write failed."}; + } +} + +class ContractFixture { +public: + ContractFixture(std::string label, const ComparisonValues& values) { + static std::atomic sequence{0U}; + root_ = binaryRoot() / "reference" / "contract-fixtures" / + (std::move(label) + "-" + + std::to_string(sequence.fetch_add(1U))); + legacy_ = root_ / "cantilever beam"; + std::error_code error; + std::filesystem::remove_all(root_, error); + error.clear(); + if (!std::filesystem::create_directories(legacy_, error) || error) { + throw std::runtime_error{"Unable to create contract fixture directory."}; + } + const auto approved = sourceRoot() / "reference" / "cantilever beam"; + for (const char* name : + {kInputName, kDisplacementName, kReactionName, kSectionName}) { + std::filesystem::copy_file( + approved / name, + legacy_ / name, + std::filesystem::copy_options::overwrite_existing); + } + results_ = root_ / "results.h5"; + writeResultsFixture(results_, legacy_ / kInputName, values); + } + + ContractFixture(const ContractFixture&) = delete; + ContractFixture& operator=(const ContractFixture&) = delete; + + ~ContractFixture() { + std::error_code ignored; + std::filesystem::remove_all(root_, ignored); + } + + const std::filesystem::path& root() const noexcept { return root_; } + const std::filesystem::path& legacy() const noexcept { return legacy_; } + const std::filesystem::path& results() const noexcept { return results_; } + +private: + std::filesystem::path root_; + std::filesystem::path legacy_; + std::filesystem::path results_; +}; + +void expectFailureCode( + const fesa::Result& result, + const std::string& expectedCode) { + ASSERT_FALSE(result.hasValue()); + ASSERT_FALSE(result.status().isOk()); + ASSERT_FALSE(result.status().diagnostics().empty()); + EXPECT_EQ(result.status().diagnostics().front().code, expectedCode); +} + +const fesa::test::RowDecision* findRow( + const fesa::test::ComparisonReport& report, + const fesa::test::ComparisonQuantity quantity, + const std::int64_t sourceNodeLabel, + const std::string& component) { + const auto found = std::find_if( + report.rows.begin(), + report.rows.end(), + [&](const fesa::test::RowDecision& row) { + return row.reference.quantity == quantity && + row.reference.sourceNodeLabel == sourceNodeLabel && + row.reference.component == component; + }); + return found == report.rows.end() ? nullptr : &*found; +} + +const fesa::test::ComponentMetrics* findMetric( + const fesa::test::ComparisonReport& report, + const fesa::test::ComparisonQuantity quantity, + const std::string& component) { + const auto found = std::find_if( + report.metrics.begin(), + report.metrics.end(), + [&](const fesa::test::ComponentMetrics& metric) { + return metric.quantity == quantity && metric.component == component; + }); + return found == report.metrics.end() ? nullptr : &*found; +} + +void expectExactRowInventory(const fesa::test::ComparisonReport& report) { + ASSERT_EQ(report.rows.size(), kExpectedRowCount); + std::size_t rowIndex = 0U; + const auto expectRow = [&](const fesa::test::ComparisonQuantity quantity, + const std::size_t node, + const std::string& component, + const std::string& unit, + const std::string& coordinateSystem, + const std::string& datasetPath) { + ASSERT_LT(rowIndex, report.rows.size()); + const auto& row = report.rows[rowIndex++]; + for (const auto* side : {&row.fesa, &row.reference}) { + EXPECT_EQ(side->modelId, "cantilever-beam-b33"); + EXPECT_EQ(side->stepName, "Step-1"); + EXPECT_EQ(side->frameIndex, 0U); + EXPECT_EQ(side->instanceName, kInstanceName); + EXPECT_EQ( + side->sourceNodeLabel, + static_cast(node + 1U)); + EXPECT_EQ(side->quantity, quantity); + EXPECT_EQ(side->component, component); + EXPECT_EQ(side->unitDimension, unit); + EXPECT_EQ(side->coordinateSystem, coordinateSystem); + EXPECT_EQ(side->hdf5DatasetPath, datasetPath); + } + }; + + const std::array displacementComponents = { + "UX", "UY", "UZ", "URX", "URY", "URZ"}; + const std::array displacementUnits = { + "length", "length", "length", "radian", "radian", "radian"}; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + for (std::size_t component = 0U; + component < displacementComponents.size(); + ++component) { + expectRow( + fesa::test::ComparisonQuantity::displacement, + node, + displacementComponents[component], + displacementUnits[component], + "global-cartesian", + "/steps/Step-1/frames/0/nodal/displacement"); + } + } + + const std::array reactionComponents = { + "RF1", "RF2", "RF3", "RM1", "RM2", "RM3"}; + const std::array reactionUnits = { + "force", + "force", + "force", + "force*length", + "force*length", + "force*length"}; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + for (std::size_t component = 0U; + component < reactionComponents.size(); + ++component) { + expectRow( + fesa::test::ComparisonQuantity::reaction, + node, + reactionComponents[component], + reactionUnits[component], + "global-cartesian", + "/steps/Step-1/frames/0/nodal/reaction"); + } + } + + const std::array sectionComponents = { + "N", "T", "My", "Mz"}; + const std::array sectionUnits = { + "force", "force*length", "force*length", "force*length"}; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + for (std::size_t component = 0U; + component < sectionComponents.size(); + ++component) { + expectRow( + fesa::test::ComparisonQuantity::sectionResultant, + node, + sectionComponents[component], + sectionUnits[component], + "beam-local", + "/steps/Step-1/frames/0/element/section_resultant"); + } + } + EXPECT_EQ(rowIndex, report.rows.size()); +} + +} // namespace + +TEST(ReferenceComparisonContract, + PrecheckRejectsMissingSchemaDuplicateAndNonfiniteRows) { + auto mismatchedValues = referenceValues(); + mismatchedValues.displacement[0U][0U] = 1.0; + + { + ContractFixture fixture{"missing-file", mismatchedValues}; + ASSERT_TRUE(std::filesystem::remove( + fixture.legacy() / kDisplacementName)); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "needs-reference-artifacts"); + } + { + ContractFixture fixture{"b31", mismatchedValues}; + auto input = readBytes(fixture.legacy() / kInputName); + replaceFirst(input, "type=B33", "type=B31"); + writeBytes(fixture.legacy() / kInputName, input); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "needs-reference-artifacts"); + } + { + ContractFixture fixture{"header", mismatchedValues}; + auto csv = readBytes(fixture.legacy() / kDisplacementName); + replaceFirst(csv, "U-U1", "U1"); + writeBytes(fixture.legacy() / kDisplacementName, csv); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "schema-mismatch"); + } + { + ContractFixture fixture{"missing-row", mismatchedValues}; + auto lines = readLines(fixture.legacy() / kDisplacementName); + ASSERT_EQ(lines.size(), kNodeCount + 1U); + lines.pop_back(); + writeLines(fixture.legacy() / kDisplacementName, lines); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "schema-mismatch"); + } + { + ContractFixture fixture{"extra-row", mismatchedValues}; + auto lines = readLines(fixture.legacy() / kDisplacementName); + ASSERT_EQ(lines.size(), kNodeCount + 1U); + std::string extra = lines.back(); + replaceFirst(extra, ",PART-1_1-1,11,", ",PART-1_1-1,12,"); + lines.push_back(std::move(extra)); + writeLines(fixture.legacy() / kDisplacementName, lines); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "schema-mismatch"); + } + { + ContractFixture fixture{"duplicate-row", mismatchedValues}; + auto lines = readLines(fixture.legacy() / kReactionName); + ASSERT_EQ(lines.size(), kNodeCount + 1U); + lines.push_back(lines[1U]); + writeLines(fixture.legacy() / kReactionName, lines); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "schema-mismatch"); + } + { + ContractFixture fixture{"nonfinite-row", mismatchedValues}; + auto csv = readBytes(fixture.legacy() / kReactionName); + replaceFirst(csv, "0.000000000E+00", "NaN"); + writeBytes(fixture.legacy() / kReactionName, csv); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "schema-mismatch"); + } + { + ContractFixture fixture{"identity", mismatchedValues}; + auto csv = readBytes(fixture.legacy() / kSectionName); + replaceFirst(csv, "PART-1_1-1", "WRONG-INSTANCE"); + writeBytes(fixture.legacy() / kSectionName, csv); + expectFailureCode( + fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()), + "schema-mismatch"); + } +} + +TEST(ReferenceComparisonContract, + AppliesAbaqusOnlyComponentScaleWithoutClampOrDrop) { + auto values = referenceValues(); + values.displacement[1U][0U] = 0.999e-9; + values.displacement[2U][0U] = 1.001e-9; + values.sectionResultants[0U][0U][2U] = 1.0e7 + 9.0; + values.sectionResultants[9U][1U][2U] = 0.0; + ContractFixture fixture{"tolerance", values}; + + auto result = fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()); + ASSERT_TRUE(result.hasValue()); + const auto& report = result.value(); + EXPECT_FALSE(report.passed); + EXPECT_EQ(report.rows.size(), kExpectedRowCount); + EXPECT_EQ(report.metrics.size(), kExpectedMetricCount); + + const auto* zero = findRow( + report, fesa::test::ComparisonQuantity::displacement, 1, "UX"); + const auto* nearZero = findRow( + report, fesa::test::ComparisonQuantity::displacement, 2, "UX"); + const auto* deliberateFailure = findRow( + report, fesa::test::ComparisonQuantity::displacement, 3, "UX"); + ASSERT_NE(zero, nullptr); + ASSERT_NE(nearZero, nullptr); + ASSERT_NE(deliberateFailure, nullptr); + EXPECT_DOUBLE_EQ(zero->reference.value, 0.0); + EXPECT_DOUBLE_EQ(zero->fesa.value, 0.0); + EXPECT_DOUBLE_EQ(zero->tolerance, 1.0e-9); + EXPECT_TRUE(zero->passed); + EXPECT_TRUE(nearZero->passed); + EXPECT_FALSE(deliberateFailure->passed); + + const auto* myMetric = findMetric( + report, fesa::test::ComparisonQuantity::sectionResultant, "My"); + ASSERT_NE(myMetric, nullptr); + EXPECT_DOUBLE_EQ(myMetric->referenceScale, 1.0e7); + const auto* scaled = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 1, "My"); + const auto* residue = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 11, "My"); + ASSERT_NE(scaled, nullptr); + ASSERT_NE(residue, nullptr); + EXPECT_DOUBLE_EQ(scaled->tolerance, 10.001); + EXPECT_DOUBLE_EQ(scaled->absoluteError, 9.0); + EXPECT_TRUE(scaled->passed); + EXPECT_DOUBLE_EQ(residue->reference.value, -1.56e-2); + EXPECT_DOUBLE_EQ(residue->fesa.value, 0.0); + EXPECT_DOUBLE_EQ(residue->absoluteError, 1.56e-2); + EXPECT_DOUBLE_EQ(residue->tolerance, 10.001); + EXPECT_TRUE(residue->passed); +} + +TEST(ReferenceComparisonContract, + ReportsEveryRowAndAggregateMetricDeterministically) { + auto values = referenceValues(); + values.displacement[0U][0U] = 0.5e-9; + values.displacement[1U][0U] = -1.0e-9; + ContractFixture fixture{"metrics", values}; + + auto result = fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()); + ASSERT_TRUE(result.hasValue()); + const auto& report = result.value(); + ASSERT_TRUE(report.passed); + expectExactRowInventory(report); + ASSERT_EQ(report.metrics.size(), kExpectedMetricCount); + EXPECT_TRUE(std::all_of( + report.rows.begin(), + report.rows.end(), + [](const fesa::test::RowDecision& row) { return row.passed; })); + + const auto* metric = findMetric( + report, fesa::test::ComparisonQuantity::displacement, "UX"); + const auto* worst = findRow( + report, fesa::test::ComparisonQuantity::displacement, 2, "UX"); + ASSERT_NE(metric, nullptr); + ASSERT_NE(worst, nullptr); + EXPECT_DOUBLE_EQ(metric->referenceScale, 0.0); + EXPECT_DOUBLE_EQ(metric->maximumAbsoluteError, 1.0e-9); + EXPECT_DOUBLE_EQ(metric->maximumNormalizedError, 1.0); + EXPECT_NEAR( + metric->rmsError, + std::sqrt(1.25 / static_cast(kNodeCount)) * 1.0e-9, + 1.0e-21); + EXPECT_NEAR(metric->normError, std::sqrt(1.25) * 1.0e-9, 1.0e-21); + EXPECT_EQ( + metric->worstRow, + static_cast(worst - report.rows.data())); + + EXPECT_FALSE(report.stressComparisonApplicable); + EXPECT_NE(report.stressComparisonReason.find("N/A"), std::string::npos); + EXPECT_NE(report.stressComparisonReason.find("HDF5"), std::string::npos); + EXPECT_DOUBLE_EQ(report.physicsEvidence.freeResidualNorm, 0.0); + EXPECT_EQ( + report.physicsEvidence.appliedForce, + (std::array{0.0, 0.0, -1.0e6})); + EXPECT_EQ( + report.physicsEvidence.reactionForce, + (std::array{0.0, 0.0, 1.0e6})); + EXPECT_EQ( + report.physicsEvidence.appliedMomentAboutOrigin, + (std::array{0.0, 1.0e7, 0.0})); + EXPECT_EQ( + report.physicsEvidence.reactionMomentAboutOrigin, + (std::array{0.0, -1.0e7, 0.0})); + EXPECT_TRUE(report.physicsEvidence.endpointConsistencyPassed); + + const auto jsonA = fixture.root() / "comparison-a.json"; + const auto jsonB = fixture.root() / "comparison-b.json"; + ASSERT_TRUE( + fesa::test::ReferenceComparison::writeDeterministicJson(report, jsonA) + .isOk()); + ASSERT_TRUE( + fesa::test::ReferenceComparison::writeDeterministicJson(report, jsonB) + .isOk()); + const std::string first = readBytes(jsonA); + EXPECT_EQ(first, readBytes(jsonB)); + for (const char* required : { + "\"rows\"", + "\"metrics\"", + "\"stress_comparison_applicable\":false", + "\"stress_comparison_reason\"", + "\"physics_evidence\"", + "\"free_residual_norm\"", + "\"applied_force\"", + "\"reaction_force\"", + "\"applied_moment_about_origin\"", + "\"reaction_moment_about_origin\"", + "\"endpoint_consistency_passed\""}) { + EXPECT_NE(first.find(required), std::string::npos) << required; + } +} + +TEST(ReferenceComparisonContract, + NormalizesEligibleStationsWithoutAveraging) { + auto values = referenceValues(); + values.sectionResultants[0U][0U] = {2.0e-4, 3.0e-4, 1.0e7, 4.0e-4}; + values.sectionResultants[0U][1U][2U] = 9.0e6 - 5.0; + values.sectionResultants[1U][0U][2U] = 9.0e6 + 5.0; + ContractFixture fixture{"stations", values}; + + auto result = fesa::test::ReferenceComparison::compare( + fixture.results(), fixture.legacy()); + ASSERT_TRUE(result.hasValue()); + const auto& report = result.value(); + ASSERT_TRUE(report.passed); + EXPECT_TRUE(report.physicsEvidence.endpointConsistencyPassed); + + const auto* n = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 1, "N"); + const auto* t = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 1, "T"); + const auto* my = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 1, "My"); + const auto* mz = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 1, "Mz"); + ASSERT_NE(n, nullptr); + ASSERT_NE(t, nullptr); + ASSERT_NE(my, nullptr); + ASSERT_NE(mz, nullptr); + EXPECT_DOUBLE_EQ(n->fesa.value, 2.0e-4); + EXPECT_DOUBLE_EQ(t->fesa.value, 3.0e-4); + EXPECT_DOUBLE_EQ(my->fesa.value, 1.0e7); + EXPECT_DOUBLE_EQ(mz->fesa.value, 4.0e-4); + + const auto* interior = findRow( + report, fesa::test::ComparisonQuantity::sectionResultant, 2, "My"); + ASSERT_NE(interior, nullptr); + EXPECT_DOUBLE_EQ(interior->fesa.value, 9.0e6 - 5.0); + EXPECT_DOUBLE_EQ(interior->reference.value, 9.0e6); + EXPECT_DOUBLE_EQ(interior->absoluteError, 5.0); + + auto mismatchValues = referenceValues(); + mismatchValues.sectionResultants[0U][1U][2U] = 9.0e6 - 6.0; + mismatchValues.sectionResultants[1U][0U][2U] = 9.0e6 + 6.0; + ContractFixture mismatch{"station-mismatch", mismatchValues}; + expectFailureCode( + fesa::test::ReferenceComparison::compare( + mismatch.results(), mismatch.legacy()), + "tolerance-failure"); +}