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 d674912..b7285e3 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 @@ -1021,3 +1021,112 @@ substitution, reconstruction, recovery, and output sequencing. - concerns: none; no critical implementation, environment, backend, or upstream-contract conflict was found. + +## Step 22 — result-recovery + +- task_id: `TASK-22` +- status: `completed` +- changed_files: `include/fesa/results/result_recovery.hpp`, + `src/fesa/results/result_recovery.cpp`, + `tests/unit/results/result_recovery_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-004`, `FESA-REQ-LS3DEB-027`, + `FESA-REQ-LS3DEB-031`, `FESA-REQ-LS3DEB-032`, + `FESA-REQ-LS3DEB-034`, `FESA-REQ-LS3DEB-035`, + `FESA-REQ-LS3DEB-042` +- test_ids: `T22-RECOVERY-001`, `T22-RECOVERY-002`, + `T22-RECOVERY-003`, `T22-RECOVERY-004`, `T22-RECOVERY-005`, + `T22-RECOVERY-006` + +| stage | exact command | exit_code | expected_or_observed_result | evidence_tail | +| --- | --- | ---: | --- | --- | +| RED-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 1 | Exactly six planned tests were registered before production and the recovery public API was absent | MSVC C1083 reported missing `fesa/results/result_recovery.hpp` from `result_recovery_test.cpp` | +| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | Recovery, station normalization, six tests, solver library, and unit executable compile and link | `result_recovery.cpp`, its test, `fesa_solver.lib`, and `fesa_unit_tests.exe` built without a FESA warning under `/W4 /WX` | +| GREEN-test | `ctest --test-dir .harness/build -C Debug -R ResultRecovery --output-on-failure` | 0 | Residual reaction, residual criterion, distinct result types, signs, stress order, and station policy pass | 6/6 exact `ResultRecovery` tests passed | +| RED-state-consistency | `cmake --build .harness/build --config Debug --target fesa_tests; ctest --test-dir .harness/build -C Debug -R ResultRecovery.EnforcesNormalizedFreeResidual --output-on-failure` | 1 | A constrained displacement inconsistent with prescribed `dc` must fail before residual evaluation | Test expected `invalid-recovery-state` but observed later `free-residual-tolerance-failure`, proving missing prevalidation | +| GREEN-state-consistency | same focused build/test command after the minimum validation change | 0 | Prescribed-order/state consistency fails before recovery and the exact zero-load/zero-displacement equilibrium remains valid | Focused test passed; subsequent exact suite passed 6/6 | +| 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 | Windows SDK 10.0.26100.0, oneMKL 2026.1 ILP64/dynamic, oneTBB, and HDF5 resolved; configure/generate completed | +| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build passes without a new FESA warning | `fesa_solver.lib` and `fesa_unit_tests.exe` built under `/W4 /WX` | +| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R ResultRecovery --output-on-failure` | 0 | Focused Step 22 suite remains green | 6/6 exact `ResultRecovery` tests passed | +| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | CTest discovers the accumulated suite and all six exact names | 68 tests discovered, including 6 `ResultRecovery` tests | +| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 68/68 tests passed | +| VERIFY-contract-scans | Exact-test-count, public-backend, forbidden shear/averaging, whitespace, and reference diff/status scans | 0 | Recovery scope and adapter/reference boundaries remain isolated | exact tests 6; public backend leaks 0; forbidden scope 0; whitespace clean; reference unchanged | + +- contract_checks: full `internalForce=K*d` and + `residual=internalForce-externalForce` are computed as private candidates; + the full residual is preserved in reaction so free residual evidence remains + visible and constrained entries are the physical support reactions. No + endpoint-action summation participates in reaction recovery. +- contract_checks: the free residual uses + `||r_f||/max(||Fint_f||,||Fext_f||)` with no unit floor. A positive physical + denominator must meet `1e-10`; exact zero numerator/denominator is accepted + as zero-load equilibrium, while nonfinite arithmetic fails closed. +- contract_checks: dimensions, finite input, stable free/constrained order, + exact prescribed displacement, active entity references, and twelve-DOF + scatters are checked before computing candidates. Internal/residual/reaction + vectors and endpoint/Gauss/stress rows replace AnalysisState only after all + elements recover finite results, so failure is atomic. +- contract_checks: active element order produces endpoint `0,1`, Gauss `1,2`, + and element/Gauss/section stress order. Equilibrium end action + `[FX,FY,FZ,MX,MY,MZ]`, endpoint `[N,T,My,Mz]`, Gauss strain/resultant, and + section-point/default-centroid `S11` remain distinct. +- contract_checks: node-station projection requires exact endpoint row + identity, a two-endpoint unloaded interior chain, identical section and + local axes, finite component tolerances `[N,T,My,Mz]`, and agreement before + selecting the smaller stable element ID. Positive-local-x section-cut rows + compare directly and are never averaged; reversed, branched, loaded, + section-jump, nonfinite, and mismatch cases fail structurally. +- generated_evidence: `.harness/build/src/fesa/Debug/fesa_solver.lib`, + `.harness/build/tests/Debug/fesa_unit_tests.exe` +- reference_diff: unchanged; `git diff --exit-code -- reference/` exit 0 +- handoff: Step 23 can serialize the stable nodal vectors and distinct + endpoint/Gauss/stress rows; Step 24 reference tooling can reuse + `normalizeSectionResultantsToNodeStations` without backend or HDF5 types. +- concerns: none; no critical implementation, environment, numerical, or + upstream-contract conflict was found. + +### Step 22 Review Fix Round 1 — full-space free-residual visibility + +- review_trigger: the approved I/O contract requires Reaction at every node + and component, including the accepted free-DOF residual; the AnalysisState + storage comment likewise makes reaction a full-index residual view. +- RED-test: `cmake --build .harness/build --config Debug --target fesa_tests; + ctest --test-dir .harness/build -C Debug -R + '^ResultRecovery.EnforcesNormalizedFreeResidual$' --output-on-failure` + built successfully and exited 1 in CTest. The accepted below-tolerance free + residual was `-1.000000082740371e-09`, while its reaction entry was `0`. +- production_fix: the private reaction candidate is now a full copy of + `residual=K*d-F`. Constrained entries remain the physical reactions, free + entries retain residual evidence, and element end actions are still never + summed into reaction. +- GREEN-focused: the same focused build/test command exited 0 and passed 1/1 + after the minimum production change. +- VERIFY-build: `cmake --build .harness/build --config Debug` exited 0 under + MSVC x64 Debug with no new FESA warning. +- VERIFY-targeted: `ctest --test-dir .harness/build -C Debug -R + ResultRecovery --output-on-failure` exited 0 and passed the unchanged exact + six names 6/6. +- VERIFY-discovery-full: `ctest --test-dir .harness/build -C Debug + --show-only=json-v1` discovered 68 tests, exactly six ResultRecovery tests; + full `ctest --test-dir .harness/build -C Debug --output-on-failure` exited 0 + and passed 68/68. +- failure_atomicity: the existing tolerance-failure branch in + `ResultRecovery.EnforcesNormalizedFreeResidual` still proves prior internal, + residual, reaction, and endpoint state remains unchanged on failure. +- supersession: this review section records the historical pre-fix finding + that implementation and report originally kept only constrained residual + entries. The final Step 22 contract above now states that reaction is the + full residual dataset; only constrained entries are physical reactions. +- phase_index: the review-fix worker intentionally left phase metadata + unchanged. After independent root verification, finalization preserved exact + `started_at=2026-08-09T20:42:15+0900`, added + `completed_at=2026-08-09T21:14:47+0900`, and updated the summary to state + the corrected full-residual reaction contract. +- evidence_hygiene: an initial parallel VERIFY/audit wrapper terminated before + product commands because its PowerShell count expression was malformed. It + is excluded from TDD and product VERIFY evidence; the recorded evidence is + the subsequent uncontended serial build, targeted, discovery, full sequence + followed by a separate successful audit command. +- concerns: none; no upstream-contract conflict or critical blocker remains. diff --git a/include/fesa/results/result_recovery.hpp b/include/fesa/results/result_recovery.hpp new file mode 100644 index 0000000..2344628 --- /dev/null +++ b/include/fesa/results/result_recovery.hpp @@ -0,0 +1,36 @@ +#pragma once + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/analysis/analysis_state.hpp" +#include "fesa/core/status.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/math/sparse_matrix.hpp" + +#include +#include + +namespace fesa { + +struct NodeStationResultRow { + SourceEntityId node; + EntityIndex representativeElement; + std::array sectionResultant; +}; + +// Recovers full-space equilibrium and the V0 beam output rows without +// exposing element or sparse-backend details to result consumers. +class ResultRecovery { +public: + static Status recover(const AnalysisModel& model, + const DofManager& dofs, + const SparseMatrix& fullStiffness, + AnalysisState& state); + + static Result> + normalizeSectionResultantsToNodeStations( + const AnalysisModel& model, + const std::vector& endpointRows, + const std::array& componentTolerances); +}; + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 0157d1a..0fe18f8 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -18,6 +18,7 @@ add_library( math/sparse_matrix.cpp math/vector.cpp model/domain.cpp + results/result_recovery.cpp solvers/linear/mkl_pardiso_solver.cpp ) diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp new file mode 100644 index 0000000..840c8f5 --- /dev/null +++ b/src/fesa/results/result_recovery.cpp @@ -0,0 +1,740 @@ +#include "fesa/results/result_recovery.hpp" + +#include "fesa/elements/euler_beam_3d.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +constexpr std::size_t kDofsPerNode = 6U; +constexpr std::size_t kElementDofCount = 12U; +constexpr double kFreeResidualTolerance = 1.0e-10; +constexpr double kAxisTolerance = 1.0e-12; + +using AxisSet = std::array, 3>; + +Status recoveryFailure(const std::string& code, + const SourceLocation& location, + const std::string& identity, + const std::string& message) { + return Status::failure( + FailureCategory::model, + {{Severity::error, + code, + location, + "RESULT_RECOVERY", + identity, + message}}); +} + +template +Result recoveryResultFailure(const std::string& code, + const SourceLocation& location, + const std::string& identity, + const std::string& message) { + return Result::failure( + recoveryFailure(code, location, identity, message)); +} + +bool sameSourceIdentity(const SourceEntityId& left, + const SourceEntityId& right) { + return left.instanceName == right.instanceName && + left.sourceLabel == right.sourceLabel && + left.sourceLabelText == right.sourceLabelText; +} + +bool finite(const std::array& values) { + return std::all_of( + values.begin(), values.end(), + [](const double value) { return std::isfinite(value); }); +} + +bool finite(const std::array& values) { + return std::all_of( + values.begin(), values.end(), + [](const double value) { return std::isfinite(value); }); +} + +bool finite(const Vector& values) { + for (std::size_t index = 0U; index < values.size(); ++index) { + if (!std::isfinite(values[index])) { + return false; + } + } + return true; +} + +double indexedNorm(const Vector& values, + const std::vector& indices) { + double result = 0.0; + for (const std::size_t index : indices) { + result = std::hypot(result, values[index]); + } + return result; +} + +bool strictlyIncreasing(const std::vector& values) { + return std::adjacent_find( + values.begin(), values.end(), + [](const std::size_t left, const std::size_t right) { + return left >= right; + }) == values.end(); +} + +Status validateRecoveryInputs(const AnalysisModel& model, + const DofManager& dofs, + const SparseMatrix& fullStiffness, + const AnalysisState& state) { + const Domain& domain = model.domain(); + if (domain.nodes().size() > + (std::numeric_limits::max)() / kDofsPerNode) { + return recoveryFailure( + "invalid-recovery-dimensions", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "The semantic node count cannot be represented in full-DOF space."); + } + const std::size_t fullCount = domain.nodes().size() * kDofsPerNode; + if (dofs.fullDofCount() != fullCount || + fullStiffness.rows() != fullCount || + fullStiffness.columns() != fullCount || + state.displacement().size() != fullCount || + state.externalForce().size() != fullCount || + state.internalForce().size() != fullCount || + state.residual().size() != fullCount || + state.reaction().size() != fullCount) { + return recoveryFailure( + "invalid-recovery-dimensions", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "Model, DOF, stiffness, and AnalysisState full-space dimensions must agree."); + } + const Status matrixStatus = fullStiffness.validate(); + if (!matrixStatus.isOk()) { + return matrixStatus; + } + if (!finite(state.displacement()) || !finite(state.externalForce())) { + return recoveryFailure( + "nonfinite-recovery-value", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "Displacement and external-force inputs must be finite."); + } + + const auto& freeDofs = dofs.freeDofs(); + const auto& constrainedDofs = dofs.constrainedDofs(); + if (freeDofs.size() != dofs.freeDofCount() || + constrainedDofs.size() != dofs.constrainedDofCount() || + dofs.prescribedValues().size() != constrainedDofs.size() || + freeDofs.size() + constrainedDofs.size() != fullCount || + !strictlyIncreasing(freeDofs) || + !strictlyIncreasing(constrainedDofs)) { + return recoveryFailure( + "invalid-recovery-order", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "Free and constrained DOFs must form stable increasing full-space orders."); + } + std::vector ownership(fullCount, 0U); + try { + for (std::size_t equation = 0U; equation < freeDofs.size(); ++equation) { + const std::size_t fullDof = freeDofs[equation]; + if (fullDof >= fullCount || ownership[fullDof] != 0U || + dofs.freeEquation(fullDof) != equation) { + return recoveryFailure( + "invalid-recovery-order", + {domain.sourcePath(), 0U}, + std::to_string(fullDof), + "Free equations must match stable full-DOF order."); + } + ownership[fullDof] = 1U; + } + for (std::size_t constrained = 0U; + constrained < constrainedDofs.size(); + ++constrained) { + const std::size_t fullDof = constrainedDofs[constrained]; + if (fullDof >= fullCount || ownership[fullDof] != 0U || + dofs.freeEquation(fullDof).has_value()) { + return recoveryFailure( + "invalid-recovery-order", + {domain.sourcePath(), 0U}, + std::to_string(fullDof), + "Constrained DOFs must be unique and absent from free equations."); + } + if (!std::isfinite(dofs.prescribedValues()[constrained]) || + state.displacement()[fullDof] != + dofs.prescribedValues()[constrained]) { + return recoveryFailure( + "invalid-recovery-state", + {domain.sourcePath(), 0U}, + std::to_string(fullDof), + "Constrained displacement must equal its prescribed value before recovery."); + } + ownership[fullDof] = 2U; + } + } catch (const std::out_of_range&) { + return recoveryFailure( + "invalid-recovery-dimensions", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "DofManager equation storage must cover every full DOF."); + } + if (std::find(ownership.begin(), ownership.end(), 0U) != ownership.end()) { + return recoveryFailure( + "invalid-recovery-order", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "Free and constrained DOFs must partition the full range."); + } + + EntityIndex previousElement = 0U; + bool firstElement = true; + for (const EntityIndex element : model.activeElements()) { + if (element >= domain.elements().size() || + (!firstElement && element <= previousElement)) { + return recoveryFailure( + "invalid-recovery-entity", + {domain.sourcePath(), 0U}, + std::to_string(element), + "Active elements must be unique in stable internal-index order."); + } + firstElement = false; + previousElement = element; + const auto& definition = domain.elements()[element]; + if (definition.nodeIndices[0U] >= domain.nodes().size() || + definition.nodeIndices[1U] >= domain.nodes().size() || + definition.materialIndex >= domain.materials().size() || + definition.sectionIndex >= domain.sections().size()) { + return recoveryFailure( + "invalid-recovery-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Active beam references must resolve before recovery."); + } + try { + const auto& scatter = dofs.elementScatter(element); + for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) { + for (std::size_t component = 0U; + component < kDofsPerNode; + ++component) { + const std::size_t expected = + static_cast(definition.nodeIndices[endpoint]) * + kDofsPerNode + + component; + if (scatter[endpoint * kDofsPerNode + component] != expected || + expected >= fullCount) { + return recoveryFailure( + "invalid-recovery-order", + definition.location, + definition.sourceId.sourceLabelText, + "Element scatter must preserve endpoint/component full-DOF order."); + } + } + } + } catch (const std::out_of_range&) { + return recoveryFailure( + "invalid-recovery-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Every active element requires one twelve-DOF scatter map."); + } + } + return Status::ok(); +} + +char asciiLower(const char value) { + if (value >= 'A' && value <= 'Z') { + return static_cast(value + ('a' - 'A')); + } + return value; +} + +bool equalName(const std::string& left, const std::string& right) { + return left.size() == right.size() && + std::equal( + left.begin(), left.end(), right.begin(), + [](const char leftValue, const char rightValue) { + return asciiLower(leftValue) == asciiLower(rightValue); + }); +} + +bool tryPositiveInteger(const std::string& text, std::int64_t& value) { + const char* const first = text.data(); + const char* const last = first + text.size(); + const auto parsed = std::from_chars(first, last, value); + return parsed.ec == std::errc{} && parsed.ptr == last && value > 0; +} + +Result> resolveLoadTarget(const Domain& domain, + const NodalLoad& load) { + std::vector sets; + for (const auto& set : domain.nodeSets()) { + if (equalName(set.name, load.target)) { + sets.push_back(&set); + } + } + std::vector nodes; + std::int64_t sourceLabel = 0; + if (tryPositiveInteger(load.target, sourceLabel)) { + for (std::size_t node = 0U; node < domain.nodes().size(); ++node) { + if (domain.nodes()[node].sourceId.sourceLabel == sourceLabel) { + nodes.push_back(static_cast(node)); + } + } + } + if (sets.size() > 1U || nodes.size() > 1U || + (!sets.empty() && !nodes.empty())) { + return recoveryResultFailure>( + "invalid-node-station-entity", + load.location, + load.target, + "A station-eligibility load target must resolve unambiguously."); + } + if (!sets.empty()) { + std::vector seen(domain.nodes().size(), 0U); + for (const EntityIndex node : sets.front()->nodeIndices) { + if (node >= domain.nodes().size() || seen[node] != 0U) { + return recoveryResultFailure>( + "invalid-node-station-entity", + load.location, + load.target, + "A station-eligibility node set must contain unique valid nodes."); + } + seen[node] = 1U; + } + return Result>::success( + sets.front()->nodeIndices); + } + if (!nodes.empty()) { + return Result>::success(std::move(nodes)); + } + return recoveryResultFailure>( + "invalid-node-station-entity", + load.location, + load.target, + "A station-eligibility load target must resolve to a node or node set."); +} + +std::array cross(const std::array& left, + const std::array& right) { + return { + left[1U] * right[2U] - left[2U] * right[1U], + left[2U] * right[0U] - left[0U] * right[2U], + left[0U] * right[1U] - left[1U] * right[0U]}; +} + +double dot(const std::array& left, + const std::array& right) { + return left[0U] * right[0U] + left[1U] * right[1U] + + left[2U] * right[2U]; +} + +double norm(const std::array& value) { + return std::hypot(value[0U], value[1U], value[2U]); +} + +std::optional localAxes(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 std::array delta = { + second[0U] - first[0U], + second[1U] - first[1U], + second[2U] - first[2U]}; + const double length = norm(delta); + if (!std::isfinite(length) || !(length > 0.0)) { + return std::nullopt; + } + const std::array ex = { + delta[0U] / length, delta[1U] / length, delta[2U] / length}; + const auto& guide = domain.sections()[element.sectionIndex].firstAxis; + const double projection = dot(guide, ex); + const std::array eyTrial = { + guide[0U] - projection * ex[0U], + guide[1U] - projection * ex[1U], + guide[2U] - projection * ex[2U]}; + const double eyNorm = norm(eyTrial); + if (!std::isfinite(eyNorm) || !(eyNorm > 0.0)) { + return std::nullopt; + } + const std::array ey = { + eyTrial[0U] / eyNorm, + eyTrial[1U] / eyNorm, + eyTrial[2U] / eyNorm}; + const std::array ez = cross(ex, ey); + const AxisSet axes = {ex, ey, ez}; + for (const auto& axis : axes) { + for (const double component : axis) { + if (!std::isfinite(component)) { + return std::nullopt; + } + } + } + return axes; +} + +bool sameAxes(const AxisSet& left, const AxisSet& right) { + for (std::size_t axis = 0U; axis < left.size(); ++axis) { + for (std::size_t component = 0U; + component < left[axis].size(); + ++component) { + if (std::abs(left[axis][component] - right[axis][component]) > + kAxisTolerance) { + return false; + } + } + } + return true; +} + +} // namespace + +Status ResultRecovery::recover(const AnalysisModel& model, + const DofManager& dofs, + const SparseMatrix& fullStiffness, + AnalysisState& state) { + const Status inputStatus = + validateRecoveryInputs(model, dofs, fullStiffness, state); + if (!inputStatus.isOk()) { + return inputStatus; + } + + Vector internalForce = fullStiffness.multiply(state.displacement()); + if (!finite(internalForce)) { + return recoveryFailure( + "nonfinite-recovery-value", + {model.domain().sourcePath(), 0U}, + model.domain().sourceContentIdentity(), + "Full stiffness multiplication must produce finite internal force."); + } + Vector residual{dofs.fullDofCount()}; + for (std::size_t fullDof = 0U; fullDof < residual.size(); ++fullDof) { + residual[fullDof] = + internalForce[fullDof] - state.externalForce()[fullDof]; + if (!std::isfinite(residual[fullDof])) { + return recoveryFailure( + "nonfinite-recovery-value", + {model.domain().sourcePath(), 0U}, + std::to_string(fullDof), + "Internal-minus-external residual must remain finite."); + } + } + + const double residualNorm = indexedNorm(residual, dofs.freeDofs()); + const double internalNorm = indexedNorm(internalForce, dofs.freeDofs()); + const double externalNorm = + indexedNorm(state.externalForce(), dofs.freeDofs()); + const double denominator = (std::max)(internalNorm, externalNorm); + if (!std::isfinite(residualNorm) || !std::isfinite(denominator)) { + return recoveryFailure( + "nonfinite-recovery-value", + {model.domain().sourcePath(), 0U}, + "free-residual", + "Free residual and its physical normalization scale must be finite."); + } + // An exact zero-load equilibrium is well-defined as zero. No unit floor is + // introduced; a nonzero residual with zero physical scale fails closed. + const double normalizedResidual = denominator == 0.0 + ? (residualNorm == 0.0 + ? 0.0 + : (std::numeric_limits::infinity)()) + : residualNorm / denominator; + if (!std::isfinite(normalizedResidual) || + normalizedResidual > kFreeResidualTolerance) { + return recoveryFailure( + "free-residual-tolerance-failure", + {model.domain().sourcePath(), 0U}, + "free-residual", + "The normalized free residual exceeds 1e-10."); + } + + // The full-space reaction dataset preserves free residual evidence while + // its constrained entries are the physical reactions from K*d-F. Element + // end actions remain distinct output and are never re-summed here. + Vector reaction = residual; + + std::vector endpointRows; + std::vector gaussRows; + std::vector stressRows; + endpointRows.reserve(model.activeElements().size() * 2U); + gaussRows.reserve(model.activeElements().size() * 2U); + const Domain& domain = model.domain(); + for (const EntityIndex elementIndex : model.activeElements()) { + const auto& definition = domain.elements()[elementIndex]; + auto beam = EulerBeam3D::create( + domain.nodes()[definition.nodeIndices[0U]], + domain.nodes()[definition.nodeIndices[1U]], + domain.sections()[definition.sectionIndex], + domain.materials()[definition.materialIndex]); + if (!beam.hasValue()) { + return beam.status(); + } + + Vector elementDisplacement{kElementDofCount}; + const auto& scatter = dofs.elementScatter(elementIndex); + for (std::size_t localDof = 0U; + localDof < kElementDofCount; + ++localDof) { + elementDisplacement[localDof] = + state.displacement()[scatter[localDof]]; + } + const BeamRecovery recovered = + beam.value().recover(elementDisplacement); + for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) { + if (!finite(recovered.equilibriumEndActions[endpoint]) || + !finite(recovered.endpointSectionResultants[endpoint])) { + return recoveryFailure( + "nonfinite-recovery-value", + definition.location, + definition.sourceId.sourceLabelText, + "Endpoint recovery values must be finite."); + } + endpointRows.push_back({ + elementIndex, + static_cast(endpoint), + domain.nodes()[definition.nodeIndices[endpoint]].sourceId, + recovered.equilibriumEndActions[endpoint], + recovered.endpointSectionResultants[endpoint]}); + } + for (std::size_t gauss = 0U; gauss < 2U; ++gauss) { + if (!finite(recovered.gaussGeneralizedStrains[gauss]) || + !finite(recovered.gaussGeneralizedResultants[gauss])) { + return recoveryFailure( + "nonfinite-recovery-value", + definition.location, + definition.sourceId.sourceLabelText, + "Gauss recovery values must be finite."); + } + gaussRows.push_back({ + elementIndex, + static_cast(gauss + 1U), + recovered.gaussGeneralizedStrains[gauss], + recovered.gaussGeneralizedResultants[gauss]}); + } + for (const auto& point : recovered.stressPoints) { + if ((point.gaussPoint != 1 && point.gaussPoint != 2) || + !std::isfinite(point.x1) || !std::isfinite(point.x2) || + !std::isfinite(point.s11)) { + return recoveryFailure( + "nonfinite-recovery-value", + definition.location, + definition.sourceId.sourceLabelText, + "Stress recovery identity and values must be finite and ordered."); + } + stressRows.push_back({ + elementIndex, + point.gaussPoint, + point.sectionPoint, + point.x1, + point.x2, + point.s11, + point.source}); + } + } + + // Commit only after all validation and element recovery succeeds so a + // failed recovery cannot leave a partially updated AnalysisState. + state.internalForce() = std::move(internalForce); + state.residual() = std::move(residual); + state.reaction() = std::move(reaction); + state.endpointResults() = std::move(endpointRows); + state.gaussResults() = std::move(gaussRows); + state.stressResults() = std::move(stressRows); + return Status::ok(); +} + +Result> +ResultRecovery::normalizeSectionResultantsToNodeStations( + const AnalysisModel& model, + const std::vector& endpointRows, + const std::array& componentTolerances) { + const Domain& domain = model.domain(); + for (const double tolerance : componentTolerances) { + if (!std::isfinite(tolerance) || tolerance < 0.0) { + return recoveryResultFailure>( + "invalid-node-station-tolerance", + {domain.sourcePath(), 0U}, + "component-tolerances", + "Node-station component tolerances must be finite and nonnegative."); + } + } + if (model.activeElements().size() > + (std::numeric_limits::max)() / 2U || + endpointRows.size() != model.activeElements().size() * 2U) { + return recoveryResultFailure>( + "invalid-node-station-shape", + {domain.sourcePath(), 0U}, + std::to_string(endpointRows.size()), + "Endpoint rows must contain exactly two rows per active element."); + } + + std::vector> rowsByNode( + domain.nodes().size()); + for (std::size_t order = 0U; + order < model.activeElements().size(); + ++order) { + const EntityIndex elementIndex = model.activeElements()[order]; + if (elementIndex >= domain.elements().size()) { + return recoveryResultFailure>( + "invalid-node-station-entity", + {domain.sourcePath(), 0U}, + std::to_string(elementIndex), + "Every active station element must be a valid stable entity."); + } + const auto& definition = domain.elements()[elementIndex]; + for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) { + const auto& row = endpointRows[order * 2U + endpoint]; + const EntityIndex nodeIndex = definition.nodeIndices[endpoint]; + if (nodeIndex >= domain.nodes().size() || + row.element != elementIndex || + row.endpoint != static_cast(endpoint) || + !sameSourceIdentity(row.node, domain.nodes()[nodeIndex].sourceId)) { + return recoveryResultFailure>( + "invalid-node-station-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Endpoint rows must preserve active element, endpoint, and source-node order."); + } + if (!finite(row.sectionResultant)) { + return recoveryResultFailure>( + "nonfinite-node-station-value", + definition.location, + definition.sourceId.sourceLabelText, + "Node-station section resultants must be finite."); + } + rowsByNode[nodeIndex].push_back(&row); + } + } + + std::vector loadedNodes(domain.nodes().size(), 0U); + if (model.activeLoads().size() != model.step().loads.size()) { + return recoveryResultFailure>( + "invalid-node-station-entity", + model.step().location, + model.step().name, + "The active load view must preserve every sole-step load."); + } + for (std::size_t order = 0U; order < model.activeLoads().size(); ++order) { + const EntityIndex loadIndex = model.activeLoads()[order]; + if (loadIndex != order || loadIndex >= model.step().loads.size()) { + return recoveryResultFailure>( + "invalid-node-station-entity", + model.step().location, + std::to_string(loadIndex), + "Active loads must remain in stable source order."); + } + const auto& load = model.step().loads[loadIndex]; + if (!std::isfinite(load.magnitude)) { + return recoveryResultFailure>( + "nonfinite-node-station-value", + load.location, + load.target, + "Station eligibility requires finite concentrated loads."); + } + auto targets = resolveLoadTarget(domain, load); + if (!targets.hasValue()) { + return Result>::failure( + targets.status()); + } + if (load.magnitude != 0.0) { + for (const EntityIndex node : targets.value()) { + loadedNodes[node] = 1U; + } + } + } + + std::vector stations; + stations.reserve(domain.nodes().size()); + for (std::size_t nodeIndex = 0U; + nodeIndex < rowsByNode.size(); + ++nodeIndex) { + const auto& incident = rowsByNode[nodeIndex]; + if (incident.empty()) { + continue; + } + if (incident.size() == 1U) { + stations.push_back({ + domain.nodes()[nodeIndex].sourceId, + incident.front()->element, + incident.front()->sectionResultant}); + continue; + } + if (incident.size() != 2U || loadedNodes[nodeIndex] != 0U) { + return recoveryResultFailure>( + "ineligible-node-station", + domain.nodes()[nodeIndex].location, + domain.nodes()[nodeIndex].sourceId.sourceLabelText, + "Interior station collapse requires exactly two unloaded endpoints."); + } + + const auto& firstElement = domain.elements()[incident[0U]->element]; + const auto& secondElement = domain.elements()[incident[1U]->element]; + const bool chainOrientation = + incident[0U]->endpoint != incident[1U]->endpoint && + ((incident[0U]->endpoint == 1 && incident[1U]->endpoint == 0) || + (incident[0U]->endpoint == 0 && incident[1U]->endpoint == 1)); + const auto firstAxes = localAxes(domain, firstElement); + const auto secondAxes = localAxes(domain, secondElement); + if (!chainOrientation || + firstElement.sectionIndex != secondElement.sectionIndex || + !firstAxes.has_value() || !secondAxes.has_value() || + !sameAxes(*firstAxes, *secondAxes)) { + return recoveryResultFailure>( + "ineligible-node-station", + domain.nodes()[nodeIndex].location, + domain.nodes()[nodeIndex].sourceId.sourceLabelText, + "Interior station endpoints require one consistent section and local-axis chain."); + } + + // Endpoint sectionResultant rows already use the positive-local-x cut + // convention. Once common orientation is proven, no outward-action + // endpoint sign is applied and the values are directly comparable. + for (std::size_t component = 0U; + component < componentTolerances.size(); + ++component) { + const double difference = std::abs( + incident[0U]->sectionResultant[component] - + incident[1U]->sectionResultant[component]); + if (!std::isfinite(difference)) { + return recoveryResultFailure>( + "nonfinite-node-station-value", + domain.nodes()[nodeIndex].location, + domain.nodes()[nodeIndex].sourceId.sourceLabelText, + "Endpoint comparison must produce a finite difference."); + } + if (difference > componentTolerances[component]) { + return recoveryResultFailure>( + "node-station-tolerance-failure", + domain.nodes()[nodeIndex].location, + domain.nodes()[nodeIndex].sourceId.sourceLabelText, + "Interior endpoint resultants disagree beyond component tolerance."); + } + } + + const EndpointResultRow* representative = + incident[0U]->element < incident[1U]->element + ? incident[0U] + : incident[1U]; + stations.push_back({ + domain.nodes()[nodeIndex].sourceId, + representative->element, + representative->sectionResultant}); + } + return Result>::success( + std::move(stations)); +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7a93343..0e8b38b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -23,6 +23,7 @@ add_executable( unit/model/domain_test.cpp unit/model/model_types_test.cpp unit/results/result_records_test.cpp + unit/results/result_recovery_test.cpp unit/solvers/linear/linear_solver_test.cpp unit/solvers/linear/mkl_pardiso_solver_test.cpp ) diff --git a/tests/unit/results/result_recovery_test.cpp b/tests/unit/results/result_recovery_test.cpp new file mode 100644 index 0000000..c7ced70 --- /dev/null +++ b/tests/unit/results/result_recovery_test.cpp @@ -0,0 +1,471 @@ +#include "fesa/results/result_recovery.hpp" + +#include "fesa/analysis/analysis_model.hpp" +#include "fesa/analysis/analysis_state.hpp" +#include "fesa/assembly/parallel_for.hpp" +#include "fesa/assembly/sparse_assembler.hpp" +#include "fesa/fem/dof_manager.hpp" +#include "fesa/model/domain.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr double kYoungsModulus = 100.0; +constexpr double kPoissonRatio = 0.25; +constexpr double kLength = 2.0; + +struct RecoveryFixture { + std::unique_ptr domain; + std::unique_ptr model; + std::unique_ptr dofs; + std::unique_ptr stiffness; +}; + +fesa::ModelDefinition makeDefinition( + const bool twoElements = false, + std::vector> sectionPoints = {}, + std::vector loads = {}, + const bool reverseSecond = false, + const bool sectionJump = false, + const bool nonzeroPrescription = true) { + const std::filesystem::path source{"models/result-recovery.inp"}; + fesa::ModelDefinition definition{}; + definition.sourcePath = source; + definition.sourceContentIdentity = "fnv1a64:0123456789abcdef"; + definition.nodes = { + {{"Beam-1", 1, "1"}, {0.0, 0.0, 0.0}, {source, 10U}}, + {{"Beam-1", 2, "2"}, {kLength, 0.0, 0.0}, {source, 11U}}}; + if (twoElements) { + definition.nodes.push_back( + {{"Beam-1", 3, "3"}, {2.0 * kLength, 0.0, 0.0}, {source, 12U}}); + } + definition.materials = { + {"Material", kYoungsModulus, kPoissonRatio, {source, 20U}}}; + definition.sections = {{ + "Section", + 2.0, + 3.0, + 0.0, + 4.0, + 5.0, + {0.0, 1.0, 0.0}, + std::move(sectionPoints), + {source, 30U}}}; + if (sectionJump) { + auto secondSection = definition.sections.front(); + secondSection.name = "Section-2"; + secondSection.area = 2.5; + secondSection.location.line = 31U; + definition.sections.push_back(std::move(secondSection)); + } + definition.elements = { + {{"Beam-1", 10, "10"}, {0U, 1U}, 0U, 0U, {source, 40U}}}; + if (twoElements) { + definition.elements.push_back({ + {"Beam-1", 20, "20"}, + reverseSecond ? std::array{2U, 1U} + : std::array{1U, 2U}, + 0U, + sectionJump ? 1U : 0U, + {source, 41U}}); + } + definition.steps = {{ + "Step-1", + {{"1", 1, 1, nonzeroPrescription ? 0.1 : 0.0, {source, 50U}}, + {"1", 2, 6, 0.0, {source, 51U}}}, + std::move(loads), + 0.1, + 1.0, + 0.01, + 1.0, + {source, 49U}}}; + return definition; +} + +RecoveryFixture makeFixture( + const bool twoElements = false, + std::vector> sectionPoints = {}, + std::vector loads = {}, + const bool reverseSecond = false, + const bool sectionJump = false, + const bool nonzeroPrescription = true) { + auto domainResult = fesa::Domain::create(makeDefinition( + twoElements, + std::move(sectionPoints), + std::move(loads), + reverseSecond, + sectionJump, + nonzeroPrescription)); + if (!domainResult.hasValue()) { + throw std::runtime_error{"Recovery fixture Domain construction failed."}; + } + auto domain = std::make_unique( + std::move(domainResult.value())); + + auto modelResult = fesa::AnalysisModel::create(*domain); + if (!modelResult.hasValue()) { + throw std::runtime_error{"Recovery fixture AnalysisModel construction failed."}; + } + auto model = std::make_unique( + std::move(modelResult.value())); + + auto dofsResult = fesa::DofManager::create(*model); + if (!dofsResult.hasValue()) { + throw std::runtime_error{"Recovery fixture DofManager construction failed."}; + } + auto dofs = std::make_unique( + std::move(dofsResult.value())); + + fesa::SerialParallelFor serial; + auto stiffnessResult = fesa::SparseAssembler::assembleStiffness( + *model, *dofs, serial); + if (!stiffnessResult.hasValue()) { + throw std::runtime_error{"Recovery fixture stiffness assembly failed."}; + } + auto stiffness = std::make_unique( + std::move(stiffnessResult.value())); + return { + std::move(domain), + std::move(model), + std::move(dofs), + std::move(stiffness)}; +} + +fesa::AnalysisState makeAxialEquilibriumState(const RecoveryFixture& fixture) { + auto state = fesa::AnalysisState::create( + *fixture.dofs, {"Step-1", 0U}); + state.displacement()[0U] = 0.1; + state.displacement()[6U] = 0.3; + const fesa::Vector internal = fixture.stiffness->multiply(state.displacement()); + for (const std::size_t fullDof : fixture.dofs->freeDofs()) { + state.externalForce()[fullDof] = internal[fullDof]; + } + return state; +} + +fesa::AnalysisState makePatchState( + const RecoveryFixture& fixture, + const double epsilon, + const double twist, + const double kappaY, + const double kappaZ) { + auto state = fesa::AnalysisState::create( + *fixture.dofs, {"Step-1", 0U}); + state.displacement()[0U] = 0.1; + state.displacement()[6U] = 0.1 + epsilon * kLength; + state.displacement()[7U] = 0.5 * kappaZ * kLength * kLength; + state.displacement()[8U] = -0.5 * kappaY * kLength * kLength; + state.displacement()[9U] = twist * kLength; + state.displacement()[10U] = kappaY * kLength; + state.displacement()[11U] = kappaZ * kLength; + state.externalForce() = fixture.stiffness->multiply(state.displacement()); + return state; +} + +void expectStatusCode(const fesa::Status& status, const std::string& code) { + ASSERT_FALSE(status.isOk()); + EXPECT_EQ(status.failureCategory(), fesa::FailureCategory::model); + ASSERT_EQ(status.diagnostics().size(), 1U); + EXPECT_EQ(status.diagnostics()[0U].code, code); +} + +void expectScaledNear( + const double actual, + const double expected, + const double relativeTolerance = 1.0e-12) { + ASSERT_TRUE(std::isfinite(actual)); + ASSERT_TRUE(std::isfinite(expected)); + EXPECT_LE( + std::abs(actual - expected), + relativeTolerance * (std::max)(std::abs(expected), 1.0)); +} + +std::vector makeStationRows( + const RecoveryFixture& fixture) { + const auto& nodes = fixture.domain->nodes(); + const auto& elements = fixture.domain->elements(); + return { + {0U, 0, nodes[elements[0U].nodeIndices[0U]].sourceId, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {1.0, 2.0, 3.0, 4.0}}, + {0U, 1, nodes[elements[0U].nodeIndices[1U]].sourceId, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {5.0, 6.0, 7.0, 8.0}}, + {1U, 0, nodes[elements[1U].nodeIndices[0U]].sourceId, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {5.0, 6.0, 7.0, 8.0}}, + {1U, 1, nodes[elements[1U].nodeIndices[1U]].sourceId, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {9.0, 10.0, 11.0, 12.0}}}; +} + +} // namespace + +TEST(ResultRecovery, ComputesResidualReactionForNonzeroPrescription) { + const auto fixture = makeFixture(); + auto state = makeAxialEquilibriumState(fixture); + + const fesa::Status status = fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state); + + ASSERT_TRUE(status.isOk()); + EXPECT_DOUBLE_EQ(state.internalForce()[0U], -20.0); + EXPECT_DOUBLE_EQ(state.internalForce()[6U], 20.0); + EXPECT_DOUBLE_EQ(state.residual()[0U], -20.0); + EXPECT_DOUBLE_EQ(state.residual()[6U], 0.0); + EXPECT_DOUBLE_EQ(state.reaction()[0U], -20.0); + for (const std::size_t fullDof : fixture.dofs->freeDofs()) { + EXPECT_DOUBLE_EQ(state.reaction()[fullDof], 0.0); + } +} + +TEST(ResultRecovery, EnforcesNormalizedFreeResidual) { + const auto fixture = makeFixture(); + + auto failed = makeAxialEquilibriumState(fixture); + failed.internalForce()[0U] = 91.0; + failed.residual()[0U] = 92.0; + failed.reaction()[0U] = 93.0; + failed.reaction()[6U] = 94.0; + failed.endpointResults().push_back({}); + failed.externalForce()[6U] += 1.0e-7; + expectStatusCode( + fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, failed), + "free-residual-tolerance-failure"); + EXPECT_DOUBLE_EQ(failed.internalForce()[0U], 91.0); + EXPECT_DOUBLE_EQ(failed.residual()[0U], 92.0); + EXPECT_DOUBLE_EQ(failed.reaction()[0U], 93.0); + EXPECT_DOUBLE_EQ(failed.reaction()[6U], 94.0); + EXPECT_EQ(failed.endpointResults().size(), 1U); + + auto thresholdPass = makeAxialEquilibriumState(fixture); + thresholdPass.externalForce()[6U] += 1.0e-10 * 20.0 * 0.5; + const fesa::Status thresholdStatus = fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + thresholdPass); + ASSERT_TRUE(thresholdStatus.isOk()); + EXPECT_NE(thresholdPass.residual()[6U], 0.0); + EXPECT_DOUBLE_EQ( + thresholdPass.reaction()[6U], thresholdPass.residual()[6U]); + + const auto zeroFixture = makeFixture(false, {}, {}, false, false, false); + auto zeroEquilibrium = fesa::AnalysisState::create( + *zeroFixture.dofs, {"Step-1", 0U}); + EXPECT_TRUE(fesa::ResultRecovery::recover( + *zeroFixture.model, + *zeroFixture.dofs, + *zeroFixture.stiffness, + zeroEquilibrium) + .isOk()); + + auto wrongPrescription = makeAxialEquilibriumState(fixture); + wrongPrescription.displacement()[0U] = 0.0; + expectStatusCode( + fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + wrongPrescription), + "invalid-recovery-state"); + + auto nonfinite = makeAxialEquilibriumState(fixture); + nonfinite.displacement()[6U] = std::numeric_limits::quiet_NaN(); + expectStatusCode( + fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, nonfinite), + "nonfinite-recovery-value"); + + const auto wrongFixture = makeFixture(true); + auto wrongState = fesa::AnalysisState::create( + *wrongFixture.dofs, {"Step-1", 0U}); + expectStatusCode( + fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, wrongState), + "invalid-recovery-dimensions"); +} + +TEST(ResultRecovery, KeepsEndActionSectionAndGaussResultsDistinct) { + const auto fixture = makeFixture(); + auto state = makePatchState(fixture, 0.02, 0.03, -0.04, 0.05); + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state) + .isOk()); + ASSERT_EQ(state.endpointResults().size(), 2U); + ASSERT_EQ(state.gaussResults().size(), 2U); + EXPECT_EQ(state.endpointResults()[0U].endpoint, 0); + EXPECT_EQ(state.endpointResults()[1U].endpoint, 1); + EXPECT_EQ(state.gaussResults()[0U].gaussPoint, 1); + EXPECT_EQ(state.gaussResults()[1U].gaussPoint, 2); + EXPECT_DOUBLE_EQ( + state.endpointResults()[0U].endAction[0U], + -state.endpointResults()[0U].sectionResultant[0U]); + EXPECT_DOUBLE_EQ( + state.endpointResults()[1U].endAction[0U], + state.endpointResults()[1U].sectionResultant[0U]); + EXPECT_DOUBLE_EQ( + state.gaussResults()[0U].generalizedResultant[0U], + state.endpointResults()[0U].sectionResultant[0U]); +} + +TEST(ResultRecovery, MatchesAxialTorsionAndTwoPlaneEndSigns) { + const auto fixture = makeFixture(); + const double epsilon = 0.02; + const double twist = -0.03; + const double kappaY = 0.04; + const double kappaZ = -0.05; + auto state = makePatchState(fixture, epsilon, twist, kappaY, kappaZ); + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state) + .isOk()); + const double shearModulus = + kYoungsModulus / (2.0 * (1.0 + kPoissonRatio)); + const std::array expected = { + kYoungsModulus * 2.0 * epsilon, + shearModulus * 5.0 * twist, + kYoungsModulus * 3.0 * kappaY, + kYoungsModulus * 4.0 * kappaZ}; + const std::array endComponents = {0U, 3U, 4U, 5U}; + for (std::size_t component = 0U; component < expected.size(); ++component) { + expectScaledNear( + state.endpointResults()[0U].sectionResultant[component], + expected[component]); + expectScaledNear( + state.endpointResults()[1U].sectionResultant[component], + expected[component]); + expectScaledNear( + state.endpointResults()[0U].endAction[endComponents[component]], + -expected[component]); + expectScaledNear( + state.endpointResults()[1U].endAction[endComponents[component]], + expected[component]); + } +} + +TEST(ResultRecovery, OrdersStressPointsAndDefaultCentroid) { + const std::vector> sectionPoints = { + {0.25, -0.5}, {-0.4, 0.3}}; + const auto fixture = makeFixture(false, sectionPoints); + auto state = makePatchState(fixture, 0.01, 0.0, 0.02, -0.03); + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state) + .isOk()); + ASSERT_EQ(state.stressResults().size(), 4U); + for (std::size_t gauss = 0U; gauss < 2U; ++gauss) { + for (std::size_t point = 0U; point < sectionPoints.size(); ++point) { + const auto& row = state.stressResults()[gauss * 2U + point]; + EXPECT_EQ(row.element, 0U); + EXPECT_EQ(row.gaussPoint, static_cast(gauss + 1U)); + EXPECT_EQ(row.sectionPoint, point + 1U); + EXPECT_DOUBLE_EQ(row.x1, sectionPoints[point][0U]); + EXPECT_DOUBLE_EQ(row.x2, sectionPoints[point][1U]); + EXPECT_EQ(row.source, "input"); + expectScaledNear( + row.s11, + kYoungsModulus * + (0.01 + row.x2 * 0.02 - row.x1 * -0.03)); + } + } + + const auto defaultFixture = makeFixture(); + auto defaultState = makePatchState(defaultFixture, 0.01, 0.0, 0.0, 0.0); + ASSERT_TRUE(fesa::ResultRecovery::recover( + *defaultFixture.model, + *defaultFixture.dofs, + *defaultFixture.stiffness, + defaultState) + .isOk()); + ASSERT_EQ(defaultState.stressResults().size(), 2U); + for (const auto& row : defaultState.stressResults()) { + EXPECT_EQ(row.sectionPoint, 0U); + EXPECT_DOUBLE_EQ(row.x1, 0.0); + EXPECT_DOUBLE_EQ(row.x2, 0.0); + EXPECT_EQ(row.source, "fesa-default"); + } +} + +TEST(ResultRecovery, RequiresInteriorEndpointConsistencyWithoutAveraging) { + const auto fixture = makeFixture(true); + const std::array tolerances = {1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6}; + auto rows = makeStationRows(fixture); + rows[2U].sectionResultant[0U] += 0.5e-6; + + auto normalized = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *fixture.model, rows, tolerances); + ASSERT_TRUE(normalized.hasValue()); + ASSERT_EQ(normalized.value().size(), 3U); + EXPECT_EQ(normalized.value()[1U].representativeElement, 0U); + EXPECT_DOUBLE_EQ(normalized.value()[1U].sectionResultant[0U], 5.0); + + rows[2U].sectionResultant[0U] = 5.0 + 2.0e-6; + auto mismatch = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *fixture.model, rows, tolerances); + ASSERT_FALSE(mismatch.hasValue()); + expectStatusCode(mismatch.status(), "node-station-tolerance-failure"); + + rows = makeStationRows(fixture); + rows[2U].sectionResultant[1U] = + std::numeric_limits::infinity(); + auto nonfinite = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *fixture.model, rows, tolerances); + ASSERT_FALSE(nonfinite.hasValue()); + expectStatusCode(nonfinite.status(), "nonfinite-node-station-value"); + + auto invalidTolerance = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *fixture.model, + makeStationRows(fixture), + {1.0e-6, -1.0, 1.0e-6, 1.0e-6}); + ASSERT_FALSE(invalidTolerance.hasValue()); + expectStatusCode( + invalidTolerance.status(), "invalid-node-station-tolerance"); + + const std::filesystem::path source{"models/result-recovery.inp"}; + const auto loadedFixture = makeFixture( + true, {}, {{"2", 2, 1.0, {source, 60U}}}); + auto loaded = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *loadedFixture.model, + makeStationRows(loadedFixture), + tolerances); + ASSERT_FALSE(loaded.hasValue()); + expectStatusCode(loaded.status(), "ineligible-node-station"); + + const auto reversedFixture = makeFixture(true, {}, {}, true); + auto reversed = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *reversedFixture.model, + makeStationRows(reversedFixture), + tolerances); + ASSERT_FALSE(reversed.hasValue()); + expectStatusCode(reversed.status(), "ineligible-node-station"); + + const auto jumpFixture = makeFixture(true, {}, {}, false, true); + auto jumped = + fesa::ResultRecovery::normalizeSectionResultantsToNodeStations( + *jumpFixture.model, + makeStationRows(jumpFixture), + tolerances); + ASSERT_FALSE(jumped.hasValue()); + expectStatusCode(jumped.status(), "ineligible-node-station"); +}