feat(cpp-object-oriented-modular-refactoring): step 17 - generic-result-recovery

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<double, 4> 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,

src/fesa/results/result_recovery.cpp

@@ -4,15 +4,17 @@
 #include <array>
 #include <cmath>
 #include <cstddef>
+#include <functional>
 #include <limits>
+#include <memory>
 #include <optional>
 #include <stdexcept>
 #include <string>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
-#include "fesa/elements/euler_beam_3d.h"
+#include "fesa/elements/element_factory.h"
-#include "fesa/elements/mitc4_shell.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<std::size_t>& 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] !=
+          full_displacement[full_dof] != dofs.PrescribedValues()[constrained]) {
-              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;
+  if (elements.size() != model.ActiveElements().size()) {
-  bool first_element = true;
+    return RecoveryFailure(
-  for (const EntityIndex element : model.ActiveBeamElements()) {
+        "invalid-recovery-entity", {domain.SourcePath(), 0U},
-    if (element >= domain.BeamElements().Size() ||
+        std::to_string(elements.size()),
-        (!first_element && element <= previous_element)) {
+        "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;
+    first_definition = false;
-    previous_element = element;
+    previous_definition = definition_index;
-    const auto& definition = domain.BeamElements()[element];
+    const auto& definition = domain.Elements()[definition_index];
-    if (definition.node_indices[0U] >= domain.Nodes().size() ||
+    const auto& layout = elements[source_order].get().DofLayout();
-        definition.node_indices[1U] >= domain.Nodes().size() ||
+    if (!SameSourceIdentity(layout.source_id, definition.SourceId()) ||
-        definition.material_index >= domain.LinearElasticMaterials().Size() ||
+        layout.node_indices != definition.NodeIndices()) {
-        definition.section_index >= domain.Sections().Size()) {
       return RecoveryFailure(
-          "invalid-recovery-entity", definition.location,
+          "invalid-recovery-entity", {domain.SourcePath(), 0U},
-          definition.source_id.source_label_text,
+          definition.SourceId().source_label_text,
-          "Active beam references must resolve before recovery.");
+          "Runtime recovery layouts must preserve active source identity and "
+          "topology order.");
     }
-    try {
+    auto scatter = dofs.ElementScatter(layout);
-      const auto& scatter = dofs.ElementScatter(element);
+    if (!scatter.HasValue()) {
-      for (std::size_t endpoint = 0U; endpoint < 2U; ++endpoint) {
+      return scatter.GetStatus();
-        for (std::size_t component = 0U; component < kDofsPerNode;
-             ++component) {
-          const std::size_t expected =
-              static_cast<std::size_t>(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.");
     }
   }
+  return Status::Ok();
+}
 
-  if (!model.ActiveBeamElements().empty() && !domain.ShellElements().Empty()) {
+/// @brief Appends one typed beam bundle without changing row identity or sign.
-    return RecoveryFailure(
+Status AppendBeamRows(const SourceLocation& location,
-        "unsupported-mixed-element-model", {domain.SourcePath(), 0U},
+                      const SourceEntityId& source_id,
-        "B33:FESA-MITC4",
+                      const EntityIndex expected_element,
-        "Result recovery does not support mixed beam and shell models.");
+                      const BeamElementResultRows& rows,
-  }
+                      std::vector<EndpointResultRow>& endpoint_rows,
-  if (domain.ShellElements().Size() >
+                      std::vector<GaussResultRow>& gauss_rows,
-      static_cast<std::size_t>((std::numeric_limits<EntityIndex>::max)())) {
+                      std::vector<StressS11Row>& stress_rows) {
-    return RecoveryFailure("invalid-recovery-dimensions",
+  for (const auto& row : rows.endpoint_rows) {
-                           {domain.SourcePath(), 0U},
+    if (row.element != expected_element) {
-                           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<EntityIndex>(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<std::size_t>(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&) {
       return RecoveryFailure(
-          "invalid-recovery-entity", definition.location,
+          "invalid-element-result-identity", location,
-          definition.source_id.source_label_text,
+          source_id.source_label_text,
-          "Every active shell requires one twenty-four-DOF scatter map.");
+          "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<EntityIndex>& 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<AxisSet> LocalAxes(const Domain& domain,
-                                 const EulerBeam3DDefinition& element) {
+                                 const EntityIndex first_node,
-  const auto& first = domain.Nodes()[element.node_indices[0U]].coordinates;
+                                 const EntityIndex second_node,
-  const auto& second = domain.Nodes()[element.node_indices[1U]].coordinates;
+                                 const std::array<double, 3>& 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] =
+    residual[full_dof] = internal_force[full_dof] - external_force[full_dof];
-        internal_force[full_dof] - state.ExternalForce()[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);
+      FreeEquationInternalTermNorms(full_stiffness, displacement, dofs);
-  const double external_norm =
+  const double external_norm = IndexedNorm(external_force, dofs.FreeDofs());
-      IndexedNorm(state.ExternalForce(), 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<EndpointResultRow> endpoint_rows;
   std::vector<GaussResultRow> gauss_rows;
   std::vector<StressS11Row> 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<int>(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<int>(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<EntityIndex> expected_shell_elements;
-  if (!domain.ShellElements().Empty()) {
+  bool has_beam_results = false;
-    if (domain.ShellElements().Size() >
+  bool has_shell_results = false;
-        (std::numeric_limits<std::size_t>::max)() / kShellLocationCount) {
+  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},
+          "invalid-element-result-identity", {domain.SourcePath(), 0U},
-          domain.SourceContentIdentity(),
+          element.DofLayout().source_id.source_label_text,
-          "The shell result-row inventory exceeds the addressable range.");
+          "A runtime result bundle must retain its element source identity.");
-    }
-    std::vector<std::optional<std::array<double, 3>>> 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;
     }
 
-    shell_candidate.rows.reserve(domain.ShellElements().Size() *
+    const Status aggregation_status = std::visit(
-                                 kShellLocationCount);
+        [&](const auto& payload) -> Status {
-    expected_shell_elements.reserve(domain.ShellElements().Size());
+          using Payload = std::decay_t<decltype(payload)>;
-    constexpr std::array<ShellMidsurfaceLocation, kShellLocationCount>
+          if constexpr (std::is_same_v<Payload, BeamElementResultRows>) {
-        locations{ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
+            has_beam_results = true;
-                  ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
+            return AppendBeamRows(
-    constexpr std::array<ShellSectionPosition, 3> positions{
+                {domain.SourcePath(), 0U}, recovered.Value().source_id,
-        ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
+                element_index, payload, endpoint_rows, gauss_rows, stress_rows);
-        ShellSectionPosition::kTop};
+          } else {
-    constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
+            has_shell_results = true;
-    for (std::size_t element_order = 0U;
+            return AppendShellRows({domain.SourcePath(), 0U},
-         element_order < domain.ShellElements().Size(); ++element_order) {
+                                   recovered.Value().source_id, element_index,
-      const EntityIndex element_index = static_cast<EntityIndex>(element_order);
+                                   payload, expected_shell_elements,
-      const auto& definition = domain.ShellElements()[element_order];
+                                   shell_candidate);
-      std::array<const Node*, 4> nodes{};
+          }
-      std::array<std::array<double, 3>, 4> directors{};
+        },
-      for (std::size_t node_position = 0U;
+        recovered.Value().payload);
-           node_position < definition.node_indices.size(); ++node_position) {
+    if (!aggregation_status.IsOk()) {
-        const EntityIndex node = definition.node_indices[node_position];
+      return aggregation_status;
-        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 evidence_status = PopulateShellGlobalEvidence(
+  }
-        domain, dofs, state.ExternalForce(), residual, normalized_residual,
+  if (has_beam_results && has_shell_results) {
-        shell_candidate);
+    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<std::unique_ptr<Element>> 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<std::vector<NodeStationResultRow>>
 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 first_axes = LocalAxes(
-    const auto second_axes = LocalAxes(domain, second_element);
+        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() ||

tests/unit/results/result_recovery_test.cpp

@@ -7,6 +7,7 @@
 #include <cstddef>
 #include <cstdint>
 #include <filesystem>
+#include <functional>
 #include <limits>
 #include <memory>
 #include <stdexcept>
@@ -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<fesa::SparseMatrix> 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<fesa::ElementStiffnessContribution> ComputeStiffness()
+      const override {
+    const std::size_t local_dof_count =
+        layout_.node_indices.size() * layout_.components_per_node.size();
+    return fesa::Result<fesa::ElementStiffnessContribution>::Success(
+        {layout_, fesa::Matrix{local_dof_count, local_dof_count}});
+  }
+
+  fesa::Result<fesa::ElementResultBundle> Recover(
+      const fesa::Vector& element_displacement) const override {
+    observed_displacement_ = element_displacement;
+    if (fail_recovery_) {
+      return fesa::Result<fesa::ElementResultBundle>::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<fesa::ElementResultBundle>::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<fesa::DofComponent> 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<fesa::ShellMidsurfaceLocation, 4> locations{
+      fesa::ShellMidsurfaceLocation::kGp1, fesa::ShellMidsurfaceLocation::kGp2,
+      fesa::ShellMidsurfaceLocation::kGp3, fesa::ShellMidsurfaceLocation::kGp4};
+  const std::array<std::array<double, 2>, 4> coordinates{
+      std::array<double, 2>{-gauss, -gauss},
+      std::array<double, 2>{gauss, -gauss}, std::array<double, 2>{gauss, gauss},
+      std::array<double, 2>{-gauss, gauss}};
+  const std::array<fesa::ShellSectionPosition, 3> positions{
+      fesa::ShellSectionPosition::kBottom, fesa::ShellSectionPosition::kMiddle,
+      fesa::ShellSectionPosition::kTop};
+  constexpr std::array<double, 3> 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<double>(point + 1U);
+    row.section_resultant[0U] = 10.0 * static_cast<double>(point + 1U);
+    for (std::size_t position = 0U; position < positions.size(); ++position) {
+      row.stress[position] = {
+          positions[position],
+          zeta[position],
+          {static_cast<double>(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<std::array<double, 2>> section_points = {},
@@ -343,6 +480,158 @@ std::vector<fesa::EndpointResultRow> 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<double>(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<double>(full_dof) + 0.25;
+    state.ExternalForce()[full_dof] = -static_cast<double>(full_dof) - 0.5;
+    state.InternalForce()[full_dof] = 100.0 + static_cast<double>(full_dof);
+    state.Residual()[full_dof] = 200.0 + static_cast<double>(full_dof);
+    state.Reaction()[full_dof] = 300.0 + static_cast<double>(full_dof);
+  }
+  auto stale_beam_bundle = MakeFakeBeamBundle(*fixture.domain, 0U, 9.0);
+  auto stale_beam_rows =
+      std::get<fesa::BeamElementResultRows>(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<fesa::ShellElementResultRows>(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);