From d711e6d4fdbc577c39dff8b8c777a152733cac35 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 16 Aug 2026 10:40:49 +0900 Subject: [PATCH] feat(cpp-object-oriented-modular-refactoring): step 17 - generic-result-recovery --- include/fesa/results/result_recovery.h | 22 +- src/fesa/results/result_recovery.cpp | 505 +++++++++----------- tests/unit/results/result_recovery_test.cpp | 289 +++++++++++ 3 files changed, 543 insertions(+), 273 deletions(-) diff --git a/include/fesa/results/result_recovery.h b/include/fesa/results/result_recovery.h index 2cb797e..d9ac6de 100644 --- a/include/fesa/results/result_recovery.h +++ b/include/fesa/results/result_recovery.h @@ -7,6 +7,7 @@ #include "fesa/analysis/analysis_model.h" #include "fesa/analysis/analysis_state.h" #include "fesa/core/status.h" +#include "fesa/elements/element.h" #include "fesa/fem/dof_manager.h" #include "fesa/math/sparse_matrix.h" @@ -19,10 +20,29 @@ struct NodeStationResultRow { std::array section_resultant; }; -/// @brief Recovers full equilibrium and active concrete element rows. +/// @brief Recovers full equilibrium and typed runtime element rows. class ResultRecovery { public: + /// @brief Aggregates runtime element bundles into one atomic state candidate. + /// @param model Non-owning active semantic model view. + /// @param elements Runtime elements in stable active source order. + /// @param dofs Owner of the matching full-space scatter and partition. + /// @param full_stiffness Assembled full-space stiffness matrix. + /// @param full_displacement Reconstructed full-space displacement candidate. + /// @param full_external_force Assembled full-space external-force candidate. + /// @param state Prior state replaced only after every bundle and global + /// evidence validates. + /// @return Success after atomic commit or a structured model failure. + static Status Recover(const AnalysisModel& model, const ElementView& elements, + const DofManager& dofs, + const SparseMatrix& full_stiffness, + const Vector& full_displacement, + const Vector& full_external_force, + AnalysisState& state); + /// @brief Builds and atomically commits a complete recovery candidate. + /// @note This compatibility facade creates runtime elements until the + /// procedure owns their lifetime directly. /// @return Success after full residual K*d-F and all result rows validate. static Status Recover(const AnalysisModel& model, const DofManager& dofs, const SparseMatrix& full_stiffness, diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp index da03a72..5bceee4 100644 --- a/src/fesa/results/result_recovery.cpp +++ b/src/fesa/results/result_recovery.cpp @@ -4,15 +4,17 @@ #include #include #include +#include #include +#include #include #include #include +#include #include #include -#include "fesa/elements/euler_beam_3d.h" -#include "fesa/elements/mitc4_shell.h" +#include "fesa/elements/element_factory.h" #include "fesa/math/vector3.h" #include "fesa/model/source_target_resolver.h" @@ -20,9 +22,6 @@ namespace fesa { namespace { constexpr std::size_t kDofsPerNode = 6U; -constexpr std::size_t kElementDofCount = 12U; -constexpr std::size_t kShellElementDofCount = 24U; -constexpr std::size_t kShellLocationCount = 4U; constexpr double kFreeResidualTolerance = 1.0e-10; constexpr double kGlobalEquilibriumTolerance = 1.0e-10; constexpr double kAxisTolerance = 1.0e-12; @@ -118,8 +117,11 @@ bool StrictlyIncreasing(const std::vector& values) { /// @brief Validates complete recovery inputs before creating candidate results. Status ValidateRecoveryInputs(const AnalysisModel& model, + const ElementView& elements, const DofManager& dofs, const SparseMatrix& full_stiffness, + const Vector& full_displacement, + const Vector& full_external_force, const AnalysisState& state) { const Domain& domain = model.GetDomain(); if (domain.Nodes().size() > @@ -133,6 +135,8 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, if (dofs.FullDofCount() != full_count || full_stiffness.Rows() != full_count || full_stiffness.Columns() != full_count || + full_displacement.Size() != full_count || + full_external_force.Size() != full_count || state.Displacement().Size() != full_count || state.ExternalForce().Size() != full_count || state.InternalForce().Size() != full_count || @@ -148,7 +152,11 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, if (!matrix_status.IsOk()) { return matrix_status; } - if (!IsFinite(state.Displacement()) || !IsFinite(state.ExternalForce())) { + const Status dof_status = dofs.ValidateInvariants(); + if (!dof_status.IsOk()) { + return dof_status; + } + if (!IsFinite(full_displacement) || !IsFinite(full_external_force)) { return RecoveryFailure( "nonfinite-recovery-value", {domain.SourcePath(), 0U}, domain.SourceContentIdentity(), @@ -191,8 +199,7 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, "Constrained DOFs must be unique and absent from free equations."); } if (!std::isfinite(dofs.PrescribedValues()[constrained]) || - state.Displacement()[full_dof] != - dofs.PrescribedValues()[constrained]) { + full_displacement[full_dof] != dofs.PrescribedValues()[constrained]) { return RecoveryFailure("invalid-recovery-state", {domain.SourcePath(), 0U}, std::to_string(full_dof), @@ -214,111 +221,125 @@ Status ValidateRecoveryInputs(const AnalysisModel& model, "Free and constrained DOFs must partition the full range."); } - EntityIndex previous_element = 0U; - bool first_element = true; - for (const EntityIndex element : model.ActiveBeamElements()) { - if (element >= domain.BeamElements().Size() || - (!first_element && element <= previous_element)) { + if (elements.size() != model.ActiveElements().size()) { + return RecoveryFailure( + "invalid-recovery-entity", {domain.SourcePath(), 0U}, + std::to_string(elements.size()), + "Runtime elements must match the active recovery inventory."); + } + EntityIndex previous_definition = 0U; + bool first_definition = true; + for (std::size_t source_order = 0U; source_order < elements.size(); + ++source_order) { + const EntityIndex definition_index = model.ActiveElements()[source_order]; + if (definition_index >= domain.Elements().Size() || + (!first_definition && definition_index <= previous_definition)) { return RecoveryFailure( "invalid-recovery-entity", {domain.SourcePath(), 0U}, - std::to_string(element), + std::to_string(definition_index), "Active elements must be unique in stable internal-index order."); } - first_element = false; - previous_element = element; - const auto& definition = domain.BeamElements()[element]; - if (definition.node_indices[0U] >= domain.Nodes().size() || - definition.node_indices[1U] >= domain.Nodes().size() || - definition.material_index >= domain.LinearElasticMaterials().Size() || - definition.section_index >= domain.Sections().Size()) { + first_definition = false; + 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()) || + layout.node_indices != definition.NodeIndices()) { return RecoveryFailure( - "invalid-recovery-entity", definition.location, - definition.source_id.source_label_text, - "Active beam references must resolve before recovery."); + "invalid-recovery-entity", {domain.SourcePath(), 0U}, + definition.SourceId().source_label_text, + "Runtime recovery layouts must preserve active source identity and " + "topology order."); } - try { - const auto& scatter = dofs.ElementScatter(element); - for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) { - for (std::size_t component = 0U; component < kDofsPerNode; - ++component) { - const std::size_t expected = - static_cast(definition.node_indices[endpoint]) * - kDofsPerNode + - component; - if (scatter[endpoint * kDofsPerNode + component] != expected || - expected >= full_count) { - return RecoveryFailure("invalid-recovery-order", - definition.location, - definition.source_id.source_label_text, - "Element scatter must preserve " - "endpoint/component full-DOF order."); - } - } - } - } catch (const std::out_of_range&) { - return RecoveryFailure( - "invalid-recovery-entity", definition.location, - definition.source_id.source_label_text, - "Every active element requires one twelve-DOF scatter map."); + auto scatter = dofs.ElementScatter(layout); + if (!scatter.HasValue()) { + return scatter.GetStatus(); } } + return Status::Ok(); +} - if (!model.ActiveBeamElements().empty() && !domain.ShellElements().Empty()) { - return RecoveryFailure( - "unsupported-mixed-element-model", {domain.SourcePath(), 0U}, - "B33:FESA-MITC4", - "Result recovery does not support mixed beam and shell models."); - } - if (domain.ShellElements().Size() > - static_cast((std::numeric_limits::max)())) { - return RecoveryFailure("invalid-recovery-dimensions", - {domain.SourcePath(), 0U}, - domain.SourceContentIdentity(), - "The shell element count cannot be represented by " - "stable element identities."); - } - for (std::size_t element_order = 0U; - element_order < domain.ShellElements().Size(); ++element_order) { - const auto& definition = domain.ShellElements()[element_order]; - if (definition.material_index >= domain.LinearElasticMaterials().Size() || - definition.section_index >= domain.ShellSections().Size()) { - return RecoveryFailure("invalid-recovery-entity", definition.location, - definition.source_id.source_label_text, - "Active shell material and section references " - "must resolve before recovery."); - } - try { - const auto& scatter = - dofs.ShellElementScatter(static_cast(element_order)); - for (std::size_t node_position = 0U; - node_position < definition.node_indices.size(); ++node_position) { - const EntityIndex node = definition.node_indices[node_position]; - if (node >= domain.Nodes().size()) { - return RecoveryFailure( - "invalid-recovery-entity", definition.location, - definition.source_id.source_label_text, - "Active shell node references must resolve before recovery."); - } - for (std::size_t component = 0U; component < kDofsPerNode; - ++component) { - const std::size_t local = node_position * kDofsPerNode + component; - const std::size_t expected = - static_cast(node) * kDofsPerNode + component; - if (scatter[local] != expected || expected >= full_count) { - return RecoveryFailure( - "invalid-recovery-order", definition.location, - definition.source_id.source_label_text, - "Shell scatter must preserve node/component full-DOF order."); - } - } - } - } catch (const std::out_of_range&) { +/// @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-recovery-entity", definition.location, - definition.source_id.source_label_text, - "Every active shell requires one twenty-four-DOF scatter map."); + "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(); } @@ -476,16 +497,18 @@ Status PopulateShellGlobalEvidence(const Domain& domain, const DofManager& dofs, } std::optional LocalAxes(const Domain& domain, - const EulerBeam3DDefinition& element) { - const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates; - const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates; + const EntityIndex first_node, + const EntityIndex second_node, + const std::array& first_axis) { + const auto& first = domain.Nodes()[first_node].coordinates; + const auto& second = domain.Nodes()[second_node].coordinates; const Vector3 delta = Vector3{second} - Vector3{first}; const double length = delta.Norm(); if (!std::isfinite(length) || !(length > 0.0)) { return std::nullopt; } const Vector3 ex = delta / length; - const Vector3 guide{domain.Sections()[element.section_index].first_axis}; + const Vector3 guide{first_axis}; const double projection = guide.Dot(ex); const Vector3 ey_trial = guide - projection * ex; const double ey_norm = ey_trial.Norm(); @@ -519,16 +542,24 @@ bool SameAxes(const AxisSet& left, const AxisSet& right) { } // namespace Status ResultRecovery::Recover(const AnalysisModel& model, + const ElementView& elements, const DofManager& dofs, const SparseMatrix& full_stiffness, + const Vector& full_displacement, + const Vector& full_external_force, AnalysisState& state) { + // Copy aliased AnalysisState inputs before candidate commit can replace the + // caller-owned state. + const Vector displacement = full_displacement; + const Vector external_force = full_external_force; const Status input_status = - ValidateRecoveryInputs(model, dofs, full_stiffness, state); + ValidateRecoveryInputs(model, elements, dofs, full_stiffness, + displacement, external_force, state); if (!input_status.IsOk()) { return input_status; } - Vector internal_force = full_stiffness.Multiply(state.Displacement()); + Vector internal_force = full_stiffness.Multiply(displacement); if (!IsFinite(internal_force)) { return RecoveryFailure( "nonfinite-recovery-value", {model.GetDomain().SourcePath(), 0U}, @@ -537,8 +568,7 @@ Status ResultRecovery::Recover(const AnalysisModel& model, } Vector residual{dofs.FullDofCount()}; for (std::size_t full_dof = 0U; full_dof < residual.Size(); ++full_dof) { - residual[full_dof] = - internal_force[full_dof] - state.ExternalForce()[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}, @@ -549,9 +579,8 @@ Status ResultRecovery::Recover(const AnalysisModel& model, const double residual_norm = IndexedNorm(residual, dofs.FreeDofs()); const auto internal_term_norms = - FreeEquationInternalTermNorms(full_stiffness, state.Displacement(), dofs); - const double external_norm = - IndexedNorm(state.ExternalForce(), dofs.FreeDofs()); + 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. @@ -586,171 +615,68 @@ Status ResultRecovery::Recover(const AnalysisModel& model, std::vector endpoint_rows; std::vector gauss_rows; std::vector stress_rows; - endpoint_rows.reserve(model.ActiveBeamElements().size() * 2U); - gauss_rows.reserve(model.ActiveBeamElements().size() * 2U); const Domain& domain = model.GetDomain(); - for (const EntityIndex element_index : model.ActiveBeamElements()) { - const auto& definition = domain.BeamElements()[element_index]; - auto beam = EulerBeam3D::Create( - domain.Nodes()[definition.node_indices[0U]], - domain.Nodes()[definition.node_indices[1U]], - domain.Sections()[definition.section_index], - domain.LinearElasticMaterials()[definition.material_index]); - if (!beam.HasValue()) { - return beam.GetStatus(); - } - - Vector element_displacement{kElementDofCount}; - const auto& scatter = dofs.ElementScatter(element_index); - for (std::size_t local_dof = 0U; local_dof < kElementDofCount; - ++local_dof) { - element_displacement[local_dof] = - state.Displacement()[scatter[local_dof]]; - } - const BeamRecovery recovered = - beam.Value().RecoverBeam(element_displacement); - for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) { - if (!IsFinite(recovered.equilibrium_end_actions[endpoint]) || - !IsFinite(recovered.endpoint_section_resultants[endpoint])) { - return RecoveryFailure("nonfinite-recovery-value", definition.location, - definition.source_id.source_label_text, - "Endpoint recovery values must be finite."); - } - endpoint_rows.push_back( - {element_index, static_cast(endpoint), - domain.Nodes()[definition.node_indices[endpoint]].source_id, - recovered.equilibrium_end_actions[endpoint], - recovered.endpoint_section_resultants[endpoint]}); - } - for (std::size_t gauss = 0U; gauss < 2U; ++gauss) { - if (!IsFinite(recovered.gauss_generalized_strains[gauss]) || - !IsFinite(recovered.gauss_generalized_resultants[gauss])) { - return RecoveryFailure("nonfinite-recovery-value", definition.location, - definition.source_id.source_label_text, - "Gauss recovery values must be finite."); - } - gauss_rows.push_back({element_index, static_cast(gauss + 1U), - recovered.gauss_generalized_strains[gauss], - recovered.gauss_generalized_resultants[gauss]}); - } - for (const auto& point : recovered.stress_points) { - if ((point.gauss_point != 1 && point.gauss_point != 2) || - !std::isfinite(point.x1) || !std::isfinite(point.x2) || - !std::isfinite(point.s11)) { - return RecoveryFailure( - "nonfinite-recovery-value", definition.location, - definition.source_id.source_label_text, - "Stress recovery identity and values must be finite and ordered."); - } - stress_rows.push_back({element_index, point.gauss_point, - point.section_point, point.x1, point.x2, point.s11, - point.source}); - } - } - ShellStateCandidate shell_candidate{}; std::vector expected_shell_elements; - if (!domain.ShellElements().Empty()) { - if (domain.ShellElements().Size() > - (std::numeric_limits::max)() / kShellLocationCount) { + bool has_beam_results = false; + bool has_shell_results = false; + for (std::size_t source_order = 0U; source_order < elements.size(); + ++source_order) { + const EntityIndex element_index = model.ActiveElements()[source_order]; + const Element& element = elements[source_order].get(); + auto scatter = dofs.ElementScatter(element.DofLayout()); + if (!scatter.HasValue()) { + return scatter.GetStatus(); + } + Vector element_displacement{scatter.Value().size()}; + for (std::size_t local_dof = 0U; local_dof < scatter.Value().size(); + ++local_dof) { + element_displacement[local_dof] = + displacement[scatter.Value()[local_dof]]; + } + auto recovered = element.Recover(element_displacement); + if (!recovered.HasValue()) { + return recovered.GetStatus(); + } + if (!SameSourceIdentity(recovered.Value().source_id, + element.DofLayout().source_id)) { return RecoveryFailure( - "invalid-recovery-dimensions", {domain.SourcePath(), 0U}, - domain.SourceContentIdentity(), - "The shell result-row inventory exceeds the addressable range."); - } - std::vector>> directors_by_node( - domain.Nodes().size()); - for (const auto& frame : domain.ShellNodeInitialFrames()) { - if (frame.node_index >= directors_by_node.size() || - directors_by_node[frame.node_index].has_value()) { - return RecoveryFailure( - "invalid-recovery-entity", {domain.SourcePath(), 0U}, - std::to_string(frame.node_index), - "Shell initial directors must map uniquely to model nodes."); - } - directors_by_node[frame.node_index] = frame.director; + "invalid-element-result-identity", {domain.SourcePath(), 0U}, + element.DofLayout().source_id.source_label_text, + "A runtime result bundle must retain its element source identity."); } - shell_candidate.rows.reserve(domain.ShellElements().Size() * - kShellLocationCount); - expected_shell_elements.reserve(domain.ShellElements().Size()); - constexpr std::array - locations{ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2, - ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4}; - constexpr std::array positions{ - ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle, - ShellSectionPosition::kTop}; - constexpr std::array zeta{-1.0, 0.0, 1.0}; - for (std::size_t element_order = 0U; - element_order < domain.ShellElements().Size(); ++element_order) { - const EntityIndex element_index = static_cast(element_order); - const auto& definition = domain.ShellElements()[element_order]; - std::array nodes{}; - std::array, 4> directors{}; - for (std::size_t node_position = 0U; - node_position < definition.node_indices.size(); ++node_position) { - const EntityIndex node = definition.node_indices[node_position]; - if (!directors_by_node[node].has_value()) { - return RecoveryFailure( - "invalid-recovery-entity", definition.location, - definition.source_id.source_label_text, - "Shell recovery requires one initial director per element node."); - } - nodes[node_position] = &domain.Nodes()[node]; - directors[node_position] = *directors_by_node[node]; - } - - auto shell = Mitc4Shell::Create( - nodes, directors, domain.ShellSections()[definition.section_index], - domain.LinearElasticMaterials()[definition.material_index]); - if (!shell.HasValue()) { - return shell.GetStatus(); - } - Vector element_displacement{kShellElementDofCount}; - const auto& scatter = dofs.ShellElementScatter(element_index); - for (std::size_t local_dof = 0U; local_dof < kShellElementDofCount; - ++local_dof) { - element_displacement[local_dof] = - state.Displacement()[scatter[local_dof]]; - } - auto recovered = shell.Value().RecoverPhysical(element_displacement); - if (!recovered.HasValue()) { - return recovered.GetStatus(); - } - const double accumulated_energy = shell_candidate.physical_strain_energy + - recovered.Value().strain_energy; - if (!std::isfinite(accumulated_energy)) { - return RecoveryFailure( - "nonfinite-recovery-value", definition.location, - definition.source_id.source_label_text, - "Source-order physical shell energy reduction must remain finite."); - } - shell_candidate.physical_strain_energy = accumulated_energy; - expected_shell_elements.push_back(element_index); - - for (std::size_t point = 0U; point < recovered.Value().points.size(); - ++point) { - const auto& physical_point = recovered.Value().points[point]; - ShellResultRow row{}; - row.element = element_index; - row.location = locations[point]; - row.natural_coordinates = physical_point.natural_coordinates; - row.local_frame = {physical_point.local_frame.e1, - physical_point.local_frame.e2, - physical_point.local_frame.e3}; - row.generalized_strain = physical_point.generalized_strain; - row.section_resultant = physical_point.section_resultant; - for (std::size_t position = 0U; position < positions.size(); - ++position) { - row.stress[position] = {positions[position], zeta[position], - physical_point.in_plane_stress[position]}; - } - shell_candidate.rows.push_back(std::move(row)); - } + const Status aggregation_status = std::visit( + [&](const auto& payload) -> Status { + using Payload = std::decay_t; + if constexpr (std::is_same_v) { + has_beam_results = true; + return AppendBeamRows( + {domain.SourcePath(), 0U}, recovered.Value().source_id, + element_index, payload, endpoint_rows, gauss_rows, stress_rows); + } else { + has_shell_results = true; + return AppendShellRows({domain.SourcePath(), 0U}, + recovered.Value().source_id, element_index, + payload, expected_shell_elements, + shell_candidate); + } + }, + recovered.Value().payload); + if (!aggregation_status.IsOk()) { + return aggregation_status; } - const Status evidence_status = PopulateShellGlobalEvidence( - domain, dofs, state.ExternalForce(), residual, normalized_residual, - shell_candidate); + } + if (has_beam_results && 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 (!evidence_status.IsOk()) { return evidence_status; } @@ -760,6 +686,8 @@ Status ResultRecovery::Recover(const AnalysisModel& model, // 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); @@ -775,6 +703,34 @@ Status ResultRecovery::Recover(const AnalysisModel& model, return Status::Ok(); } +Status ResultRecovery::Recover(const AnalysisModel& model, + const DofManager& dofs, + const SparseMatrix& full_stiffness, + AnalysisState& state) { + const Domain& domain = model.GetDomain(); + std::vector> owned_elements; + ElementView elements; + owned_elements.reserve(model.ActiveElements().size()); + elements.reserve(model.ActiveElements().size()); + const ElementFactory factory; + for (const EntityIndex element_index : model.ActiveElements()) { + if (element_index >= domain.Elements().Size()) { + return RecoveryFailure( + "invalid-recovery-entity", {domain.SourcePath(), 0U}, + std::to_string(element_index), + "An active recovery element is outside the Domain."); + } + auto candidate = factory.Create(domain.Elements()[element_index], domain); + if (!candidate.HasValue()) { + return candidate.GetStatus(); + } + owned_elements.push_back(std::move(candidate.Value())); + elements.push_back(std::cref(*owned_elements.back())); + } + return Recover(model, elements, dofs, full_stiffness, state.Displacement(), + state.ExternalForce(), state); +} + Result> ResultRecovery::NormalizeSectionResultantsToNodeStations( const AnalysisModel& model, @@ -893,8 +849,13 @@ ResultRecovery::NormalizeSectionResultantsToNodeStations( incident[0U]->endpoint != incident[1U]->endpoint && ((incident[0U]->endpoint == 1 && incident[1U]->endpoint == 0) || (incident[0U]->endpoint == 0 && incident[1U]->endpoint == 1)); - const auto first_axes = LocalAxes(domain, first_element); - const auto second_axes = LocalAxes(domain, second_element); + const auto first_axes = LocalAxes( + domain, first_element.node_indices[0U], first_element.node_indices[1U], + domain.Sections()[first_element.section_index].first_axis); + const auto second_axes = + LocalAxes(domain, second_element.node_indices[0U], + second_element.node_indices[1U], + domain.Sections()[second_element.section_index].first_axis); if (!chain_orientation || first_element.section_index != second_element.section_index || !first_axes.has_value() || !second_axes.has_value() || diff --git a/tests/unit/results/result_recovery_test.cpp b/tests/unit/results/result_recovery_test.cpp index 8a6b268..baee64d 100644 --- a/tests/unit/results/result_recovery_test.cpp +++ b/tests/unit/results/result_recovery_test.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -19,6 +20,7 @@ #include "fesa/assembly/load_assembler.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" @@ -43,6 +45,141 @@ struct ShellRecoveryFixture { std::unique_ptr stiffness; }; +class FakeRecoveryElement final : public fesa::Element { + public: + FakeRecoveryElement(fesa::ElementDofLayout layout, + fesa::ElementResultBundle bundle, + const bool fail_recovery = false) + : layout_{std::move(layout)}, + bundle_{std::move(bundle)}, + fail_recovery_{fail_recovery} {} + + const fesa::ElementDofLayout& DofLayout() const noexcept override { + return layout_; + } + + fesa::Result ComputeStiffness() + const override { + const std::size_t local_dof_count = + layout_.node_indices.size() * layout_.components_per_node.size(); + return fesa::Result::Success( + {layout_, fesa::Matrix{local_dof_count, local_dof_count}}); + } + + fesa::Result Recover( + const fesa::Vector& element_displacement) const override { + observed_displacement_ = element_displacement; + if (fail_recovery_) { + return fesa::Result::Failure( + fesa::Status::Failure( + fesa::FailureCategory::kModel, + {{fesa::Severity::kError, + "fake-recovery-failure", + {}, + "RESULT_RECOVERY", + layout_.source_id.source_label_text, + "The fake runtime element rejected recovery."}})); + } + return fesa::Result::Success(bundle_); + } + + const fesa::Vector& ObservedDisplacement() const noexcept { + return observed_displacement_; + } + + private: + fesa::ElementDofLayout layout_; + fesa::ElementResultBundle bundle_; + bool fail_recovery_; + mutable fesa::Vector observed_displacement_{0U}; +}; + +std::vector FullNodeComponents() { + return {fesa::DofComponent::kUx, fesa::DofComponent::kUy, + fesa::DofComponent::kUz, fesa::DofComponent::kUrx, + fesa::DofComponent::kUry, fesa::DofComponent::kUrz}; +} + +fesa::ElementDofLayout MakeRuntimeLayout(const fesa::Domain& domain, + const fesa::EntityIndex element) { + const auto& definition = domain.Elements()[element]; + return {definition.SourceId(), definition.NodeIndices(), + FullNodeComponents()}; +} + +fesa::ElementResultBundle MakeFakeBeamBundle(const fesa::Domain& domain, + const fesa::EntityIndex element, + const double value) { + const auto& definition = domain.Elements()[element]; + fesa::BeamElementResultRows rows{}; + rows.endpoint_rows = {{element, + 0, + domain.Nodes()[definition.NodeIndices()[0U]].source_id, + {-value, 0.0, 0.0, 0.0, 0.0, 0.0}, + {value, 2.0 * value, 3.0 * value, 4.0 * value}}, + {element, + 1, + domain.Nodes()[definition.NodeIndices()[1U]].source_id, + {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, 0.0, 0.0, 0.0}, + {value, 2.0 * value, 3.0 * value, 4.0 * value}}, + {element, + 2, + {value + 0.5, 0.0, 0.0, 0.0}, + {value, 2.0 * value, 3.0 * value, 4.0 * value}}}; + rows.stress_rows = {{element, 1, 0U, 0.0, 0.0, 5.0 * value, "fake"}, + {element, 2, 0U, 0.0, 0.0, 6.0 * value, "fake"}}; + return {definition.SourceId(), std::move(rows)}; +} + +fesa::ElementResultBundle MakeFakeShellBundle(const fesa::Domain& domain, + 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}}}; + row.generalized_strain[0U] = static_cast(point + 1U); + row.section_resultant[0U] = 10.0 * static_cast(point + 1U); + for (std::size_t position = 0U; position < positions.size(); ++position) { + row.stress[position] = { + positions[position], + zeta[position], + {static_cast(point + position + 1U), 0.0, 0.0}}; + } + rows.rows.push_back(std::move(row)); + } + return {domain.Elements()[element].SourceId(), std::move(rows)}; +} + +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]); + } +} + fesa::ModelDefinition MakeDefinition( const bool two_elements = false, std::vector> section_points = {}, @@ -343,6 +480,158 @@ std::vector MakeStationRows( } // namespace +// C-RECOVERY-001 +TEST(ResultRecovery, AggregatesFakeBeamBundlesInStableRuntimeOrder) { + const auto fixture = MakeFixture(true); + const auto input = MakeAxialEquilibriumState(fixture); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + + FakeRecoveryElement first{MakeRuntimeLayout(*fixture.domain, 0U), + MakeFakeBeamBundle(*fixture.domain, 0U, 1.0)}; + FakeRecoveryElement second{MakeRuntimeLayout(*fixture.domain, 1U), + MakeFakeBeamBundle(*fixture.domain, 1U, 2.0)}; + const fesa::ElementView elements{std::cref(first), std::cref(second)}; + + const fesa::Status status = fesa::ResultRecovery::Recover( + *fixture.model, elements, *fixture.dofs, *fixture.stiffness, + input.Displacement(), input.ExternalForce(), state); + + ASSERT_TRUE(status.IsOk()); + ExpectVectorEqual(state.Displacement(), input.Displacement()); + ExpectVectorEqual(state.ExternalForce(), input.ExternalForce()); + ASSERT_EQ(first.ObservedDisplacement().Size(), 12U); + EXPECT_DOUBLE_EQ(first.ObservedDisplacement()[0U], 0.1); + EXPECT_DOUBLE_EQ(first.ObservedDisplacement()[6U], 0.3); + ASSERT_EQ(second.ObservedDisplacement().Size(), 12U); + EXPECT_DOUBLE_EQ(second.ObservedDisplacement()[0U], 0.3); + EXPECT_DOUBLE_EQ(second.ObservedDisplacement()[6U], 0.0); + + ASSERT_EQ(state.EndpointResults().size(), 4U); + ASSERT_EQ(state.GaussResults().size(), 4U); + ASSERT_EQ(state.StressResults().size(), 4U); + for (std::size_t element = 0U; element < 2U; ++element) { + const double value = static_cast(element + 1U); + const auto& negative_endpoint = state.EndpointResults()[2U * element]; + const auto& positive_endpoint = state.EndpointResults()[2U * element + 1U]; + EXPECT_EQ(negative_endpoint.element, element); + EXPECT_EQ(negative_endpoint.endpoint, 0); + EXPECT_DOUBLE_EQ(negative_endpoint.end_action[0U], -value); + EXPECT_DOUBLE_EQ(negative_endpoint.section_resultant[0U], value); + EXPECT_EQ(positive_endpoint.element, element); + EXPECT_EQ(positive_endpoint.endpoint, 1); + EXPECT_DOUBLE_EQ(positive_endpoint.end_action[0U], value); + EXPECT_DOUBLE_EQ(positive_endpoint.section_resultant[0U], value); + EXPECT_EQ(state.GaussResults()[2U * element].element, element); + EXPECT_EQ(state.GaussResults()[2U * element].gauss_point, 1); + EXPECT_EQ(state.GaussResults()[2U * element + 1U].element, element); + EXPECT_EQ(state.GaussResults()[2U * element + 1U].gauss_point, 2); + EXPECT_EQ(state.StressResults()[2U * element].element, element); + EXPECT_EQ(state.StressResults()[2U * element + 1U].element, element); + } + EXPECT_TRUE(state.ShellResults().empty()); +} + +// C-RECOVERY-001 +TEST(ResultRecovery, AggregatesFakeShellBundleAndPhysicalEnergy) { + const auto fixture = MakeShellFixture(MakeShellDefinition()); + const auto input = MakeShellPhysicalState(fixture); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + constexpr double physical_energy = 37.5; + FakeRecoveryElement shell{ + MakeRuntimeLayout(*fixture.domain, 0U), + MakeFakeShellBundle(*fixture.domain, 0U, physical_energy)}; + const fesa::ElementView elements{std::cref(shell)}; + + const fesa::Status status = fesa::ResultRecovery::Recover( + *fixture.model, elements, *fixture.dofs, *fixture.stiffness, + input.Displacement(), input.ExternalForce(), state); + + ASSERT_TRUE(status.IsOk()); + ExpectVectorEqual(state.Displacement(), input.Displacement()); + ExpectVectorEqual(state.ExternalForce(), input.ExternalForce()); + EXPECT_TRUE(state.EndpointResults().empty()); + EXPECT_TRUE(state.GaussResults().empty()); + EXPECT_TRUE(state.StressResults().empty()); + ASSERT_EQ(state.ShellResults().size(), 4U); + EXPECT_EQ(state.ShellResults()[0U].location, + fesa::ShellMidsurfaceLocation::kGp1); + EXPECT_EQ(state.ShellResults()[3U].location, + fesa::ShellMidsurfaceLocation::kGp4); + EXPECT_DOUBLE_EQ(state.ShellResults()[0U].generalized_strain[0U], 1.0); + EXPECT_DOUBLE_EQ(state.ShellResults()[3U].section_resultant[0U], 40.0); + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), physical_energy); +} + +// C-RECOVERY-001 +TEST(ResultRecovery, FailingFakeBundleRollsBackTheWholeState) { + const auto fixture = MakeFixture(true); + const auto input = MakeAxialEquilibriumState(fixture); + 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); + } + auto stale_beam_bundle = MakeFakeBeamBundle(*fixture.domain, 0U, 9.0); + auto stale_beam_rows = + std::get(stale_beam_bundle.payload); + state.EndpointResults() = std::move(stale_beam_rows.endpoint_rows); + state.GaussResults() = std::move(stale_beam_rows.gauss_rows); + state.StressResults() = std::move(stale_beam_rows.stress_rows); + auto stale_shell_bundle = MakeFakeShellBundle(*fixture.domain, 0U, 71.0); + const auto& stale_shell_rows = + std::get(stale_shell_bundle.payload); + fesa::ShellStateCandidate stale_shell_candidate{}; + stale_shell_candidate.rows = stale_shell_rows.rows; + stale_shell_candidate.physical_strain_energy = + stale_shell_rows.physical_strain_energy; + stale_shell_candidate.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + stale_shell_candidate.verification_metrics = {1.0e-11, 2.0e-11, 3.0e-11}; + ASSERT_TRUE( + state.CommitShellResults({0U}, std::move(stale_shell_candidate)).IsOk()); + const fesa::AnalysisState prior = state; + + FakeRecoveryElement first{MakeRuntimeLayout(*fixture.domain, 0U), + MakeFakeBeamBundle(*fixture.domain, 0U, 1.0)}; + FakeRecoveryElement failing{MakeRuntimeLayout(*fixture.domain, 1U), + MakeFakeBeamBundle(*fixture.domain, 1U, 2.0), + true}; + const fesa::ElementView elements{std::cref(first), std::cref(failing)}; + + const fesa::Status status = fesa::ResultRecovery::Recover( + *fixture.model, elements, *fixture.dofs, *fixture.stiffness, + input.Displacement(), input.ExternalForce(), state); + + ExpectStatusCode(status, "fake-recovery-failure"); + 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()); + EXPECT_EQ(state.Identity().step_name, prior.Identity().step_name); + EXPECT_EQ(state.Identity().frame_index, prior.Identity().frame_index); + ASSERT_EQ(state.EndpointResults().size(), prior.EndpointResults().size()); + EXPECT_EQ(state.EndpointResults()[0U].element, + prior.EndpointResults()[0U].element); + EXPECT_EQ(state.EndpointResults()[0U].end_action, + prior.EndpointResults()[0U].end_action); + ASSERT_EQ(state.GaussResults().size(), prior.GaussResults().size()); + EXPECT_EQ(state.GaussResults()[0U].generalized_resultant, + prior.GaussResults()[0U].generalized_resultant); + ASSERT_EQ(state.StressResults().size(), prior.StressResults().size()); + EXPECT_DOUBLE_EQ(state.StressResults()[0U].s11, + prior.StressResults()[0U].s11); + ASSERT_EQ(state.ShellResults().size(), prior.ShellResults().size()); + EXPECT_EQ(state.ShellResults()[0U].generalized_strain, + prior.ShellResults()[0U].generalized_strain); + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), prior.PhysicalStrainEnergy()); + EXPECT_EQ(state.Equilibrium(), prior.Equilibrium()); + EXPECT_EQ(state.VerificationMetrics(), prior.VerificationMetrics()); +} + TEST(ResultRecovery, ComputesResidualReactionForNonzeroPrescription) { const auto fixture = MakeFixture(); auto state = MakeAxialEquilibriumState(fixture);