From f60d4edc261e862febe739b1dc31471360931fe1 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 16 Aug 2026 12:52:11 +0900 Subject: [PATCH] feat(cpp-object-oriented-modular-refactoring): step 22 - result-recovery-modules --- src/fesa/CMakeLists.txt | 5 + src/fesa/results/analysis_state_commit.cpp | 86 +++ src/fesa/results/analysis_state_commit.h | 16 + src/fesa/results/beam_result_recovery.cpp | 106 ++++ src/fesa/results/beam_result_recovery.h | 23 + .../results/global_equilibrium_recovery.cpp | 269 ++++++++ .../results/global_equilibrium_recovery.h | 37 ++ src/fesa/results/recovery_candidate.cpp | 39 ++ src/fesa/results/recovery_candidate.h | 65 ++ src/fesa/results/result_recovery.cpp | 417 ++---------- src/fesa/results/shell_result_recovery.cpp | 65 ++ src/fesa/results/shell_result_recovery.h | 20 + tests/CMakeLists.txt | 1 + .../result_recovery_components_test.cpp | 593 ++++++++++++++++++ 14 files changed, 1373 insertions(+), 369 deletions(-) create mode 100644 src/fesa/results/analysis_state_commit.cpp create mode 100644 src/fesa/results/analysis_state_commit.h create mode 100644 src/fesa/results/beam_result_recovery.cpp create mode 100644 src/fesa/results/beam_result_recovery.h create mode 100644 src/fesa/results/global_equilibrium_recovery.cpp create mode 100644 src/fesa/results/global_equilibrium_recovery.h create mode 100644 src/fesa/results/recovery_candidate.cpp create mode 100644 src/fesa/results/recovery_candidate.h create mode 100644 src/fesa/results/shell_result_recovery.cpp create mode 100644 src/fesa/results/shell_result_recovery.h create mode 100644 tests/unit/results/result_recovery_components_test.cpp diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index b5c799a..cbb43b6 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -37,7 +37,12 @@ add_library( model/source_target_resolver.cpp properties/general_beam_section.cpp properties/shell_section.cpp + results/analysis_state_commit.cpp + results/beam_result_recovery.cpp + results/global_equilibrium_recovery.cpp + results/recovery_candidate.cpp results/result_recovery.cpp + results/shell_result_recovery.cpp solvers/linear/mkl_pardiso_solver.cpp ) diff --git a/src/fesa/results/analysis_state_commit.cpp b/src/fesa/results/analysis_state_commit.cpp new file mode 100644 index 0000000..df6ed96 --- /dev/null +++ b/src/fesa/results/analysis_state_commit.cpp @@ -0,0 +1,86 @@ +#include "results/analysis_state_commit.h" + +#include +#include + +namespace fesa::results_internal { + +Status CommitRecoveryCandidate(RecoveryCandidate candidate, + AnalysisState& state) { + const std::size_t full_count = state.Displacement().Size(); + if (candidate.displacement.Size() != full_count || + candidate.external_force.Size() != full_count || + candidate.internal_force.Size() != full_count || + candidate.residual.Size() != full_count || + candidate.reaction.Size() != full_count) { + return RecoveryFailure( + "invalid-recovery-dimensions", {}, state.Identity().step_name, + "Recovery candidate vectors must match the AnalysisState dimension."); + } + if (!IsFiniteVector(candidate.displacement) || + !IsFiniteVector(candidate.external_force) || + !IsFiniteVector(candidate.internal_force) || + !IsFiniteVector(candidate.residual) || + !IsFiniteVector(candidate.reaction)) { + return RecoveryFailure( + "nonfinite-recovery-value", {}, state.Identity().step_name, + "Recovery candidate vectors must contain only finite values."); + } + + const bool has_any_beam_rows = !candidate.endpoint_rows.empty() || + !candidate.gauss_rows.empty() || + !candidate.stress_rows.empty(); + const bool has_complete_beam_rows = !candidate.endpoint_rows.empty() && + !candidate.gauss_rows.empty() && + !candidate.stress_rows.empty(); + if (candidate.has_beam_results != has_complete_beam_rows || + (has_any_beam_rows && !has_complete_beam_rows)) { + return RecoveryFailure( + "invalid-recovery-result-family", {}, state.Identity().step_name, + "Beam recovery flags require endpoint, Gauss, and stress rows as one " + "complete result family."); + } + + const std::array zero_equilibrium{}; + const std::array zero_metrics{}; + const bool has_any_shell_results = + !candidate.expected_shell_elements.empty() || + !candidate.shell.rows.empty() || + candidate.shell.physical_strain_energy != 0.0 || + candidate.shell.equilibrium != zero_equilibrium || + candidate.shell.verification_metrics != zero_metrics; + const bool has_complete_shell_inventory = + !candidate.expected_shell_elements.empty() && + !candidate.shell.rows.empty(); + if (candidate.has_shell_results != has_complete_shell_inventory || + (has_any_shell_results && !has_complete_shell_inventory)) { + return RecoveryFailure( + "invalid-recovery-result-family", {}, state.Identity().step_name, + "Shell recovery flags require expected elements and physical rows as " + "one complete result family."); + } + if (candidate.has_beam_results && candidate.has_shell_results) { + return RecoveryFailure( + "unsupported-mixed-element-model", {}, "mixed-runtime-results", + "Result recovery does not support mixed beam and shell models."); + } + + AnalysisState candidate_state = state; + candidate_state.Displacement() = std::move(candidate.displacement); + candidate_state.ExternalForce() = std::move(candidate.external_force); + candidate_state.InternalForce() = std::move(candidate.internal_force); + candidate_state.Residual() = std::move(candidate.residual); + candidate_state.Reaction() = std::move(candidate.reaction); + candidate_state.EndpointResults() = std::move(candidate.endpoint_rows); + candidate_state.GaussResults() = std::move(candidate.gauss_rows); + candidate_state.StressResults() = std::move(candidate.stress_rows); + const Status shell_commit_status = candidate_state.CommitShellResults( + candidate.expected_shell_elements, std::move(candidate.shell)); + if (!shell_commit_status.IsOk()) { + return shell_commit_status; + } + state = std::move(candidate_state); + return Status::Ok(); +} + +} // namespace fesa::results_internal diff --git a/src/fesa/results/analysis_state_commit.h b/src/fesa/results/analysis_state_commit.h new file mode 100644 index 0000000..bb32464 --- /dev/null +++ b/src/fesa/results/analysis_state_commit.h @@ -0,0 +1,16 @@ +#ifndef FESA_RESULTS_ANALYSIS_STATE_COMMIT_H_ +#define FESA_RESULTS_ANALYSIS_STATE_COMMIT_H_ + +#include "fesa/analysis/analysis_state.h" +#include "fesa/core/status.h" +#include "results/recovery_candidate.h" + +namespace fesa::results_internal { + +/// @brief Validates a complete recovery candidate and atomically updates state. +Status CommitRecoveryCandidate(RecoveryCandidate candidate, + AnalysisState& state); + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_ANALYSIS_STATE_COMMIT_H_ diff --git a/src/fesa/results/beam_result_recovery.cpp b/src/fesa/results/beam_result_recovery.cpp new file mode 100644 index 0000000..8db716f --- /dev/null +++ b/src/fesa/results/beam_result_recovery.cpp @@ -0,0 +1,106 @@ +#include "results/beam_result_recovery.h" + +#include +#include +#include +#include +#include + +namespace fesa::results_internal { + +Status BeamResultRecovery(const SourceLocation& location, + const SourceEntityId& source_id, + const EntityIndex expected_element, + const std::array& expected_nodes, + const BeamElementResultRows& rows, + RecoveryCandidate& candidate) { + if (rows.endpoint_rows.size() != expected_nodes.size()) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Beam endpoint rows must preserve endpoint 0 and 1 source order."); + } + for (std::size_t endpoint = 0U; endpoint < expected_nodes.size(); + ++endpoint) { + const auto& row = rows.endpoint_rows[endpoint]; + if (row.element != expected_element || + row.endpoint != static_cast(endpoint) || + !SameSourceIdentity(row.node, expected_nodes[endpoint])) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Beam endpoint rows must preserve active element, endpoint, and " + "source-node order."); + } + if (!IsFiniteArray(row.end_action) || + !IsFiniteArray(row.section_resultant)) { + return RecoveryFailure( + "nonfinite-recovery-value", location, source_id.source_label_text, + "Beam endpoint actions and section resultants must remain finite."); + } + } + constexpr std::array expected_gauss_points{1, 2}; + if (rows.gauss_rows.size() != expected_gauss_points.size()) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Beam Gauss rows must preserve GP1 and GP2 source order."); + } + for (std::size_t point = 0U; point < expected_gauss_points.size(); ++point) { + const auto& row = rows.gauss_rows[point]; + if (row.element != expected_element || + row.gauss_point != expected_gauss_points[point]) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Beam Gauss rows must preserve active element and GP1 through GP2 " + "order."); + } + if (!IsFiniteArray(row.generalized_strain) || + !IsFiniteArray(row.generalized_resultant)) { + return RecoveryFailure( + "nonfinite-recovery-value", location, source_id.source_label_text, + "Beam generalized recovery values must remain finite."); + } + } + if (rows.stress_rows.empty() || rows.stress_rows.front().gauss_point != 1 || + rows.stress_rows.back().gauss_point != 2) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Beam stress rows must preserve GP and section-point source order."); + } + for (std::size_t order = 0U; order < rows.stress_rows.size(); ++order) { + const auto& row = rows.stress_rows[order]; + if (row.element != expected_element || row.gauss_point < 1 || + row.gauss_point > 2 || + (order != 0U && std::tie(rows.stress_rows[order - 1U].gauss_point, + rows.stress_rows[order - 1U].section_point) >= + std::tie(row.gauss_point, row.section_point))) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Beam stress rows must preserve active element, GP, and " + "section-point order."); + } + if (!std::isfinite(row.x1) || !std::isfinite(row.x2) || + !std::isfinite(row.s11)) { + return RecoveryFailure( + "nonfinite-recovery-value", location, source_id.source_label_text, + "Beam section-point stress values must remain finite."); + } + } + + candidate.has_beam_results = true; + candidate.endpoint_rows.insert(candidate.endpoint_rows.end(), + rows.endpoint_rows.begin(), + rows.endpoint_rows.end()); + candidate.gauss_rows.insert(candidate.gauss_rows.end(), + rows.gauss_rows.begin(), rows.gauss_rows.end()); + candidate.stress_rows.insert(candidate.stress_rows.end(), + rows.stress_rows.begin(), + rows.stress_rows.end()); + return Status::Ok(); +} + +} // namespace fesa::results_internal diff --git a/src/fesa/results/beam_result_recovery.h b/src/fesa/results/beam_result_recovery.h new file mode 100644 index 0000000..45938ac --- /dev/null +++ b/src/fesa/results/beam_result_recovery.h @@ -0,0 +1,23 @@ +#ifndef FESA_RESULTS_BEAM_RESULT_RECOVERY_H_ +#define FESA_RESULTS_BEAM_RESULT_RECOVERY_H_ + +#include + +#include "fesa/core/source_identity.h" +#include "fesa/core/status.h" +#include "fesa/elements/element.h" +#include "results/recovery_candidate.h" + +namespace fesa::results_internal { + +/// @brief Validates and appends distinct beam recovery row collections. +Status BeamResultRecovery(const SourceLocation& location, + const SourceEntityId& source_id, + EntityIndex expected_element, + const std::array& expected_nodes, + const BeamElementResultRows& rows, + RecoveryCandidate& candidate); + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_BEAM_RESULT_RECOVERY_H_ diff --git a/src/fesa/results/global_equilibrium_recovery.cpp b/src/fesa/results/global_equilibrium_recovery.cpp new file mode 100644 index 0000000..1fa096b --- /dev/null +++ b/src/fesa/results/global_equilibrium_recovery.cpp @@ -0,0 +1,269 @@ +#include "results/global_equilibrium_recovery.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/math/vector3.h" + +namespace fesa::results_internal { +namespace { + +constexpr std::size_t kDofsPerNode = 6U; +constexpr double kFreeResidualTolerance = 1.0e-10; +constexpr double kGlobalEquilibriumTolerance = 1.0e-10; + +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; +} + +/// @brief Computes separate Kff*df and Kfc*dc norms in stable CSR order. +std::array FreeEquationInternalTermNorms( + const SparseMatrix& stiffness, const Vector& displacement, + const DofManager& dofs) { + std::vector free_column(dofs.FullDofCount(), 0U); + for (const std::size_t full_dof : dofs.FreeDofs()) { + free_column[full_dof] = 1U; + } + + double free_term_norm = 0.0; + double constrained_term_norm = 0.0; + for (const std::size_t row : dofs.FreeDofs()) { + double free_term = 0.0; + double constrained_term = 0.0; + for (std::size_t position = stiffness.RowOffsets()[row]; + position < stiffness.RowOffsets()[row + 1U]; ++position) { + const std::size_t column = stiffness.ColumnIndices()[position]; + const double contribution = + stiffness.Values()[position] * displacement[column]; + if (free_column[column] != 0U) { + free_term += contribution; + } else { + constrained_term += contribution; + } + } + free_term_norm = std::hypot(free_term_norm, free_term); + constrained_term_norm = std::hypot(constrained_term_norm, constrained_term); + } + return {free_term_norm, constrained_term_norm}; +} + +bool AccumulateVectorAndScale(std::array& total, double& scale, + const std::array& contribution) { + const Vector3 contribution_vector{contribution}; + const double magnitude = contribution_vector.Norm(); + const double accumulated_scale = scale + magnitude; + if (!contribution_vector.IsFinite() || !std::isfinite(magnitude) || + !std::isfinite(accumulated_scale)) { + return false; + } + for (std::size_t component = 0U; component < contribution.size(); + ++component) { + const double accumulated = total[component] + contribution[component]; + if (!std::isfinite(accumulated)) { + return false; + } + total[component] = accumulated; + } + scale = accumulated_scale; + return true; +} + +double NormalizedBalance(const std::array& balance, + const double scale) { + const double balance_norm = Vector3{balance}.Norm(); + if (!std::isfinite(balance_norm) || !std::isfinite(scale)) { + return (std::numeric_limits::infinity)(); + } + if (scale == 0.0) { + return balance_norm == 0.0 ? 0.0 + : (std::numeric_limits::infinity)(); + } + return balance_norm / scale; +} + +} // namespace + +Result GlobalEquilibriumRecovery( + const Domain& domain, const DofManager& dofs, + const SparseMatrix& full_stiffness, const Vector& full_displacement, + const Vector& full_external_force) { + const std::size_t full_count = dofs.FullDofCount(); + if (full_stiffness.Rows() != full_count || + full_stiffness.Columns() != full_count || + full_displacement.Size() != full_count || + full_external_force.Size() != full_count || + domain.Nodes().size() * kDofsPerNode != full_count) { + return RecoveryResultFailure( + "invalid-recovery-dimensions", {domain.SourcePath(), 0U}, + domain.SourceContentIdentity(), + "Model, DOF, stiffness, and recovery vector dimensions must agree."); + } + const Status matrix_status = full_stiffness.Validate(); + if (!matrix_status.IsOk()) { + return Result::Failure(matrix_status); + } + if (!IsFiniteVector(full_displacement) || + !IsFiniteVector(full_external_force)) { + return RecoveryResultFailure( + "nonfinite-recovery-value", {domain.SourcePath(), 0U}, + domain.SourceContentIdentity(), + "Displacement and external-force inputs must be finite."); + } + + GlobalEquilibriumResult result{full_stiffness.Multiply(full_displacement), + Vector{full_count}, Vector{full_count}, 0.0}; + if (!IsFiniteVector(result.internal_force)) { + return RecoveryResultFailure( + "nonfinite-recovery-value", {domain.SourcePath(), 0U}, + domain.SourceContentIdentity(), + "Full stiffness multiplication must produce finite internal force."); + } + for (std::size_t full_dof = 0U; full_dof < result.residual.Size(); + ++full_dof) { + result.residual[full_dof] = + result.internal_force[full_dof] - full_external_force[full_dof]; + if (!std::isfinite(result.residual[full_dof])) { + return RecoveryResultFailure( + "nonfinite-recovery-value", {domain.SourcePath(), 0U}, + std::to_string(full_dof), + "Internal-minus-external residual must remain finite."); + } + } + + const double residual_norm = IndexedNorm(result.residual, dofs.FreeDofs()); + const auto internal_term_norms = + FreeEquationInternalTermNorms(full_stiffness, full_displacement, dofs); + const double external_norm = + IndexedNorm(full_external_force, dofs.FreeDofs()); + // Normalize against Kff*df, Kfc*dc, and Ff separately; the already-cancelled + // K*d term would hide prescribed-displacement scaling. + const double denominator = + (std::max)({internal_term_norms[0U], internal_term_norms[1U], + external_norm}); + if (!std::isfinite(residual_norm) || !std::isfinite(denominator)) { + return RecoveryResultFailure( + "nonfinite-recovery-value", {domain.SourcePath(), 0U}, "free-residual", + "Free residual and its physical normalization scale must be finite."); + } + result.normalized_free_residual = + denominator == 0.0 + ? (residual_norm == 0.0 ? 0.0 + : (std::numeric_limits::infinity)()) + : residual_norm / denominator; + if (!std::isfinite(result.normalized_free_residual) || + result.normalized_free_residual > kFreeResidualTolerance) { + return RecoveryResultFailure( + "free-residual-tolerance-failure", {domain.SourcePath(), 0U}, + "free-residual", "The normalized free residual exceeds 1e-10."); + } + + // The state reaction vector intentionally preserves free residual evidence; + // shell global equilibrium below separately filters constrained reactions. + result.reaction = result.residual; + return Result::Success(std::move(result)); +} + +Status RecoverShellGlobalEvidence(const Domain& domain, const DofManager& dofs, + const Vector& external_force, + const GlobalEquilibriumResult& global, + RecoveryCandidate& candidate) { + std::vector constrained(dofs.FullDofCount(), 0U); + for (const std::size_t full_dof : dofs.ConstrainedDofs()) { + constrained[full_dof] = 1U; + } + + std::array applied_force{}; + std::array reaction_force{}; + std::array applied_moment{}; + std::array reaction_moment{}; + double applied_force_scale = 0.0; + double reaction_force_scale = 0.0; + double applied_moment_scale = 0.0; + double reaction_moment_scale = 0.0; + for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { + const std::size_t offset = node * kDofsPerNode; + std::array nodal_applied_force{}; + std::array nodal_reaction_force{}; + std::array nodal_applied_moment{}; + std::array nodal_reaction_moment{}; + for (std::size_t component = 0U; component < 3U; ++component) { + nodal_applied_force[component] = external_force[offset + component]; + nodal_applied_moment[component] = external_force[offset + 3U + component]; + if (constrained[offset + component] != 0U) { + nodal_reaction_force[component] = global.residual[offset + component]; + } + if (constrained[offset + 3U + component] != 0U) { + nodal_reaction_moment[component] = + global.residual[offset + 3U + component]; + } + } + + const Vector3 position{domain.Nodes()[node].coordinates}; + const Vector3 applied_force_moment = + position.Cross(Vector3{nodal_applied_force}); + const Vector3 reaction_force_moment = + position.Cross(Vector3{nodal_reaction_force}); + for (std::size_t component = 0U; component < 3U; ++component) { + nodal_applied_moment[component] += applied_force_moment[component]; + nodal_reaction_moment[component] += reaction_force_moment[component]; + } + if (!AccumulateVectorAndScale(applied_force, applied_force_scale, + nodal_applied_force) || + !AccumulateVectorAndScale(reaction_force, reaction_force_scale, + nodal_reaction_force) || + !AccumulateVectorAndScale(applied_moment, applied_moment_scale, + nodal_applied_moment) || + !AccumulateVectorAndScale(reaction_moment, reaction_moment_scale, + nodal_reaction_moment)) { + return RecoveryFailure("nonfinite-recovery-value", + domain.Nodes()[node].location, + domain.Nodes()[node].source_id.source_label_text, + "Global force and moment evidence must remain " + "finite in source-node order."); + } + } + + const Vector3 force_balance = + Vector3{applied_force} + Vector3{reaction_force}; + const Vector3 moment_balance = + Vector3{applied_moment} + Vector3{reaction_moment}; + for (std::size_t component = 0U; component < 3U; ++component) { + candidate.shell.equilibrium[component] = force_balance[component]; + candidate.shell.equilibrium[3U + component] = moment_balance[component]; + } + const double force_metric = + NormalizedBalance(force_balance.Components(), + (std::max)(applied_force_scale, reaction_force_scale)); + const double moment_metric = NormalizedBalance( + moment_balance.Components(), + (std::max)(applied_moment_scale, reaction_moment_scale)); + candidate.shell.verification_metrics = {global.normalized_free_residual, + force_metric, moment_metric}; + if (!IsFiniteArray(candidate.shell.equilibrium) || + !Vector3{candidate.shell.verification_metrics}.IsFinite()) { + return RecoveryFailure("nonfinite-recovery-value", + {domain.SourcePath(), 0U}, "global-equilibrium", + "Global equilibrium values and their physical " + "normalization scales must be finite."); + } + if (force_metric > kGlobalEquilibriumTolerance || + moment_metric > kGlobalEquilibriumTolerance) { + return RecoveryFailure( + "global-equilibrium-tolerance-failure", {domain.SourcePath(), 0U}, + "global-equilibrium", + "Normalized global force or moment balance exceeds 1e-10."); + } + return Status::Ok(); +} + +} // namespace fesa::results_internal diff --git a/src/fesa/results/global_equilibrium_recovery.h b/src/fesa/results/global_equilibrium_recovery.h new file mode 100644 index 0000000..4119c39 --- /dev/null +++ b/src/fesa/results/global_equilibrium_recovery.h @@ -0,0 +1,37 @@ +#ifndef FESA_RESULTS_GLOBAL_EQUILIBRIUM_RECOVERY_H_ +#define FESA_RESULTS_GLOBAL_EQUILIBRIUM_RECOVERY_H_ + +#include "fesa/core/status.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/math/vector.h" +#include "fesa/model/domain.h" +#include "results/recovery_candidate.h" + +namespace fesa::results_internal { + +/// @brief Owns full residual, reaction evidence, and free equilibrium scale. +struct GlobalEquilibriumResult { + Vector internal_force; + Vector residual; + Vector reaction; + double normalized_free_residual{0.0}; +}; + +/// @brief Computes K*d-F, full reaction evidence, and free residual metric. +/// @return The complete global residual candidate or a structured failure. +Result GlobalEquilibriumRecovery( + const Domain& domain, const DofManager& dofs, + const SparseMatrix& full_stiffness, const Vector& full_displacement, + const Vector& full_external_force); + +/// @brief Adds shell global force and moment equilibrium evidence. +/// @note Moments are taken about the fixed global origin. +Status RecoverShellGlobalEvidence(const Domain& domain, const DofManager& dofs, + const Vector& external_force, + const GlobalEquilibriumResult& global, + RecoveryCandidate& candidate); + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_GLOBAL_EQUILIBRIUM_RECOVERY_H_ diff --git a/src/fesa/results/recovery_candidate.cpp b/src/fesa/results/recovery_candidate.cpp new file mode 100644 index 0000000..d16300c --- /dev/null +++ b/src/fesa/results/recovery_candidate.cpp @@ -0,0 +1,39 @@ +#include "results/recovery_candidate.h" + +#include +#include + +namespace fesa::results_internal { + +RecoveryCandidate::RecoveryCandidate(const std::size_t full_dof_count) + : displacement{full_dof_count}, + external_force{full_dof_count}, + internal_force{full_dof_count}, + residual{full_dof_count}, + reaction{full_dof_count} {} + +Status RecoveryFailure(const std::string& code, const SourceLocation& location, + const std::string& identity, + const std::string& message) { + return Status::Failure(FailureCategory::kModel, + {{Severity::kError, code, location, "RESULT_RECOVERY", + identity, message}}); +} + +bool SameSourceIdentity(const SourceEntityId& left, + const SourceEntityId& right) { + return left.instance_name == right.instance_name && + left.source_label == right.source_label && + left.source_label_text == right.source_label_text; +} + +bool IsFiniteVector(const Vector& values) { + for (std::size_t index = 0U; index < values.Size(); ++index) { + if (!std::isfinite(values[index])) { + return false; + } + } + return true; +} + +} // namespace fesa::results_internal diff --git a/src/fesa/results/recovery_candidate.h b/src/fesa/results/recovery_candidate.h new file mode 100644 index 0000000..ab11fcd --- /dev/null +++ b/src/fesa/results/recovery_candidate.h @@ -0,0 +1,65 @@ +#ifndef FESA_RESULTS_RECOVERY_CANDIDATE_H_ +#define FESA_RESULTS_RECOVERY_CANDIDATE_H_ + +#include +#include +#include +#include +#include +#include + +#include "fesa/core/source_identity.h" +#include "fesa/core/status.h" +#include "fesa/math/vector.h" +#include "fesa/results/result_records.h" + +namespace fesa::results_internal { + +/// @brief Owns every mutable recovery value before final state commit. +struct RecoveryCandidate { + explicit RecoveryCandidate(std::size_t full_dof_count); + + Vector displacement; + Vector external_force; + Vector internal_force; + Vector residual; + Vector reaction; + std::vector endpoint_rows; + std::vector gauss_rows; + std::vector stress_rows; + ShellStateCandidate shell; + std::vector expected_shell_elements; + bool has_beam_results{false}; + bool has_shell_results{false}; +}; + +/// @brief Creates a structured result-recovery model failure. +Status RecoveryFailure(const std::string& code, const SourceLocation& location, + const std::string& identity, const std::string& message); + +/// @brief Creates a failed result-recovery Result with one diagnostic. +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)); +} + +/// @brief Compares preserved source identity fields exactly. +bool SameSourceIdentity(const SourceEntityId& left, + const SourceEntityId& right); + +/// @brief Tests a fixed-size recovery component array for finite values. +template +bool IsFiniteArray(const std::array& values) { + return std::all_of(values.begin(), values.end(), + [](const double value) { return std::isfinite(value); }); +} + +/// @brief Tests an owned Vector for finite values. +bool IsFiniteVector(const Vector& values); + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_RECOVERY_CANDIDATE_H_ diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp index 4000db6..854c79d 100644 --- a/src/fesa/results/result_recovery.cpp +++ b/src/fesa/results/result_recovery.cpp @@ -18,13 +18,16 @@ #include "fesa/loads/concentrated_nodal_load.h" #include "fesa/math/vector3.h" #include "fesa/model/source_target_resolver.h" +#include "results/analysis_state_commit.h" +#include "results/beam_result_recovery.h" +#include "results/global_equilibrium_recovery.h" +#include "results/recovery_candidate.h" +#include "results/shell_result_recovery.h" namespace fesa { namespace { constexpr std::size_t kDofsPerNode = 6U; -constexpr double kFreeResidualTolerance = 1.0e-10; -constexpr double kGlobalEquilibriumTolerance = 1.0e-10; constexpr double kAxisTolerance = 1.0e-12; using AxisSet = std::array, 3>; @@ -45,66 +48,9 @@ Result RecoveryResultFailure(const std::string& code, return Result::Failure(RecoveryFailure(code, location, identity, message)); } -bool SameSourceIdentity(const SourceEntityId& left, - const SourceEntityId& right) { - return left.instance_name == right.instance_name && - left.source_label == right.source_label && - left.source_label_text == right.source_label_text; -} - template bool IsFinite(const std::array& values) { - return std::all_of(values.begin(), values.end(), - [](const double value) { return std::isfinite(value); }); -} - -bool IsFinite(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; -} - -/// @brief Computes separate Kff*df and Kfc*dc norms in stable CSR order. -std::array FreeEquationInternalTermNorms( - const SparseMatrix& stiffness, const Vector& displacement, - const DofManager& dofs) { - std::vector free_column(dofs.FullDofCount(), 0U); - for (const std::size_t full_dof : dofs.FreeDofs()) { - free_column[full_dof] = 1U; - } - - double free_term_norm = 0.0; - double constrained_term_norm = 0.0; - for (const std::size_t row : dofs.FreeDofs()) { - double free_term = 0.0; - double constrained_term = 0.0; - for (std::size_t position = stiffness.RowOffsets()[row]; - position < stiffness.RowOffsets()[row + 1U]; ++position) { - const std::size_t column = stiffness.ColumnIndices()[position]; - const double contribution = - stiffness.Values()[position] * displacement[column]; - if (free_column[column] != 0U) { - free_term += contribution; - } else { - constrained_term += contribution; - } - } - free_term_norm = std::hypot(free_term_norm, free_term); - constrained_term_norm = std::hypot(constrained_term_norm, constrained_term); - } - return {free_term_norm, constrained_term_norm}; + return results_internal::IsFiniteArray(values); } /// @brief Checks that equation-space indices preserve stable full-DOF order. @@ -157,7 +103,8 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, if (!dof_status.IsOk()) { return dof_status; } - if (!IsFinite(full_displacement) || !IsFinite(full_external_force)) { + if (!results_internal::IsFiniteVector(full_displacement) || + !results_internal::IsFiniteVector(full_external_force)) { return RecoveryFailure( "nonfinite-recovery-value", {domain.SourcePath(), 0U}, domain.SourceContentIdentity(), @@ -244,7 +191,8 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, previous_definition = definition_index; const auto& definition = domain.Elements()[definition_index]; const auto& layout = elements[source_order].get().DofLayout(); - if (!SameSourceIdentity(layout.source_id, definition.SourceId()) || + if (!results_internal::SameSourceIdentity(layout.source_id, + definition.SourceId()) || layout.node_indices != definition.NodeIndices()) { return RecoveryFailure( "invalid-recovery-entity", {domain.SourcePath(), 0U}, @@ -260,90 +208,6 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, return Status::Ok(); } -/// @brief Appends one typed beam bundle without changing row identity or sign. -Status AppendBeamRows(const SourceLocation& location, - const SourceEntityId& source_id, - const EntityIndex expected_element, - const BeamElementResultRows& rows, - std::vector& endpoint_rows, - std::vector& gauss_rows, - std::vector& stress_rows) { - for (const auto& row : rows.endpoint_rows) { - if (row.element != expected_element) { - return RecoveryFailure( - "invalid-element-result-identity", location, - source_id.source_label_text, - "Beam endpoint rows must retain their active element identity."); - } - if (!IsFinite(row.end_action) || !IsFinite(row.section_resultant)) { - return RecoveryFailure( - "nonfinite-recovery-value", location, source_id.source_label_text, - "Beam endpoint actions and section resultants must remain finite."); - } - endpoint_rows.push_back(row); - } - for (const auto& row : rows.gauss_rows) { - if (row.element != expected_element) { - return RecoveryFailure( - "invalid-element-result-identity", location, - source_id.source_label_text, - "Beam Gauss rows must retain their active element identity."); - } - if (!IsFinite(row.generalized_strain) || - !IsFinite(row.generalized_resultant)) { - return RecoveryFailure( - "nonfinite-recovery-value", location, source_id.source_label_text, - "Beam generalized recovery values must remain finite."); - } - gauss_rows.push_back(row); - } - for (const auto& row : rows.stress_rows) { - if (row.element != expected_element) { - return RecoveryFailure( - "invalid-element-result-identity", location, - source_id.source_label_text, - "Beam stress rows must retain their active element identity."); - } - if (!std::isfinite(row.x1) || !std::isfinite(row.x2) || - !std::isfinite(row.s11)) { - return RecoveryFailure( - "nonfinite-recovery-value", location, source_id.source_label_text, - "Beam section-point stress values must remain finite."); - } - stress_rows.push_back(row); - } - return Status::Ok(); -} - -/// @brief Appends one physical shell bundle in stable runtime order. -Status AppendShellRows(const SourceLocation& location, - const SourceEntityId& source_id, - const EntityIndex expected_element, - const ShellElementResultRows& rows, - std::vector& expected_elements, - ShellStateCandidate& candidate) { - const double accumulated_energy = - candidate.physical_strain_energy + rows.physical_strain_energy; - if (!std::isfinite(rows.physical_strain_energy) || - !std::isfinite(accumulated_energy)) { - return RecoveryFailure( - "nonfinite-recovery-value", location, source_id.source_label_text, - "Source-order physical shell energy reduction must remain finite."); - } - for (const auto& row : rows.rows) { - if (row.element != expected_element) { - return RecoveryFailure( - "invalid-element-result-identity", location, - source_id.source_label_text, - "Shell rows must retain their active element identity."); - } - candidate.rows.push_back(row); - } - expected_elements.push_back(expected_element); - candidate.physical_strain_energy = accumulated_energy; - return Status::Ok(); -} - Result> ResolveLoadTarget( const SourceTargetResolver& resolver, const Domain& domain, const ConcentratedNodalLoad& load) { @@ -371,134 +235,6 @@ Result> ResolveLoadTarget( return Result>::Success(std::move(nodes)); } -bool AccumulateVectorAndScale(std::array& total, double& scale, - const std::array& contribution) { - const Vector3 contribution_vector{contribution}; - const double magnitude = contribution_vector.Norm(); - const double accumulated_scale = scale + magnitude; - if (!contribution_vector.IsFinite() || !std::isfinite(magnitude) || - !std::isfinite(accumulated_scale)) { - return false; - } - for (std::size_t component = 0U; component < contribution.size(); - ++component) { - const double accumulated = total[component] + contribution[component]; - if (!std::isfinite(accumulated)) { - return false; - } - total[component] = accumulated; - } - scale = accumulated_scale; - return true; -} - -double NormalizedBalance(const std::array& balance, - const double scale) { - const double balance_norm = Vector3{balance}.Norm(); - if (!std::isfinite(balance_norm) || !std::isfinite(scale)) { - return (std::numeric_limits::infinity)(); - } - if (scale == 0.0) { - return balance_norm == 0.0 ? 0.0 - : (std::numeric_limits::infinity)(); - } - return balance_norm / scale; -} - -/// @brief Computes global shell evidence in stable source order. -Status PopulateShellGlobalEvidence(const Domain& domain, const DofManager& dofs, - const Vector& external_force, - const Vector& residual, - const double normalized_residual, - ShellStateCandidate& candidate) { - std::vector constrained(dofs.FullDofCount(), 0U); - for (const std::size_t full_dof : dofs.ConstrainedDofs()) { - constrained[full_dof] = 1U; - } - - std::array applied_force{}; - std::array reaction_force{}; - std::array applied_moment{}; - std::array reaction_moment{}; - double applied_force_scale = 0.0; - double reaction_force_scale = 0.0; - double applied_moment_scale = 0.0; - double reaction_moment_scale = 0.0; - for (std::size_t node = 0U; node < domain.Nodes().size(); ++node) { - const std::size_t offset = node * kDofsPerNode; - std::array nodal_applied_force{}; - std::array nodal_reaction_force{}; - std::array nodal_applied_moment{}; - std::array nodal_reaction_moment{}; - for (std::size_t component = 0U; component < 3U; ++component) { - nodal_applied_force[component] = external_force[offset + component]; - nodal_applied_moment[component] = external_force[offset + 3U + component]; - if (constrained[offset + component] != 0U) { - nodal_reaction_force[component] = residual[offset + component]; - } - if (constrained[offset + 3U + component] != 0U) { - nodal_reaction_moment[component] = residual[offset + 3U + component]; - } - } - - const Vector3 position{domain.Nodes()[node].coordinates}; - const Vector3 applied_force_moment = - position.Cross(Vector3{nodal_applied_force}); - const Vector3 reaction_force_moment = - position.Cross(Vector3{nodal_reaction_force}); - for (std::size_t component = 0U; component < 3U; ++component) { - nodal_applied_moment[component] += applied_force_moment[component]; - nodal_reaction_moment[component] += reaction_force_moment[component]; - } - if (!AccumulateVectorAndScale(applied_force, applied_force_scale, - nodal_applied_force) || - !AccumulateVectorAndScale(reaction_force, reaction_force_scale, - nodal_reaction_force) || - !AccumulateVectorAndScale(applied_moment, applied_moment_scale, - nodal_applied_moment) || - !AccumulateVectorAndScale(reaction_moment, reaction_moment_scale, - nodal_reaction_moment)) { - return RecoveryFailure("nonfinite-recovery-value", - domain.Nodes()[node].location, - domain.Nodes()[node].source_id.source_label_text, - "Global force and moment evidence must remain " - "finite in source-node order."); - } - } - - const Vector3 force_balance = - Vector3{applied_force} + Vector3{reaction_force}; - const Vector3 moment_balance = - Vector3{applied_moment} + Vector3{reaction_moment}; - for (std::size_t component = 0U; component < 3U; ++component) { - candidate.equilibrium[component] = force_balance[component]; - candidate.equilibrium[3U + component] = moment_balance[component]; - } - const double force_metric = - NormalizedBalance(force_balance.Components(), - (std::max)(applied_force_scale, reaction_force_scale)); - const double moment_metric = NormalizedBalance( - moment_balance.Components(), - (std::max)(applied_moment_scale, reaction_moment_scale)); - candidate.verification_metrics = {normalized_residual, force_metric, - moment_metric}; - if (!IsFinite(candidate.equilibrium) || - !Vector3{candidate.verification_metrics}.IsFinite()) { - return RecoveryFailure("nonfinite-recovery-value", - {domain.SourcePath(), 0U}, "global-equilibrium", - "Global equilibrium values and their physical " - "normalization scales must be finite."); - } - if (force_metric > kGlobalEquilibriumTolerance || - moment_metric > kGlobalEquilibriumTolerance) { - return RecoveryFailure( - "global-equilibrium-tolerance-failure", {domain.SourcePath(), 0U}, - "global-equilibrium", - "Normalized global force or moment balance exceeds 1e-10."); - } - return Status::Ok(); -} - std::optional LocalAxes(const Domain& domain, const EntityIndex first_node, const EntityIndex second_node, @@ -562,67 +298,19 @@ Status ResultRecovery::Recover(const AnalysisModel& model, return input_status; } - Vector internal_force = full_stiffness.Multiply(displacement); - if (!IsFinite(internal_force)) { - return RecoveryFailure( - "nonfinite-recovery-value", {model.GetDomain().SourcePath(), 0U}, - model.GetDomain().SourceContentIdentity(), - "Full stiffness multiplication must produce finite internal force."); - } - Vector residual{dofs.FullDofCount()}; - for (std::size_t full_dof = 0U; full_dof < residual.Size(); ++full_dof) { - residual[full_dof] = internal_force[full_dof] - external_force[full_dof]; - if (!std::isfinite(residual[full_dof])) { - return RecoveryFailure( - "nonfinite-recovery-value", {model.GetDomain().SourcePath(), 0U}, - std::to_string(full_dof), - "Internal-minus-external residual must remain finite."); - } - } - - const double residual_norm = IndexedNorm(residual, dofs.FreeDofs()); - const auto internal_term_norms = - FreeEquationInternalTermNorms(full_stiffness, displacement, dofs); - const double external_norm = IndexedNorm(external_force, dofs.FreeDofs()); - // Normalize against the three terms of - // Kff*df + Kfc*dc - Ff. Using only the already-cancelled K*d term would - // classify prescribed-only equilibrium roundoff as a unit residual. - const double denominator = - (std::max)({internal_term_norms[0U], internal_term_norms[1U], - external_norm}); - if (!std::isfinite(residual_norm) || !std::isfinite(denominator)) { - return RecoveryFailure( - "nonfinite-recovery-value", {model.GetDomain().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 normalized_residual = - denominator == 0.0 - ? (residual_norm == 0.0 ? 0.0 - : (std::numeric_limits::infinity)()) - : residual_norm / denominator; - if (!std::isfinite(normalized_residual) || - normalized_residual > kFreeResidualTolerance) { - return RecoveryFailure( - "free-residual-tolerance-failure", {model.GetDomain().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 endpoint_rows; - std::vector gauss_rows; - std::vector stress_rows; const Domain& domain = model.GetDomain(); - ShellStateCandidate shell_candidate{}; - std::vector expected_shell_elements; - bool has_beam_results = false; - bool has_shell_results = false; + auto global = results_internal::GlobalEquilibriumRecovery( + domain, dofs, full_stiffness, displacement, external_force); + if (!global.HasValue()) { + return global.GetStatus(); + } + + results_internal::RecoveryCandidate candidate{dofs.FullDofCount()}; + candidate.displacement = displacement; + candidate.external_force = external_force; + candidate.internal_force = global.Value().internal_force; + candidate.residual = global.Value().residual; + candidate.reaction = global.Value().reaction; for (std::size_t source_order = 0U; source_order < elements.size(); ++source_order) { const EntityIndex element_index = model.ActiveElements()[source_order]; @@ -641,8 +329,8 @@ Status ResultRecovery::Recover(const AnalysisModel& model, if (!recovered.HasValue()) { return recovered.GetStatus(); } - if (!SameSourceIdentity(recovered.Value().source_id, - element.DofLayout().source_id)) { + if (!results_internal::SameSourceIdentity(recovered.Value().source_id, + element.DofLayout().source_id)) { return RecoveryFailure( "invalid-element-result-identity", {domain.SourcePath(), 0U}, element.DofLayout().source_id.source_label_text, @@ -653,16 +341,25 @@ Status ResultRecovery::Recover(const AnalysisModel& model, [&](const auto& payload) -> Status { using Payload = std::decay_t; if constexpr (std::is_same_v) { - has_beam_results = true; - return AppendBeamRows( + const auto& node_indices = element.DofLayout().node_indices; + if (node_indices.size() != 2U || + node_indices[0U] >= domain.Nodes().size() || + node_indices[1U] >= domain.Nodes().size()) { + return results_internal::RecoveryFailure( + "invalid-element-result-identity", {domain.SourcePath(), 0U}, + recovered.Value().source_id.source_label_text, + "Beam recovery requires exactly two valid source nodes."); + } + const std::array expected_nodes{ + domain.Nodes()[node_indices[0U]].source_id, + domain.Nodes()[node_indices[1U]].source_id}; + return results_internal::BeamResultRecovery( {domain.SourcePath(), 0U}, recovered.Value().source_id, - element_index, payload, endpoint_rows, gauss_rows, stress_rows); + element_index, expected_nodes, payload, candidate); } else { - has_shell_results = true; - return AppendShellRows({domain.SourcePath(), 0U}, - recovered.Value().source_id, element_index, - payload, expected_shell_elements, - shell_candidate); + return results_internal::ShellResultRecovery( + {domain.SourcePath(), 0U}, recovered.Value().source_id, + element_index, payload, candidate); } }, recovered.Value().payload); @@ -670,40 +367,21 @@ Status ResultRecovery::Recover(const AnalysisModel& model, return aggregation_status; } } - if (has_beam_results && has_shell_results) { + if (candidate.has_beam_results && candidate.has_shell_results) { return RecoveryFailure( "unsupported-mixed-element-model", {domain.SourcePath(), 0U}, "mixed-runtime-results", "Result recovery does not support mixed beam and shell models."); } - if (has_shell_results) { - const Status evidence_status = - PopulateShellGlobalEvidence(domain, dofs, external_force, residual, - normalized_residual, shell_candidate); + if (candidate.has_shell_results) { + const Status evidence_status = results_internal::RecoverShellGlobalEvidence( + domain, dofs, external_force, global.Value(), candidate); if (!evidence_status.IsOk()) { return evidence_status; } } - // Build and validate a complete candidate state first. This preserves the - // prior full residual, beam rows, and shell rows if any later shell or - // candidate-inventory validation fails. - AnalysisState candidate_state = state; - candidate_state.Displacement() = displacement; - candidate_state.ExternalForce() = external_force; - candidate_state.InternalForce() = std::move(internal_force); - candidate_state.Residual() = std::move(residual); - candidate_state.Reaction() = std::move(reaction); - candidate_state.EndpointResults() = std::move(endpoint_rows); - candidate_state.GaussResults() = std::move(gauss_rows); - candidate_state.StressResults() = std::move(stress_rows); - const Status shell_commit_status = candidate_state.CommitShellResults( - expected_shell_elements, std::move(shell_candidate)); - if (!shell_commit_status.IsOk()) { - return shell_commit_status; - } - state = std::move(candidate_state); - return Status::Ok(); + return results_internal::CommitRecoveryCandidate(std::move(candidate), state); } Status ResultRecovery::Recover(const AnalysisModel& model, @@ -774,7 +452,8 @@ ResultRecovery::NormalizeSectionResultantsToNodeStations( const EntityIndex node_index = definition.node_indices[endpoint]; if (node_index >= domain.Nodes().size() || row.element != element_index || row.endpoint != static_cast(endpoint) || - !SameSourceIdentity(row.node, domain.Nodes()[node_index].source_id)) { + !results_internal::SameSourceIdentity( + row.node, domain.Nodes()[node_index].source_id)) { return RecoveryResultFailure>( "invalid-node-station-entity", definition.location, definition.source_id.source_label_text, diff --git a/src/fesa/results/shell_result_recovery.cpp b/src/fesa/results/shell_result_recovery.cpp new file mode 100644 index 0000000..bd47684 --- /dev/null +++ b/src/fesa/results/shell_result_recovery.cpp @@ -0,0 +1,65 @@ +#include "results/shell_result_recovery.h" + +#include +#include +#include + +namespace fesa::results_internal { +namespace { + +bool IsFiniteShellRow(const ShellResultRow& row) { + if (!IsFiniteArray(row.natural_coordinates) || + !IsFiniteArray(row.generalized_strain) || + !IsFiniteArray(row.section_resultant)) { + return false; + } + for (const auto& axis : row.local_frame) { + if (!IsFiniteArray(axis)) { + return false; + } + } + return std::all_of(row.stress.begin(), row.stress.end(), + [](const ShellSectionStressRow& stress) { + return std::isfinite(stress.zeta) && + IsFiniteArray(stress.components); + }); +} + +} // namespace + +Status ShellResultRecovery(const SourceLocation& location, + const SourceEntityId& source_id, + const EntityIndex expected_element, + const ShellElementResultRows& rows, + RecoveryCandidate& candidate) { + const double accumulated_energy = + candidate.shell.physical_strain_energy + rows.physical_strain_energy; + if (!std::isfinite(rows.physical_strain_energy) || + !std::isfinite(accumulated_energy)) { + return RecoveryFailure( + "nonfinite-recovery-value", location, source_id.source_label_text, + "Source-order physical shell energy reduction must remain finite."); + } + for (const auto& row : rows.rows) { + if (row.element != expected_element) { + return RecoveryFailure( + "invalid-element-result-identity", location, + source_id.source_label_text, + "Shell rows must retain their active element identity."); + } + if (!IsFiniteShellRow(row)) { + return RecoveryFailure( + "nonfinite-recovery-value", location, source_id.source_label_text, + "Shell physical recovery rows must remain finite."); + } + } + + candidate.has_shell_results = true; + candidate.shell.rows.insert(candidate.shell.rows.end(), rows.rows.begin(), + rows.rows.end()); + candidate.expected_shell_elements.push_back(expected_element); + candidate.shell.physical_strain_energy = accumulated_energy; + return Status::Ok(); +} + +} // namespace fesa::results_internal diff --git a/src/fesa/results/shell_result_recovery.h b/src/fesa/results/shell_result_recovery.h new file mode 100644 index 0000000..183eba4 --- /dev/null +++ b/src/fesa/results/shell_result_recovery.h @@ -0,0 +1,20 @@ +#ifndef FESA_RESULTS_SHELL_RESULT_RECOVERY_H_ +#define FESA_RESULTS_SHELL_RESULT_RECOVERY_H_ + +#include "fesa/core/source_identity.h" +#include "fesa/core/status.h" +#include "fesa/elements/element.h" +#include "results/recovery_candidate.h" + +namespace fesa::results_internal { + +/// @brief Validates and appends physical shell rows and energy evidence. +Status ShellResultRecovery(const SourceLocation& location, + const SourceEntityId& source_id, + EntityIndex expected_element, + const ShellElementResultRows& rows, + RecoveryCandidate& candidate); + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_SHELL_RESULT_RECOVERY_H_ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3ce9a0e..824d5fe 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -37,6 +37,7 @@ add_executable( unit/materials/material_test.cpp unit/properties/element_property_test.cpp unit/results/result_records_test.cpp + unit/results/result_recovery_components_test.cpp unit/results/result_recovery_test.cpp unit/results/results_writer_test.cpp unit/solvers/linear/linear_solver_test.cpp diff --git a/tests/unit/results/result_recovery_components_test.cpp b/tests/unit/results/result_recovery_components_test.cpp new file mode 100644 index 0000000..91bcc79 --- /dev/null +++ b/tests/unit/results/result_recovery_components_test.cpp @@ -0,0 +1,593 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/analysis/analysis_state.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/assembly/sparse_assembler.h" +#include "fesa/elements/element.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/vector3.h" +#include "fesa/model/domain.h" +#include "results/analysis_state_commit.h" +#include "results/beam_result_recovery.h" +#include "results/global_equilibrium_recovery.h" +#include "results/recovery_candidate.h" +#include "results/shell_result_recovery.h" + +namespace { + +struct ShellFixture { + std::unique_ptr domain; + std::unique_ptr model; + std::unique_ptr dofs; + std::unique_ptr stiffness; +}; + +fesa::SourceEntityId MakeSourceId(const std::string& text) { + return {"Part-1", 10, text}; +} + +fesa::BeamElementResultRows MakeBeamRows(const fesa::EntityIndex element, + const double value) { + fesa::BeamElementResultRows rows{}; + rows.endpoint_rows = {{element, + 0, + MakeSourceId("1"), + {-value, 0.0, 0.0, 0.0, 0.0, 0.0}, + {value, 2.0 * value, 3.0 * value, 4.0 * value}}, + {element, + 1, + MakeSourceId("2"), + {value, 0.0, 0.0, 0.0, 0.0, 0.0}, + {value, 2.0 * value, 3.0 * value, 4.0 * value}}}; + rows.gauss_rows = {{element, + 1, + {value, value + 1.0, value + 2.0, value + 3.0}, + {value, 2.0 * value, 3.0 * value, 4.0 * value}}, + {element, + 2, + {value + 4.0, value + 5.0, value + 6.0, value + 7.0}, + {value + 8.0, value + 9.0, value + 10.0, value + 11.0}}}; + rows.stress_rows = {{element, 1, 0U, 0.0, 0.0, 5.0 * value, "fake"}, + {element, 2, 1U, -0.25, 0.5, 6.0 * value, "input"}}; + return rows; +} + +fesa::ShellElementResultRows MakeShellRows(const fesa::EntityIndex element, + const double physical_energy) { + const double gauss = 1.0 / std::sqrt(3.0); + const std::array locations{ + fesa::ShellMidsurfaceLocation::kGp1, fesa::ShellMidsurfaceLocation::kGp2, + fesa::ShellMidsurfaceLocation::kGp3, fesa::ShellMidsurfaceLocation::kGp4}; + const std::array, 4> coordinates{ + std::array{-gauss, -gauss}, + std::array{gauss, -gauss}, std::array{gauss, gauss}, + std::array{-gauss, gauss}}; + const std::array positions{ + fesa::ShellSectionPosition::kBottom, fesa::ShellSectionPosition::kMiddle, + fesa::ShellSectionPosition::kTop}; + constexpr std::array zeta{-1.0, 0.0, 1.0}; + + fesa::ShellElementResultRows rows{}; + rows.physical_strain_energy = physical_energy; + for (std::size_t point = 0U; point < locations.size(); ++point) { + fesa::ShellResultRow row{}; + row.element = element; + row.location = locations[point]; + row.natural_coordinates = coordinates[point]; + row.local_frame = {{{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + for (std::size_t component = 0U; component < row.generalized_strain.size(); + ++component) { + row.generalized_strain[component] = + static_cast(point + component + 1U); + row.section_resultant[component] = + 10.0 * static_cast(point + component + 1U); + } + for (std::size_t position = 0U; position < positions.size(); ++position) { + row.stress[position] = {positions[position], + zeta[position], + {static_cast(point + position + 1U), + static_cast(point + position + 2U), + static_cast(point + position + 3U)}}; + } + rows.rows.push_back(std::move(row)); + } + return rows; +} + +fesa::ModelDefinition MakeShellDefinition() { + const std::filesystem::path source{"models/recovery-components-shell.inp"}; + fesa::ModelDefinition definition{}; + definition.source_path = source; + definition.source_content_identity = "fnv1a64:1234567890abcdef"; + definition.nodes = {{{"Shell-1", 1, "1"}, {-1.0, -1.0, 0.0}, {source, 10U}}, + {{"Shell-1", 2, "2"}, {1.0, -1.0, 0.0}, {source, 11U}}, + {{"Shell-1", 3, "3"}, {1.0, 1.0, 0.0}, {source, 12U}}, + {{"Shell-1", 4, "4"}, {-1.0, 1.0, 0.0}, {source, 13U}}}; + definition.materials = {{"Material", 120.0, 0.25, {source, 20U}}}; + definition.shell_sections = {{"ShellSection", 2.0, 0U, {source, 30U}}}; + definition.shell_elements = {{{"Shell-1", 10, "10"}, + fesa::ShellSourceElementType::kS4, + {0U, 1U, 2U, 3U}, + 0U, + 0U, + {source, 40U}}}; + for (std::size_t node = 0U; node < definition.nodes.size(); ++node) { + definition.shell_node_initial_frames.push_back( + {static_cast(node), + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}}); + } + definition.steps = {{"Step-1", {}, {}, 0.1, 1.0, 0.01, 1.0, {source, 50U}}}; + return definition; +} + +ShellFixture MakeShellFixture(fesa::ModelDefinition definition) { + auto domain_result = fesa::Domain::Create(std::move(definition)); + if (!domain_result.HasValue()) { + throw std::runtime_error{"Shell Domain construction failed."}; + } + auto domain = + std::make_unique(std::move(domain_result.Value())); + + auto model_result = fesa::AnalysisModel::Create(*domain); + if (!model_result.HasValue()) { + throw std::runtime_error{"Shell AnalysisModel construction failed."}; + } + auto model = + std::make_unique(std::move(model_result.Value())); + + auto dofs_result = fesa::DofManager::Create(*model); + if (!dofs_result.HasValue()) { + throw std::runtime_error{"Shell DofManager construction failed."}; + } + auto dofs = + std::make_unique(std::move(dofs_result.Value())); + + fesa::SerialParallelFor serial; + auto stiffness_result = + fesa::SparseAssembler::AssembleStiffness(*model, *dofs, serial); + if (!stiffness_result.HasValue()) { + throw std::runtime_error{"Shell stiffness assembly failed."}; + } + auto stiffness = + std::make_unique(std::move(stiffness_result.Value())); + return {std::move(domain), std::move(model), std::move(dofs), + std::move(stiffness)}; +} + +fesa::AnalysisState MakeShellPhysicalState(const ShellFixture& fixture) { + constexpr std::array generalized{0.1, -0.05, 0.2, 0.3, + -0.15, 0.25, 0.4, -0.3}; + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; node < fixture.domain->Nodes().size(); ++node) { + const double x = fixture.domain->Nodes()[node].coordinates[0U]; + const double y = fixture.domain->Nodes()[node].coordinates[1U]; + const std::size_t offset = 6U * node; + state.Displacement()[offset] = + generalized[0U] * x + 0.5 * generalized[2U] * y; + state.Displacement()[offset + 1U] = + generalized[1U] * y + 0.5 * generalized[2U] * x; + state.Displacement()[offset + 2U] = generalized[6U] * x + + generalized[7U] * y - + 0.5 * generalized[5U] * x * y; + state.Displacement()[offset + 3U] = + -generalized[4U] * y - 0.5 * generalized[5U] * x; + state.Displacement()[offset + 4U] = + generalized[3U] * x + 0.5 * generalized[5U] * y; + } + 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.Category(), fesa::FailureCategory::kModel); + ASSERT_EQ(status.Diagnostics().size(), 1U); + EXPECT_EQ(status.Diagnostics()[0U].code, code); +} + +void ExpectVectorEqual(const fesa::Vector& actual, + const fesa::Vector& expected) { + ASSERT_EQ(actual.Size(), expected.Size()); + for (std::size_t index = 0U; index < actual.Size(); ++index) { + EXPECT_DOUBLE_EQ(actual[index], expected[index]); + } +} + +void ExpectAnalysisStateEqual(const fesa::AnalysisState& actual, + const fesa::AnalysisState& expected) { + EXPECT_EQ(actual.Identity().step_name, expected.Identity().step_name); + EXPECT_EQ(actual.Identity().frame_index, expected.Identity().frame_index); + ExpectVectorEqual(actual.Displacement(), expected.Displacement()); + ExpectVectorEqual(actual.ExternalForce(), expected.ExternalForce()); + ExpectVectorEqual(actual.InternalForce(), expected.InternalForce()); + ExpectVectorEqual(actual.Residual(), expected.Residual()); + ExpectVectorEqual(actual.Reaction(), expected.Reaction()); + ASSERT_EQ(actual.EndpointResults().size(), expected.EndpointResults().size()); + for (std::size_t row = 0U; row < actual.EndpointResults().size(); ++row) { + EXPECT_EQ(actual.EndpointResults()[row].element, + expected.EndpointResults()[row].element); + EXPECT_EQ(actual.EndpointResults()[row].endpoint, + expected.EndpointResults()[row].endpoint); + EXPECT_EQ(actual.EndpointResults()[row].node.instance_name, + expected.EndpointResults()[row].node.instance_name); + EXPECT_EQ(actual.EndpointResults()[row].node.source_label, + expected.EndpointResults()[row].node.source_label); + EXPECT_EQ(actual.EndpointResults()[row].node.source_label_text, + expected.EndpointResults()[row].node.source_label_text); + EXPECT_EQ(actual.EndpointResults()[row].end_action, + expected.EndpointResults()[row].end_action); + EXPECT_EQ(actual.EndpointResults()[row].section_resultant, + expected.EndpointResults()[row].section_resultant); + } + ASSERT_EQ(actual.GaussResults().size(), expected.GaussResults().size()); + for (std::size_t row = 0U; row < actual.GaussResults().size(); ++row) { + EXPECT_EQ(actual.GaussResults()[row].element, + expected.GaussResults()[row].element); + EXPECT_EQ(actual.GaussResults()[row].gauss_point, + expected.GaussResults()[row].gauss_point); + EXPECT_EQ(actual.GaussResults()[row].generalized_strain, + expected.GaussResults()[row].generalized_strain); + EXPECT_EQ(actual.GaussResults()[row].generalized_resultant, + expected.GaussResults()[row].generalized_resultant); + } + ASSERT_EQ(actual.StressResults().size(), expected.StressResults().size()); + for (std::size_t row = 0U; row < actual.StressResults().size(); ++row) { + EXPECT_EQ(actual.StressResults()[row].element, + expected.StressResults()[row].element); + EXPECT_EQ(actual.StressResults()[row].gauss_point, + expected.StressResults()[row].gauss_point); + EXPECT_EQ(actual.StressResults()[row].section_point, + expected.StressResults()[row].section_point); + EXPECT_DOUBLE_EQ(actual.StressResults()[row].x1, + expected.StressResults()[row].x1); + EXPECT_DOUBLE_EQ(actual.StressResults()[row].x2, + expected.StressResults()[row].x2); + EXPECT_DOUBLE_EQ(actual.StressResults()[row].s11, + expected.StressResults()[row].s11); + EXPECT_EQ(actual.StressResults()[row].source, + expected.StressResults()[row].source); + } + ASSERT_EQ(actual.ShellResults().size(), expected.ShellResults().size()); + for (std::size_t row = 0U; row < actual.ShellResults().size(); ++row) { + EXPECT_EQ(actual.ShellResults()[row].element, + expected.ShellResults()[row].element); + EXPECT_EQ(actual.ShellResults()[row].location, + expected.ShellResults()[row].location); + EXPECT_EQ(actual.ShellResults()[row].natural_coordinates, + expected.ShellResults()[row].natural_coordinates); + EXPECT_EQ(actual.ShellResults()[row].local_frame, + expected.ShellResults()[row].local_frame); + EXPECT_EQ(actual.ShellResults()[row].generalized_strain, + expected.ShellResults()[row].generalized_strain); + EXPECT_EQ(actual.ShellResults()[row].section_resultant, + expected.ShellResults()[row].section_resultant); + } + EXPECT_DOUBLE_EQ(actual.PhysicalStrainEnergy(), + expected.PhysicalStrainEnergy()); + EXPECT_EQ(actual.Equilibrium(), expected.Equilibrium()); + EXPECT_EQ(actual.VerificationMetrics(), expected.VerificationMetrics()); +} + +} // namespace + +// C-MODULE-002 +TEST(ResultRecoveryComponents, + GlobalEquilibriumUsesFullResidualAndGlobalOriginMoment) { + auto centered_definition = MakeShellDefinition(); + auto translated_definition = centered_definition; + constexpr fesa::Vector3 translation{7.0, 11.0, 0.0}; + for (auto& node : translated_definition.nodes) { + for (std::size_t component = 0U; + component < translation.Components().size(); ++component) { + node.coordinates[component] += translation[component]; + } + } + const auto centered_fixture = + MakeShellFixture(std::move(centered_definition)); + const auto translated_fixture = + MakeShellFixture(std::move(translated_definition)); + auto centered = MakeShellPhysicalState(centered_fixture); + auto translated = MakeShellPhysicalState(translated_fixture); + centered.ExternalForce()[0U] += 1.0e-9; + translated.ExternalForce()[0U] += 1.0e-9; + + auto centered_global = fesa::results_internal::GlobalEquilibriumRecovery( + *centered_fixture.domain, *centered_fixture.dofs, + *centered_fixture.stiffness, centered.Displacement(), + centered.ExternalForce()); + ASSERT_TRUE(centered_global.HasValue()); + auto translated_global = fesa::results_internal::GlobalEquilibriumRecovery( + *translated_fixture.domain, *translated_fixture.dofs, + *translated_fixture.stiffness, translated.Displacement(), + translated.ExternalForce()); + ASSERT_TRUE(translated_global.HasValue()); + EXPECT_NEAR(centered_global.Value().residual[0U], -1.0e-9, 1.0e-14); + EXPECT_DOUBLE_EQ(centered_global.Value().reaction[0U], + centered_global.Value().residual[0U]); + + fesa::results_internal::RecoveryCandidate centered_candidate{ + centered_fixture.dofs->FullDofCount()}; + const auto centered_status = + fesa::results_internal::RecoverShellGlobalEvidence( + *centered_fixture.domain, *centered_fixture.dofs, + centered.ExternalForce(), centered_global.Value(), + centered_candidate); + ASSERT_TRUE(centered_status.IsOk()); + fesa::results_internal::RecoveryCandidate translated_candidate{ + translated_fixture.dofs->FullDofCount()}; + const auto translated_status = + fesa::results_internal::RecoverShellGlobalEvidence( + *translated_fixture.domain, *translated_fixture.dofs, + translated.ExternalForce(), translated_global.Value(), + translated_candidate); + ASSERT_TRUE(translated_status.IsOk()); + + std::array centered_force{}; + for (std::size_t component = 0U; component < centered_force.size(); + ++component) { + centered_force[component] = centered_candidate.shell.equilibrium[component]; + EXPECT_NEAR(translated_candidate.shell.equilibrium[component], + centered_force[component], 1.0e-12); + } + const fesa::Vector3 translated_moment_delta = + translation.Cross(fesa::Vector3{centered_force}); + for (std::size_t component = 0U; component < centered_force.size(); + ++component) { + EXPECT_NEAR(translated_candidate.shell.equilibrium[3U + component] - + centered_candidate.shell.equilibrium[3U + component], + translated_moment_delta[component], 5.0e-12); + } +} + +// C-MODULE-002 +TEST(ResultRecoveryComponents, + BeamResultRecoveryKeepsEndpointGaussAndStressIdentities) { + constexpr fesa::EntityIndex element = 7U; + const std::array expected_nodes{MakeSourceId("1"), + MakeSourceId("2")}; + fesa::results_internal::RecoveryCandidate candidate{0U}; + const auto status = fesa::results_internal::BeamResultRecovery( + {}, MakeSourceId("10"), element, expected_nodes, + MakeBeamRows(element, 2.0), candidate); + + ASSERT_TRUE(status.IsOk()); + ASSERT_EQ(candidate.endpoint_rows.size(), 2U); + ASSERT_EQ(candidate.gauss_rows.size(), 2U); + ASSERT_EQ(candidate.stress_rows.size(), 2U); + EXPECT_TRUE(candidate.has_beam_results); + EXPECT_EQ(candidate.endpoint_rows[0U].endpoint, 0); + EXPECT_EQ(candidate.endpoint_rows[1U].endpoint, 1); + EXPECT_DOUBLE_EQ(candidate.endpoint_rows[0U].end_action[0U], -2.0); + EXPECT_DOUBLE_EQ(candidate.endpoint_rows[0U].section_resultant[0U], 2.0); + EXPECT_DOUBLE_EQ(candidate.endpoint_rows[1U].end_action[0U], 2.0); + EXPECT_DOUBLE_EQ(candidate.endpoint_rows[1U].section_resultant[0U], 2.0); + EXPECT_EQ(candidate.gauss_rows[0U].gauss_point, 1); + EXPECT_EQ(candidate.gauss_rows[1U].gauss_point, 2); + EXPECT_EQ(candidate.stress_rows[0U].section_point, 0U); + EXPECT_EQ(candidate.stress_rows[1U].section_point, 1U); + + auto invalid = MakeBeamRows(element, 3.0); + invalid.stress_rows[1U].s11 = std::numeric_limits::quiet_NaN(); + fesa::results_internal::RecoveryCandidate rejected{0U}; + ExpectStatusCode( + fesa::results_internal::BeamResultRecovery( + {}, MakeSourceId("10"), element, expected_nodes, invalid, rejected), + "nonfinite-recovery-value"); + EXPECT_TRUE(rejected.endpoint_rows.empty()); + EXPECT_TRUE(rejected.gauss_rows.empty()); + EXPECT_TRUE(rejected.stress_rows.empty()); +} + +// C-MODULE-002 +TEST(ResultRecoveryComponents, + BeamResultRecoveryRejectsMalformedB33IdentityBeforeAppending) { + constexpr fesa::EntityIndex element = 7U; + const std::array expected_nodes{MakeSourceId("1"), + MakeSourceId("2")}; + std::vector> invalid; + + auto wrong_endpoint = MakeBeamRows(element, 1.0); + wrong_endpoint.endpoint_rows[1U].endpoint = 2; + invalid.emplace_back("wrong endpoint", std::move(wrong_endpoint)); + auto swapped_endpoints = MakeBeamRows(element, 1.0); + std::swap(swapped_endpoints.endpoint_rows[0U], + swapped_endpoints.endpoint_rows[1U]); + invalid.emplace_back("swapped endpoints", std::move(swapped_endpoints)); + auto missing_endpoint = MakeBeamRows(element, 1.0); + missing_endpoint.endpoint_rows.pop_back(); + invalid.emplace_back("missing endpoint", std::move(missing_endpoint)); + auto duplicate_endpoint = MakeBeamRows(element, 1.0); + duplicate_endpoint.endpoint_rows[1U] = duplicate_endpoint.endpoint_rows[0U]; + invalid.emplace_back("duplicate endpoint", std::move(duplicate_endpoint)); + auto wrong_source_node = MakeBeamRows(element, 1.0); + wrong_source_node.endpoint_rows[0U].node = MakeSourceId("99"); + invalid.emplace_back("wrong source node", std::move(wrong_source_node)); + + auto wrong_gauss = MakeBeamRows(element, 1.0); + wrong_gauss.gauss_rows[1U].gauss_point = 3; + invalid.emplace_back("wrong GP", std::move(wrong_gauss)); + auto swapped_gauss = MakeBeamRows(element, 1.0); + std::swap(swapped_gauss.gauss_rows[0U], swapped_gauss.gauss_rows[1U]); + invalid.emplace_back("swapped GP", std::move(swapped_gauss)); + auto missing_gauss = MakeBeamRows(element, 1.0); + missing_gauss.gauss_rows.pop_back(); + invalid.emplace_back("missing GP", std::move(missing_gauss)); + auto duplicate_gauss = MakeBeamRows(element, 1.0); + duplicate_gauss.gauss_rows[1U] = duplicate_gauss.gauss_rows[0U]; + invalid.emplace_back("duplicate GP", std::move(duplicate_gauss)); + + auto malformed_stress_order = MakeBeamRows(element, 1.0); + std::swap(malformed_stress_order.stress_rows[0U], + malformed_stress_order.stress_rows[1U]); + invalid.emplace_back("malformed stress order", + std::move(malformed_stress_order)); + + for (const auto& [description, rows] : invalid) { + SCOPED_TRACE(description); + fesa::results_internal::RecoveryCandidate candidate{0U}; + ExpectStatusCode( + fesa::results_internal::BeamResultRecovery( + {}, MakeSourceId("10"), element, expected_nodes, rows, candidate), + "invalid-element-result-identity"); + EXPECT_TRUE(candidate.endpoint_rows.empty()); + EXPECT_TRUE(candidate.gauss_rows.empty()); + EXPECT_TRUE(candidate.stress_rows.empty()); + EXPECT_FALSE(candidate.has_beam_results); + } +} + +// C-MODULE-002 +TEST(ResultRecoveryComponents, + ShellResultRecoveryKeepsGpStressOrderAndPhysicalEnergy) { + constexpr fesa::EntityIndex element = 4U; + fesa::results_internal::RecoveryCandidate candidate{0U}; + const auto status = fesa::results_internal::ShellResultRecovery( + {}, MakeSourceId("20"), element, MakeShellRows(element, 12.5), candidate); + + ASSERT_TRUE(status.IsOk()); + EXPECT_TRUE(candidate.has_shell_results); + EXPECT_EQ(candidate.expected_shell_elements, + (std::vector{element})); + ASSERT_EQ(candidate.shell.rows.size(), 4U); + EXPECT_EQ(candidate.shell.rows[0U].location, + fesa::ShellMidsurfaceLocation::kGp1); + EXPECT_EQ(candidate.shell.rows[3U].location, + fesa::ShellMidsurfaceLocation::kGp4); + EXPECT_EQ(candidate.shell.rows[0U].stress[0U].position, + fesa::ShellSectionPosition::kBottom); + EXPECT_EQ(candidate.shell.rows[0U].stress[1U].position, + fesa::ShellSectionPosition::kMiddle); + EXPECT_EQ(candidate.shell.rows[0U].stress[2U].position, + fesa::ShellSectionPosition::kTop); + EXPECT_DOUBLE_EQ(candidate.shell.physical_strain_energy, 12.5); + + auto invalid = MakeShellRows(element, 12.5); + invalid.rows[2U].element = element + 1U; + fesa::results_internal::RecoveryCandidate rejected{0U}; + ExpectStatusCode(fesa::results_internal::ShellResultRecovery( + {}, MakeSourceId("20"), element, invalid, rejected), + "invalid-element-result-identity"); + EXPECT_TRUE(rejected.shell.rows.empty()); + EXPECT_TRUE(rejected.expected_shell_elements.empty()); + EXPECT_DOUBLE_EQ(rejected.shell.physical_strain_energy, 0.0); +} + +// C-MODULE-002 +TEST(ResultRecoveryComponents, + CommitRecoveryCandidateRollsBackLaterInvalidShellBundle) { + const auto fixture = MakeShellFixture(MakeShellDefinition()); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + for (std::size_t full_dof = 0U; full_dof < fixture.dofs->FullDofCount(); + ++full_dof) { + state.Displacement()[full_dof] = static_cast(full_dof) + 0.25; + state.ExternalForce()[full_dof] = -static_cast(full_dof) - 0.5; + state.InternalForce()[full_dof] = 100.0 + static_cast(full_dof); + state.Residual()[full_dof] = 200.0 + static_cast(full_dof); + state.Reaction()[full_dof] = 300.0 + static_cast(full_dof); + } + state.EndpointResults() = MakeBeamRows(0U, 9.0).endpoint_rows; + auto prior_shell_rows = MakeShellRows(0U, 71.0); + fesa::ShellStateCandidate prior_shell{}; + prior_shell.rows = prior_shell_rows.rows; + prior_shell.physical_strain_energy = prior_shell_rows.physical_strain_energy; + prior_shell.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + prior_shell.verification_metrics = {1.0e-11, 2.0e-11, 3.0e-11}; + ASSERT_TRUE(state.CommitShellResults({0U}, std::move(prior_shell)).IsOk()); + const fesa::AnalysisState prior = state; + + fesa::results_internal::RecoveryCandidate candidate{ + fixture.dofs->FullDofCount()}; + for (std::size_t full_dof = 0U; full_dof < fixture.dofs->FullDofCount(); + ++full_dof) { + candidate.displacement[full_dof] = 10.0 + static_cast(full_dof); + candidate.external_force[full_dof] = 20.0 + static_cast(full_dof); + candidate.internal_force[full_dof] = 30.0 + static_cast(full_dof); + candidate.residual[full_dof] = 40.0 + static_cast(full_dof); + candidate.reaction[full_dof] = 50.0 + static_cast(full_dof); + } + candidate.has_shell_results = true; + candidate.expected_shell_elements = {0U}; + auto invalid_shell_rows = MakeShellRows(0U, 12.5); + candidate.shell.rows = invalid_shell_rows.rows; + candidate.shell.rows[3U].location = fesa::ShellMidsurfaceLocation::kGp1; + candidate.shell.physical_strain_energy = + invalid_shell_rows.physical_strain_energy; + candidate.shell.equilibrium = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + candidate.shell.verification_metrics = {0.0, 0.0, 0.0}; + + const auto status = fesa::results_internal::CommitRecoveryCandidate( + std::move(candidate), state); + + ExpectStatusCode(status, "invalid-shell-state-inventory"); + ExpectVectorEqual(state.Displacement(), prior.Displacement()); + ExpectVectorEqual(state.ExternalForce(), prior.ExternalForce()); + ExpectVectorEqual(state.InternalForce(), prior.InternalForce()); + ExpectVectorEqual(state.Residual(), prior.Residual()); + ExpectVectorEqual(state.Reaction(), prior.Reaction()); + ASSERT_EQ(state.EndpointResults().size(), prior.EndpointResults().size()); + EXPECT_EQ(state.EndpointResults()[0U].end_action, + prior.EndpointResults()[0U].end_action); + ASSERT_EQ(state.ShellResults().size(), prior.ShellResults().size()); + EXPECT_EQ(state.ShellResults()[0U].location, + prior.ShellResults()[0U].location); + EXPECT_EQ(state.ShellResults()[0U].section_resultant, + prior.ShellResults()[0U].section_resultant); + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), prior.PhysicalStrainEnergy()); + EXPECT_EQ(state.Equilibrium(), prior.Equilibrium()); + EXPECT_EQ(state.VerificationMetrics(), prior.VerificationMetrics()); +} + +// C-MODULE-002 +TEST(ResultRecoveryComponents, + CommitRecoveryCandidateRejectsNonfiniteVectorAndPreservesWholeState) { + const auto fixture = MakeShellFixture(MakeShellDefinition()); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Prior-Step", 3U}); + for (std::size_t full_dof = 0U; full_dof < fixture.dofs->FullDofCount(); + ++full_dof) { + state.Displacement()[full_dof] = static_cast(full_dof) + 0.25; + state.ExternalForce()[full_dof] = -static_cast(full_dof) - 0.5; + state.InternalForce()[full_dof] = 100.0 + static_cast(full_dof); + state.Residual()[full_dof] = 200.0 + static_cast(full_dof); + state.Reaction()[full_dof] = 300.0 + static_cast(full_dof); + } + const auto prior_beam = MakeBeamRows(0U, 9.0); + state.EndpointResults() = prior_beam.endpoint_rows; + state.GaussResults() = prior_beam.gauss_rows; + state.StressResults() = prior_beam.stress_rows; + auto prior_shell_rows = MakeShellRows(0U, 71.0); + fesa::ShellStateCandidate prior_shell{}; + prior_shell.rows = prior_shell_rows.rows; + prior_shell.physical_strain_energy = prior_shell_rows.physical_strain_energy; + prior_shell.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + prior_shell.verification_metrics = {1.0e-11, 2.0e-11, 3.0e-11}; + ASSERT_TRUE(state.CommitShellResults({0U}, std::move(prior_shell)).IsOk()); + const fesa::AnalysisState prior = state; + + fesa::results_internal::RecoveryCandidate candidate{ + fixture.dofs->FullDofCount()}; + candidate.displacement[0U] = 10.0; + candidate.external_force[0U] = 20.0; + candidate.internal_force[0U] = 30.0; + candidate.residual[0U] = 40.0; + candidate.reaction[0U] = std::numeric_limits::quiet_NaN(); + + const auto status = fesa::results_internal::CommitRecoveryCandidate( + std::move(candidate), state); + + ExpectStatusCode(status, "nonfinite-recovery-value"); + ExpectAnalysisStateEqual(state, prior); +}