feat(linear-static-3d-euler-beam): step 22 - result-recovery
This commit is contained in:
KOKO\Mimi
@@ -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.
|
||||
|
||||
@@ -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 <array>
|
||||
#include <vector>
|
||||
|
||||
namespace fesa {
|
||||
|
||||
struct NodeStationResultRow {
|
||||
SourceEntityId node;
|
||||
EntityIndex representativeElement;
|
||||
std::array<double, 4> 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<std::vector<NodeStationResultRow>>
|
||||
normalizeSectionResultantsToNodeStations(
|
||||
const AnalysisModel& model,
|
||||
const std::vector<EndpointResultRow>& endpointRows,
|
||||
const std::array<double, 4>& componentTolerances);
|
||||
};
|
||||
|
||||
} // namespace fesa
|
||||
@@ -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
|
||||
)
|
||||
|
||||
|
||||
@@ -0,0 +1,740 @@
|
||||
#include "fesa/results/result_recovery.hpp"
|
||||
|
||||
#include "fesa/elements/euler_beam_3d.hpp"
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <charconv>
|
||||
#include <cmath>
|
||||
#include <cstddef>
|
||||
#include <cstdint>
|
||||
#include <limits>
|
||||
#include <optional>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <system_error>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
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<std::array<double, 3>, 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<class T>
|
||||
Result<T> recoveryResultFailure(const std::string& code,
|
||||
const SourceLocation& location,
|
||||
const std::string& identity,
|
||||
const std::string& message) {
|
||||
return Result<T>::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<double, 4>& values) {
|
||||
return std::all_of(
|
||||
values.begin(), values.end(),
|
||||
[](const double value) { return std::isfinite(value); });
|
||||
}
|
||||
|
||||
bool finite(const std::array<double, 6>& 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<std::size_t>& 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<std::size_t>& 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<std::size_t>::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<unsigned char> 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<std::size_t>(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<char>(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<std::vector<EntityIndex>> resolveLoadTarget(const Domain& domain,
|
||||
const NodalLoad& load) {
|
||||
std::vector<const NodeSet*> sets;
|
||||
for (const auto& set : domain.nodeSets()) {
|
||||
if (equalName(set.name, load.target)) {
|
||||
sets.push_back(&set);
|
||||
}
|
||||
}
|
||||
std::vector<EntityIndex> 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<EntityIndex>(node));
|
||||
}
|
||||
}
|
||||
}
|
||||
if (sets.size() > 1U || nodes.size() > 1U ||
|
||||
(!sets.empty() && !nodes.empty())) {
|
||||
return recoveryResultFailure<std::vector<EntityIndex>>(
|
||||
"invalid-node-station-entity",
|
||||
load.location,
|
||||
load.target,
|
||||
"A station-eligibility load target must resolve unambiguously.");
|
||||
}
|
||||
if (!sets.empty()) {
|
||||
std::vector<unsigned char> seen(domain.nodes().size(), 0U);
|
||||
for (const EntityIndex node : sets.front()->nodeIndices) {
|
||||
if (node >= domain.nodes().size() || seen[node] != 0U) {
|
||||
return recoveryResultFailure<std::vector<EntityIndex>>(
|
||||
"invalid-node-station-entity",
|
||||
load.location,
|
||||
load.target,
|
||||
"A station-eligibility node set must contain unique valid nodes.");
|
||||
}
|
||||
seen[node] = 1U;
|
||||
}
|
||||
return Result<std::vector<EntityIndex>>::success(
|
||||
sets.front()->nodeIndices);
|
||||
}
|
||||
if (!nodes.empty()) {
|
||||
return Result<std::vector<EntityIndex>>::success(std::move(nodes));
|
||||
}
|
||||
return recoveryResultFailure<std::vector<EntityIndex>>(
|
||||
"invalid-node-station-entity",
|
||||
load.location,
|
||||
load.target,
|
||||
"A station-eligibility load target must resolve to a node or node set.");
|
||||
}
|
||||
|
||||
std::array<double, 3> cross(const std::array<double, 3>& left,
|
||||
const std::array<double, 3>& 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<double, 3>& left,
|
||||
const std::array<double, 3>& right) {
|
||||
return left[0U] * right[0U] + left[1U] * right[1U] +
|
||||
left[2U] * right[2U];
|
||||
}
|
||||
|
||||
double norm(const std::array<double, 3>& value) {
|
||||
return std::hypot(value[0U], value[1U], value[2U]);
|
||||
}
|
||||
|
||||
std::optional<AxisSet> 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<double, 3> 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<double, 3> 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<double, 3> 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<double, 3> ey = {
|
||||
eyTrial[0U] / eyNorm,
|
||||
eyTrial[1U] / eyNorm,
|
||||
eyTrial[2U] / eyNorm};
|
||||
const std::array<double, 3> 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<double>::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<EndpointResultRow> endpointRows;
|
||||
std::vector<GaussResultRow> gaussRows;
|
||||
std::vector<StressS11Row> 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<int>(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<int>(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<std::vector<NodeStationResultRow>>
|
||||
ResultRecovery::normalizeSectionResultantsToNodeStations(
|
||||
const AnalysisModel& model,
|
||||
const std::vector<EndpointResultRow>& endpointRows,
|
||||
const std::array<double, 4>& componentTolerances) {
|
||||
const Domain& domain = model.domain();
|
||||
for (const double tolerance : componentTolerances) {
|
||||
if (!std::isfinite(tolerance) || tolerance < 0.0) {
|
||||
return recoveryResultFailure<std::vector<NodeStationResultRow>>(
|
||||
"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<std::size_t>::max)() / 2U ||
|
||||
endpointRows.size() != model.activeElements().size() * 2U) {
|
||||
return recoveryResultFailure<std::vector<NodeStationResultRow>>(
|
||||
"invalid-node-station-shape",
|
||||
{domain.sourcePath(), 0U},
|
||||
std::to_string(endpointRows.size()),
|
||||
"Endpoint rows must contain exactly two rows per active element.");
|
||||
}
|
||||
|
||||
std::vector<std::vector<const EndpointResultRow*>> 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<std::vector<NodeStationResultRow>>(
|
||||
"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<int>(endpoint) ||
|
||||
!sameSourceIdentity(row.node, domain.nodes()[nodeIndex].sourceId)) {
|
||||
return recoveryResultFailure<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>(
|
||||
"nonfinite-node-station-value",
|
||||
definition.location,
|
||||
definition.sourceId.sourceLabelText,
|
||||
"Node-station section resultants must be finite.");
|
||||
}
|
||||
rowsByNode[nodeIndex].push_back(&row);
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<unsigned char> loadedNodes(domain.nodes().size(), 0U);
|
||||
if (model.activeLoads().size() != model.step().loads.size()) {
|
||||
return recoveryResultFailure<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>(
|
||||
"nonfinite-node-station-value",
|
||||
load.location,
|
||||
load.target,
|
||||
"Station eligibility requires finite concentrated loads.");
|
||||
}
|
||||
auto targets = resolveLoadTarget(domain, load);
|
||||
if (!targets.hasValue()) {
|
||||
return Result<std::vector<NodeStationResultRow>>::failure(
|
||||
targets.status());
|
||||
}
|
||||
if (load.magnitude != 0.0) {
|
||||
for (const EntityIndex node : targets.value()) {
|
||||
loadedNodes[node] = 1U;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<NodeStationResultRow> 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<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>(
|
||||
"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<std::vector<NodeStationResultRow>>::success(
|
||||
std::move(stations));
|
||||
}
|
||||
|
||||
} // namespace fesa
|
||||
@@ -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
|
||||
)
|
||||
|
||||
@@ -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 <gtest/gtest.h>
|
||||
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <cstddef>
|
||||
#include <cstdint>
|
||||
#include <filesystem>
|
||||
#include <limits>
|
||||
#include <memory>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
namespace {
|
||||
|
||||
constexpr double kYoungsModulus = 100.0;
|
||||
constexpr double kPoissonRatio = 0.25;
|
||||
constexpr double kLength = 2.0;
|
||||
|
||||
struct RecoveryFixture {
|
||||
std::unique_ptr<fesa::Domain> domain;
|
||||
std::unique_ptr<fesa::AnalysisModel> model;
|
||||
std::unique_ptr<fesa::DofManager> dofs;
|
||||
std::unique_ptr<fesa::SparseMatrix> stiffness;
|
||||
};
|
||||
|
||||
fesa::ModelDefinition makeDefinition(
|
||||
const bool twoElements = false,
|
||||
std::vector<std::array<double, 2>> sectionPoints = {},
|
||||
std::vector<fesa::NodalLoad> 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<fesa::EntityIndex, 2>{2U, 1U}
|
||||
: std::array<fesa::EntityIndex, 2>{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<std::array<double, 2>> sectionPoints = {},
|
||||
std::vector<fesa::NodalLoad> 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<fesa::Domain>(
|
||||
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<fesa::AnalysisModel>(
|
||||
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<fesa::DofManager>(
|
||||
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<fesa::SparseMatrix>(
|
||||
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<fesa::EndpointResultRow> 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<double>::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<double, 4> expected = {
|
||||
kYoungsModulus * 2.0 * epsilon,
|
||||
shearModulus * 5.0 * twist,
|
||||
kYoungsModulus * 3.0 * kappaY,
|
||||
kYoungsModulus * 4.0 * kappaZ};
|
||||
const std::array<std::size_t, 4> 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<std::array<double, 2>> 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<int>(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<double, 4> 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<double>::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");
|
||||
}
|
||||
Reference in New Issue
Block a user