feat(linear-static-mitc4-shell): step 10 - shell-result-recovery

include/fesa/results/result_recovery.hpp

@@ -17,7 +17,7 @@ struct NodeStationResultRow {
     std::array<double, 4> sectionResultant;
 };
 
-// Recovers full-space equilibrium and the V0 beam output rows without
+// Recovers full-space equilibrium and the active concrete element rows without
 // exposing element or sparse-backend details to result consumers.
 class ResultRecovery {
 public:

src/fesa/results/result_recovery.cpp

@@ -1,6 +1,7 @@
 #include "fesa/results/result_recovery.hpp"
 
 #include "fesa/elements/euler_beam_3d.hpp"
+#include "fesa/elements/mitc4_shell.hpp"
 
 #include <algorithm>
 #include <array>
@@ -21,7 +22,10 @@ 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;
 
 using AxisSet = std::array<std::array<double, 3>, 3>;
@@ -56,13 +60,8 @@ bool sameSourceIdentity(const SourceEntityId& left,
            left.sourceLabelText == right.sourceLabelText;
 }
 
-bool finite(const std::array<double, 4>& values) {
+template<std::size_t Size>
-    return std::all_of(
+bool finite(const std::array<double, Size>& values) {
-        values.begin(), values.end(),
-        [](const double value) { return std::isfinite(value); });
-}
-
-bool finite(const std::array<double, 6>& values) {
     return std::all_of(
         values.begin(), values.end(),
         [](const double value) { return std::isfinite(value); });
@@ -252,6 +251,74 @@ Status validateRecoveryInputs(const AnalysisModel& model,
                 "Every active element requires one twelve-DOF scatter map.");
         }
     }
+
+    if (!model.activeElements().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::size_t>(
+            (std::numeric_limits<EntityIndex>::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 elementOrder = 0U;
+         elementOrder < domain.shellElements().size();
+         ++elementOrder) {
+        const auto& definition = domain.shellElements()[elementOrder];
+        if (definition.materialIndex >= domain.materials().size() ||
+            definition.sectionIndex >= domain.shellSections().size()) {
+            return recoveryFailure(
+                "invalid-recovery-entity",
+                definition.location,
+                definition.sourceId.sourceLabelText,
+                "Active shell material and section references must resolve before recovery.");
+        }
+        try {
+            const auto& scatter = dofs.shellElementScatter(
+                static_cast<EntityIndex>(elementOrder));
+            for (std::size_t nodePosition = 0U;
+                 nodePosition < definition.nodeIndices.size();
+                 ++nodePosition) {
+                const EntityIndex node = definition.nodeIndices[nodePosition];
+                if (node >= domain.nodes().size()) {
+                    return recoveryFailure(
+                        "invalid-recovery-entity",
+                        definition.location,
+                        definition.sourceId.sourceLabelText,
+                        "Active shell node references must resolve before recovery.");
+                }
+                for (std::size_t component = 0U;
+                     component < kDofsPerNode;
+                     ++component) {
+                    const std::size_t local =
+                        nodePosition * kDofsPerNode + component;
+                    const std::size_t expected =
+                        static_cast<std::size_t>(node) * kDofsPerNode +
+                        component;
+                    if (scatter[local] != expected || expected >= fullCount) {
+                        return recoveryFailure(
+                            "invalid-recovery-order",
+                            definition.location,
+                            definition.sourceId.sourceLabelText,
+                            "Shell scatter must preserve node/component full-DOF order.");
+                    }
+                }
+            }
+        } catch (const std::out_of_range&) {
+            return recoveryFailure(
+                "invalid-recovery-entity",
+                definition.location,
+                definition.sourceId.sourceLabelText,
+                "Every active shell requires one twenty-four-DOF scatter map.");
+        }
+    }
     return Status::ok();
 }
 
@@ -346,6 +413,144 @@ double norm(const std::array<double, 3>& value) {
     return std::hypot(value[0U], value[1U], value[2U]);
 }
 
+bool accumulateVectorAndScale(
+    std::array<double, 3>& total,
+    double& scale,
+    const std::array<double, 3>& contribution) {
+    const double magnitude = norm(contribution);
+    const double accumulatedScale = scale + magnitude;
+    if (!finite(contribution) || !std::isfinite(magnitude) ||
+        !std::isfinite(accumulatedScale)) {
+        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 = accumulatedScale;
+    return true;
+}
+
+double normalizedBalance(
+    const std::array<double, 3>& balance,
+    const double scale) {
+    const double balanceNorm = norm(balance);
+    if (!std::isfinite(balanceNorm) || !std::isfinite(scale)) {
+        return (std::numeric_limits<double>::infinity)();
+    }
+    if (scale == 0.0) {
+        return balanceNorm == 0.0
+            ? 0.0
+            : (std::numeric_limits<double>::infinity)();
+    }
+    return balanceNorm / scale;
+}
+
+Status populateShellGlobalEvidence(
+    const Domain& domain,
+    const DofManager& dofs,
+    const Vector& externalForce,
+    const Vector& residual,
+    const double normalizedResidual,
+    ShellStateCandidate& candidate) {
+    std::vector<unsigned char> constrained(dofs.fullDofCount(), 0U);
+    for (const std::size_t fullDof : dofs.constrainedDofs()) {
+        constrained[fullDof] = 1U;
+    }
+
+    std::array<double, 3> appliedForce{};
+    std::array<double, 3> reactionForce{};
+    std::array<double, 3> appliedMoment{};
+    std::array<double, 3> reactionMoment{};
+    double appliedForceScale = 0.0;
+    double reactionForceScale = 0.0;
+    double appliedMomentScale = 0.0;
+    double reactionMomentScale = 0.0;
+    for (std::size_t node = 0U; node < domain.nodes().size(); ++node) {
+        const std::size_t offset = node * kDofsPerNode;
+        std::array<double, 3> nodalAppliedForce{};
+        std::array<double, 3> nodalReactionForce{};
+        std::array<double, 3> nodalAppliedMoment{};
+        std::array<double, 3> nodalReactionMoment{};
+        for (std::size_t component = 0U; component < 3U; ++component) {
+            nodalAppliedForce[component] = externalForce[offset + component];
+            nodalAppliedMoment[component] =
+                externalForce[offset + 3U + component];
+            if (constrained[offset + component] != 0U) {
+                nodalReactionForce[component] = residual[offset + component];
+            }
+            if (constrained[offset + 3U + component] != 0U) {
+                nodalReactionMoment[component] =
+                    residual[offset + 3U + component];
+            }
+        }
+
+        const auto appliedForceMoment =
+            cross(domain.nodes()[node].coordinates, nodalAppliedForce);
+        const auto reactionForceMoment =
+            cross(domain.nodes()[node].coordinates, nodalReactionForce);
+        for (std::size_t component = 0U; component < 3U; ++component) {
+            nodalAppliedMoment[component] += appliedForceMoment[component];
+            nodalReactionMoment[component] += reactionForceMoment[component];
+        }
+        if (!accumulateVectorAndScale(
+                appliedForce, appliedForceScale, nodalAppliedForce) ||
+            !accumulateVectorAndScale(
+                reactionForce, reactionForceScale, nodalReactionForce) ||
+            !accumulateVectorAndScale(
+                appliedMoment, appliedMomentScale, nodalAppliedMoment) ||
+            !accumulateVectorAndScale(
+                reactionMoment, reactionMomentScale, nodalReactionMoment)) {
+            return recoveryFailure(
+                "nonfinite-recovery-value",
+                domain.nodes()[node].location,
+                domain.nodes()[node].sourceId.sourceLabelText,
+                "Global force and moment evidence must remain finite in source-node order.");
+        }
+    }
+
+    std::array<double, 3> forceBalance{};
+    std::array<double, 3> momentBalance{};
+    for (std::size_t component = 0U; component < 3U; ++component) {
+        forceBalance[component] =
+            appliedForce[component] + reactionForce[component];
+        momentBalance[component] =
+            appliedMoment[component] + reactionMoment[component];
+        candidate.equilibrium[component] = forceBalance[component];
+        candidate.equilibrium[3U + component] = momentBalance[component];
+    }
+    const double forceMetric = normalizedBalance(
+        forceBalance,
+        (std::max)(appliedForceScale, reactionForceScale));
+    const double momentMetric = normalizedBalance(
+        momentBalance,
+        (std::max)(appliedMomentScale, reactionMomentScale));
+    candidate.verificationMetrics = {
+        normalizedResidual, forceMetric, momentMetric};
+    if (!finite(candidate.equilibrium) ||
+        !finite(candidate.verificationMetrics)) {
+        return recoveryFailure(
+            "nonfinite-recovery-value",
+            {domain.sourcePath(), 0U},
+            "global-equilibrium",
+            "Global equilibrium values and their physical normalization scales must be finite.");
+    }
+    if (forceMetric > kGlobalEquilibriumTolerance ||
+        momentMetric > 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<AxisSet> localAxes(const Domain& domain,
                                  const EulerBeam3DDefinition& element) {
     const auto& first = domain.nodes()[element.nodeIndices[0U]].coordinates;
@@ -545,14 +750,156 @@ Status ResultRecovery::recover(const AnalysisModel& model,
         }
     }
 
-    // Commit only after all validation and element recovery succeeds so a
+    ShellStateCandidate shellCandidate{};
-    // failed recovery cannot leave a partially updated AnalysisState.
+    std::vector<EntityIndex> expectedShellElements;
-    state.internalForce() = std::move(internalForce);
+    if (!domain.shellElements().empty()) {
-    state.residual() = std::move(residual);
+        if (domain.shellElements().size() >
-    state.reaction() = std::move(reaction);
+            (std::numeric_limits<std::size_t>::max)() /
-    state.endpointResults() = std::move(endpointRows);
+                kShellLocationCount) {
-    state.gaussResults() = std::move(gaussRows);
+            return recoveryFailure(
-    state.stressResults() = std::move(stressRows);
+                "invalid-recovery-dimensions",
+                {domain.sourcePath(), 0U},
+                domain.sourceContentIdentity(),
+                "The shell result-row inventory exceeds the addressable range.");
+        }
+        std::vector<std::optional<std::array<double, 3>>> directorsByNode(
+            domain.nodes().size());
+        for (const auto& frame : domain.shellNodeInitialFrames()) {
+            if (frame.nodeIndex >= directorsByNode.size() ||
+                directorsByNode[frame.nodeIndex].has_value()) {
+                return recoveryFailure(
+                    "invalid-recovery-entity",
+                    {domain.sourcePath(), 0U},
+                    std::to_string(frame.nodeIndex),
+                    "Shell initial directors must map uniquely to model nodes.");
+            }
+            directorsByNode[frame.nodeIndex] = frame.director;
+        }
+
+        shellCandidate.rows.reserve(
+            domain.shellElements().size() * kShellLocationCount);
+        expectedShellElements.reserve(domain.shellElements().size());
+        constexpr std::array<ShellMidsurfaceLocation, kShellLocationCount>
+            locations{
+                ShellMidsurfaceLocation::gp1,
+                ShellMidsurfaceLocation::gp2,
+                ShellMidsurfaceLocation::gp3,
+                ShellMidsurfaceLocation::gp4};
+        constexpr std::array<ShellSectionPosition, 3> positions{
+            ShellSectionPosition::bottom,
+            ShellSectionPosition::middle,
+            ShellSectionPosition::top};
+        constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
+        for (std::size_t elementOrder = 0U;
+             elementOrder < domain.shellElements().size();
+             ++elementOrder) {
+            const EntityIndex elementIndex =
+                static_cast<EntityIndex>(elementOrder);
+            const auto& definition = domain.shellElements()[elementOrder];
+            std::array<const Node*, 4> nodes{};
+            std::array<std::array<double, 3>, 4> directors{};
+            for (std::size_t nodePosition = 0U;
+                 nodePosition < definition.nodeIndices.size();
+                 ++nodePosition) {
+                const EntityIndex node = definition.nodeIndices[nodePosition];
+                if (!directorsByNode[node].has_value()) {
+                    return recoveryFailure(
+                        "invalid-recovery-entity",
+                        definition.location,
+                        definition.sourceId.sourceLabelText,
+                        "Shell recovery requires one initial director per element node.");
+                }
+                nodes[nodePosition] = &domain.nodes()[node];
+                directors[nodePosition] = *directorsByNode[node];
+            }
+
+            auto shell = Mitc4Shell::create(
+                nodes,
+                directors,
+                domain.shellSections()[definition.sectionIndex],
+                domain.materials()[definition.materialIndex]);
+            if (!shell.hasValue()) {
+                return shell.status();
+            }
+            Vector elementDisplacement{kShellElementDofCount};
+            const auto& scatter = dofs.shellElementScatter(elementIndex);
+            for (std::size_t localDof = 0U;
+                 localDof < kShellElementDofCount;
+                 ++localDof) {
+                elementDisplacement[localDof] =
+                    state.displacement()[scatter[localDof]];
+            }
+            auto recovered = shell.value().recoverPhysical(
+                elementDisplacement);
+            if (!recovered.hasValue()) {
+                return recovered.status();
+            }
+            const double accumulatedEnergy =
+                shellCandidate.physicalStrainEnergy +
+                recovered.value().strainEnergy;
+            if (!std::isfinite(accumulatedEnergy)) {
+                return recoveryFailure(
+                    "nonfinite-recovery-value",
+                    definition.location,
+                    definition.sourceId.sourceLabelText,
+                    "Source-order physical shell energy reduction must remain finite.");
+            }
+            shellCandidate.physicalStrainEnergy = accumulatedEnergy;
+            expectedShellElements.push_back(elementIndex);
+
+            for (std::size_t point = 0U;
+                 point < recovered.value().points.size();
+                 ++point) {
+                const auto& physicalPoint = recovered.value().points[point];
+                ShellResultRow row{};
+                row.element = elementIndex;
+                row.location = locations[point];
+                row.naturalCoordinates = physicalPoint.naturalCoordinates;
+                row.localFrame = {
+                    physicalPoint.localFrame.e1,
+                    physicalPoint.localFrame.e2,
+                    physicalPoint.localFrame.e3};
+                row.generalizedStrain = physicalPoint.generalizedStrain;
+                row.sectionResultant = physicalPoint.sectionResultant;
+                for (std::size_t position = 0U;
+                     position < positions.size();
+                     ++position) {
+                    row.stress[position] = {
+                        positions[position],
+                        zeta[position],
+                        physicalPoint.inPlaneStress[position]};
+                }
+                shellCandidate.rows.push_back(std::move(row));
+            }
+        }
+        const Status evidenceStatus = populateShellGlobalEvidence(
+            domain,
+            dofs,
+            state.externalForce(),
+            residual,
+            normalizedResidual,
+            shellCandidate);
+        if (!evidenceStatus.isOk()) {
+            return evidenceStatus;
+        }
+    }
+
+    // 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 candidateState = state;
+    candidateState.internalForce() = std::move(internalForce);
+    candidateState.residual() = std::move(residual);
+    candidateState.reaction() = std::move(reaction);
+    candidateState.endpointResults() = std::move(endpointRows);
+    candidateState.gaussResults() = std::move(gaussRows);
+    candidateState.stressResults() = std::move(stressRows);
+    const Status shellCommitStatus = candidateState.commitShellResults(
+        expectedShellElements, std::move(shellCandidate));
+    if (!shellCommitStatus.isOk()) {
+        return shellCommitStatus;
+    }
+    state = std::move(candidateState);
     return Status::ok();
 }
 

tests/unit/results/result_recovery_test.cpp

@@ -2,6 +2,7 @@
 
 #include "fesa/analysis/analysis_model.hpp"
 #include "fesa/analysis/analysis_state.hpp"
+#include "fesa/assembly/load_assembler.hpp"
 #include "fesa/assembly/parallel_for.hpp"
 #include "fesa/assembly/sparse_assembler.hpp"
 #include "fesa/fem/dof_manager.hpp"
@@ -34,6 +35,13 @@ struct RecoveryFixture {
     std::unique_ptr<fesa::SparseMatrix> stiffness;
 };
 
+struct ShellRecoveryFixture {
+    std::unique_ptr<fesa::Domain> domain;
+    std::unique_ptr<fesa::AnalysisModel> model;
+    std::unique_ptr<fesa::DofManager> dofs;
+    std::unique_ptr<fesa::SparseMatrix> stiffness;
+};
+
 fesa::ModelDefinition makeDefinition(
     const bool twoElements = false,
     std::vector<std::array<double, 2>> sectionPoints = {},
@@ -144,6 +152,148 @@ RecoveryFixture makeFixture(
         std::move(stiffness)};
 }
 
+fesa::ModelDefinition makeShellDefinition(
+    const bool twoElements = false,
+    const bool constrainAll = false,
+    std::vector<fesa::NodalLoad> loads = {}) {
+    const std::filesystem::path source{"models/shell-result-recovery.inp"};
+    fesa::ModelDefinition definition{};
+    definition.sourcePath = source;
+    definition.sourceContentIdentity = "fnv1a64:fedcba9876543210";
+    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"}, {3.0, -1.0, 0.0}, {source, 12U}},
+        {{"Shell-1", 4, "4"}, {-1.0, 1.0, 0.0}, {source, 13U}},
+        {{"Shell-1", 5, "5"}, {1.0, 1.0, 0.0}, {source, 14U}},
+        {{"Shell-1", 6, "6"}, {3.0, 1.0, 0.0}, {source, 15U}}};
+    if (!twoElements) {
+        definition.nodes.erase(
+            definition.nodes.begin() + 2,
+            definition.nodes.begin() + 3);
+        definition.nodes.erase(definition.nodes.begin() + 4);
+    }
+    definition.materials = {
+        {"Material", 120.0, 0.25, {source, 20U}}};
+    definition.shellSections = {
+        {"ShellSection", 2.0, 0U, {source, 30U}}};
+    definition.shellElements = {{
+        {"Shell-1", 10, "10"},
+        fesa::ShellSourceElementType::s4,
+        {0U, 1U, twoElements ? 4U : 3U, twoElements ? 3U : 2U},
+        0U,
+        0U,
+        {source, 40U}}};
+    if (twoElements) {
+        definition.shellElements.push_back({
+            {"Shell-1", 20, "20"},
+            fesa::ShellSourceElementType::s4r,
+            {1U, 2U, 5U, 4U},
+            0U,
+            0U,
+            {source, 41U}});
+    }
+    for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
+        definition.shellNodeInitialFrames.push_back({
+            static_cast<fesa::EntityIndex>(node),
+            {0.0, 0.0, 1.0},
+            {1.0, 0.0, 0.0},
+            {0.0, 1.0, 0.0}});
+    }
+    if (constrainAll) {
+        std::vector<fesa::EntityIndex> allNodes;
+        allNodes.reserve(definition.nodes.size());
+        for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
+            allNodes.push_back(static_cast<fesa::EntityIndex>(node));
+        }
+        definition.nodeSets.push_back(
+            {"All", {}, std::move(allNodes), {source, 50U}});
+    }
+    definition.steps = {{
+        "Step-1",
+        constrainAll
+            ? std::vector<fesa::BoundaryCondition>{
+                  {"All", 1, 6, 0.0, {source, 60U}}}
+            : std::vector<fesa::BoundaryCondition>{},
+        std::move(loads),
+        0.1,
+        1.0,
+        0.01,
+        1.0,
+        {source, 59U}}};
+    return definition;
+}
+
+ShellRecoveryFixture makeShellFixture(fesa::ModelDefinition definition) {
+    auto domainResult = fesa::Domain::create(std::move(definition));
+    if (!domainResult.hasValue()) {
+        throw std::runtime_error{
+            "Shell recovery fixture Domain construction failed."};
+    }
+    auto domain = std::make_unique<fesa::Domain>(
+        std::move(domainResult.value()));
+
+    auto modelResult = fesa::AnalysisModel::create(*domain);
+    if (!modelResult.hasValue()) {
+        throw std::runtime_error{
+            "Shell recovery fixture AnalysisModel construction failed."};
+    }
+    auto model = std::make_unique<fesa::AnalysisModel>(
+        std::move(modelResult.value()));
+
+    auto dofsResult = fesa::DofManager::create(*model);
+    if (!dofsResult.hasValue()) {
+        throw std::runtime_error{
+            "Shell recovery fixture DofManager construction failed."};
+    }
+    auto dofs = std::make_unique<fesa::DofManager>(
+        std::move(dofsResult.value()));
+
+    fesa::SerialParallelFor serial;
+    auto stiffnessResult = fesa::SparseAssembler::assembleStiffness(
+        *model, *dofs, serial);
+    if (!stiffnessResult.hasValue()) {
+        throw std::runtime_error{
+            "Shell recovery fixture stiffness assembly failed."};
+    }
+    auto stiffness = std::make_unique<fesa::SparseMatrix>(
+        std::move(stiffnessResult.value()));
+    return {
+        std::move(domain),
+        std::move(model),
+        std::move(dofs),
+        std::move(stiffness)};
+}
+
+fesa::AnalysisState makeShellPhysicalState(
+    const ShellRecoveryFixture& fixture) {
+    constexpr std::array<double, 8> 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;
+}
+
 fesa::AnalysisState makeAxialEquilibriumState(const RecoveryFixture& fixture) {
     auto state = fesa::AnalysisState::create(
         *fixture.dofs, {"Step-1", 0U});
@@ -217,6 +367,11 @@ std::vector<fesa::EndpointResultRow> makeStationRows(
 TEST(ResultRecovery, ComputesResidualReactionForNonzeroPrescription) {
     const auto fixture = makeFixture();
     auto state = makeAxialEquilibriumState(fixture);
+    fesa::ShellStateCandidate staleShellEvidence{};
+    staleShellEvidence.physicalStrainEnergy = 123.0;
+    staleShellEvidence.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+    staleShellEvidence.verificationMetrics = {1.0e-11, 2.0e-11, 3.0e-11};
+    ASSERT_TRUE(state.commitShellResults({}, staleShellEvidence).isOk());
 
     const fesa::Status status = fesa::ResultRecovery::recover(
         *fixture.model, *fixture.dofs, *fixture.stiffness, state);
@@ -230,6 +385,14 @@ TEST(ResultRecovery, ComputesResidualReactionForNonzeroPrescription) {
     for (const std::size_t fullDof : fixture.dofs->freeDofs()) {
         EXPECT_DOUBLE_EQ(state.reaction()[fullDof], 0.0);
     }
+    EXPECT_TRUE(state.shellResults().empty());
+    EXPECT_DOUBLE_EQ(state.physicalStrainEnergy(), 0.0);
+    EXPECT_EQ(
+        state.equilibrium(),
+        (std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
+    EXPECT_EQ(
+        state.verificationMetrics(),
+        (std::array<double, 3>{0.0, 0.0, 0.0}));
 }
 
 TEST(ResultRecovery, EnforcesNormalizedFreeResidual) {
@@ -469,3 +632,359 @@ TEST(ResultRecovery, RequiresInteriorEndpointConsistencyWithoutAveraging) {
     ASSERT_FALSE(jumped.hasValue());
     expectStatusCode(jumped.status(), "ineligible-node-station");
 }
+
+// MITC4-REC-001
+TEST(ResultRecovery, RecoversShellRowsInStableElementAndGpOrder) {
+    const auto fixture = makeShellFixture(makeShellDefinition(true));
+    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];
+        state.displacement()[node * 6U] = 0.1 * x + 0.1 * y;
+        state.displacement()[node * 6U + 1U] = -0.05 * y + 0.1 * x;
+    }
+    state.externalForce() =
+        fixture.stiffness->multiply(state.displacement());
+
+    const auto status = fesa::ResultRecovery::recover(
+        *fixture.model, *fixture.dofs, *fixture.stiffness, state);
+
+    ASSERT_TRUE(status.isOk());
+    ASSERT_EQ(state.shellResults().size(), 8U);
+    const double gauss = 1.0 / std::sqrt(3.0);
+    const std::array<fesa::ShellMidsurfaceLocation, 4> locations{
+        fesa::ShellMidsurfaceLocation::gp1,
+        fesa::ShellMidsurfaceLocation::gp2,
+        fesa::ShellMidsurfaceLocation::gp3,
+        fesa::ShellMidsurfaceLocation::gp4};
+    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}};
+    constexpr std::array<double, 8> expectedStrain{
+        0.1, -0.05, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0};
+    constexpr std::array<double, 8> expectedResultant{
+        22.4, -6.4, 19.2, 0.0, 0.0, 0.0, 0.0, 0.0};
+    const std::array<std::array<double, 3>, 3> expectedFrame{{
+        {1.0, 0.0, 0.0},
+        {0.0, 1.0, 0.0},
+        {0.0, 0.0, 1.0}}};
+    for (std::size_t element = 0U; element < 2U; ++element) {
+        for (std::size_t point = 0U; point < locations.size(); ++point) {
+            const auto& row = state.shellResults()[element * 4U + point];
+            EXPECT_EQ(row.element, element);
+            EXPECT_EQ(row.location, locations[point]);
+            EXPECT_EQ(row.naturalCoordinates, coordinates[point]);
+            EXPECT_EQ(row.localFrame, expectedFrame);
+            for (std::size_t component = 0U;
+                 component < expectedStrain.size();
+                 ++component) {
+                EXPECT_NEAR(
+                    row.generalizedStrain[component],
+                    expectedStrain[component],
+                    1.0e-12);
+                EXPECT_NEAR(
+                    row.sectionResultant[component],
+                    expectedResultant[component],
+                    1.0e-12);
+            }
+        }
+    }
+}
+
+// MITC4-REC-002
+TEST(ResultRecovery, RecoversDirectBottomMiddleTopShellStress) {
+    const auto fixture = makeShellFixture(makeShellDefinition());
+    auto state = makeShellPhysicalState(fixture);
+
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *fixture.model,
+                    *fixture.dofs,
+                    *fixture.stiffness,
+                    state)
+                    .isOk());
+    constexpr std::array<fesa::ShellSectionPosition, 3> positions{
+        fesa::ShellSectionPosition::bottom,
+        fesa::ShellSectionPosition::middle,
+        fesa::ShellSectionPosition::top};
+    constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
+    constexpr std::array<std::array<double, 3>, 3> expectedStress{{
+        {-22.4, 6.4, -2.4},
+        {11.2, -3.2, 9.6},
+        {44.8, -12.8, 21.6}}};
+    ASSERT_EQ(state.shellResults().size(), 4U);
+    for (const auto& row : state.shellResults()) {
+        for (std::size_t position = 0U;
+             position < positions.size();
+             ++position) {
+            EXPECT_EQ(row.stress[position].position, positions[position]);
+            EXPECT_DOUBLE_EQ(row.stress[position].zeta, zeta[position]);
+            for (std::size_t component = 0U;
+                 component < expectedStress[position].size();
+                 ++component) {
+                EXPECT_NEAR(
+                    row.stress[position].components[component],
+                    expectedStress[position][component],
+                    1.0e-12);
+            }
+        }
+    }
+}
+
+// MITC4-REC-003
+TEST(ResultRecovery, SumsOnlyPhysicalShellEnergyInSourceOrder) {
+    const auto fixture = makeShellFixture(makeShellDefinition());
+    auto state = makeShellPhysicalState(fixture);
+    constexpr std::array<double, 4> drill{1.0, -1.0, 1.0, -1.0};
+    for (std::size_t node = 0U; node < drill.size(); ++node) {
+        state.displacement()[node * 6U + 5U] = drill[node];
+    }
+    state.externalForce() =
+        fixture.stiffness->multiply(state.displacement());
+    const double stabilizedEnergy = 0.5 * state.displacement().dot(
+        fixture.stiffness->multiply(state.displacement()));
+
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *fixture.model,
+                    *fixture.dofs,
+                    *fixture.stiffness,
+                    state)
+                    .isOk());
+    EXPECT_NEAR(state.physicalStrainEnergy(), 72.16, 1.0e-12);
+    EXPECT_GT(stabilizedEnergy, state.physicalStrainEnergy());
+}
+
+// MITC4-REC-004
+TEST(ResultRecovery, KeepsFullResidualAndComputesGlobalShellEquilibrium) {
+    const std::filesystem::path source{"models/shell-result-recovery.inp"};
+    const std::vector<fesa::NodalLoad> loads{
+        {"1", 1, 5.0, {source, 70U}},
+        {"5", 1, -5.0, {source, 71U}},
+        {"1", 4, 2.0, {source, 72U}},
+        {"5", 4, -2.0, {source, 73U}}};
+    const auto fixture = makeShellFixture(
+        makeShellDefinition(false, true, loads));
+    auto state = fesa::AnalysisState::create(
+        *fixture.dofs, {"Step-1", 0U});
+    auto fullLoad = fesa::LoadAssembler::assembleFullNodalLoad(
+        *fixture.model, *fixture.dofs);
+    ASSERT_TRUE(fullLoad.hasValue());
+    state.externalForce() = std::move(fullLoad.value());
+
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *fixture.model,
+                    *fixture.dofs,
+                    *fixture.stiffness,
+                    state)
+                    .isOk());
+    ASSERT_EQ(state.shellResults().size(), 4U);
+    for (std::size_t fullDof = 0U;
+         fullDof < fixture.dofs->fullDofCount();
+         ++fullDof) {
+        EXPECT_DOUBLE_EQ(state.internalForce()[fullDof], 0.0);
+        EXPECT_DOUBLE_EQ(
+            state.residual()[fullDof], -state.externalForce()[fullDof]);
+        EXPECT_DOUBLE_EQ(state.reaction()[fullDof], state.residual()[fullDof]);
+    }
+    EXPECT_EQ(
+        state.equilibrium(),
+        (std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
+    EXPECT_EQ(
+        state.verificationMetrics(),
+        (std::array<double, 3>{0.0, 0.0, 0.0}));
+
+    const auto freeFixture = makeShellFixture(makeShellDefinition());
+    auto perturbed = makeShellPhysicalState(freeFixture);
+    perturbed.externalForce()[0U] += 1.0e-9;
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *freeFixture.model,
+                    *freeFixture.dofs,
+                    *freeFixture.stiffness,
+                    perturbed)
+                    .isOk());
+    for (const double metric : perturbed.verificationMetrics()) {
+        EXPECT_GT(metric, 0.0);
+        EXPECT_LE(metric, 1.0e-10);
+    }
+}
+
+// MITC4-REC-004
+TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) {
+    auto centeredDefinition = makeShellDefinition();
+    auto translatedDefinition = centeredDefinition;
+    constexpr std::array<double, 3> translation{7.0, 11.0, 0.0};
+    for (auto& node : translatedDefinition.nodes) {
+        for (std::size_t component = 0U;
+             component < translation.size();
+             ++component) {
+            node.coordinates[component] += translation[component];
+        }
+    }
+    const auto centeredFixture = makeShellFixture(
+        std::move(centeredDefinition));
+    const auto translatedFixture = makeShellFixture(
+        std::move(translatedDefinition));
+    auto centered = makeShellPhysicalState(centeredFixture);
+    auto translated = makeShellPhysicalState(translatedFixture);
+    centered.externalForce()[0U] += 1.0e-9;
+    translated.externalForce()[0U] += 1.0e-9;
+
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *centeredFixture.model,
+                    *centeredFixture.dofs,
+                    *centeredFixture.stiffness,
+                    centered)
+                    .isOk());
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *translatedFixture.model,
+                    *translatedFixture.dofs,
+                    *translatedFixture.stiffness,
+                    translated)
+                    .isOk());
+    std::array<double, 3> centeredForce{};
+    for (std::size_t component = 0U; component < 3U; ++component) {
+        centeredForce[component] = centered.equilibrium()[component];
+        EXPECT_NEAR(
+            translated.equilibrium()[component],
+            centeredForce[component],
+            1.0e-12);
+    }
+    const std::array<double, 3> translatedMomentDelta{
+        translation[1U] * centeredForce[2U] -
+            translation[2U] * centeredForce[1U],
+        translation[2U] * centeredForce[0U] -
+            translation[0U] * centeredForce[2U],
+        translation[0U] * centeredForce[1U] -
+            translation[1U] * centeredForce[0U]};
+    for (std::size_t component = 0U; component < 3U; ++component) {
+        EXPECT_NEAR(
+            translated.equilibrium()[3U + component] -
+                centered.equilibrium()[3U + component],
+            translatedMomentDelta[component],
+            5.0e-12);
+    }
+    EXPECT_GT(std::abs(translatedMomentDelta[2U]), 1.0e-9);
+}
+
+// MITC4-REC-004
+TEST(ResultRecovery, UsesScaleAwareShellMetricsAndRejectsExcess) {
+    const auto fixture = makeShellFixture(makeShellDefinition());
+    auto subunit = makeShellPhysicalState(fixture);
+    auto large = makeShellPhysicalState(fixture);
+    constexpr double subunitScale = 1.0e-6;
+    constexpr double largeScale = 1.0e6;
+    subunit.displacement().scale(subunitScale);
+    subunit.externalForce().scale(subunitScale);
+    subunit.externalForce()[0U] += 1.0e-9 * subunitScale;
+    large.displacement().scale(largeScale);
+    large.externalForce().scale(largeScale);
+    large.externalForce()[0U] += 1.0e-9 * largeScale;
+
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *fixture.model,
+                    *fixture.dofs,
+                    *fixture.stiffness,
+                    subunit)
+                    .isOk());
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *fixture.model,
+                    *fixture.dofs,
+                    *fixture.stiffness,
+                    large)
+                    .isOk());
+    for (std::size_t metric = 0U; metric < 3U; ++metric) {
+        EXPECT_GT(subunit.verificationMetrics()[metric], 0.0);
+        EXPECT_GT(large.verificationMetrics()[metric], 0.0);
+        EXPECT_NEAR(
+            subunit.verificationMetrics()[metric],
+            large.verificationMetrics()[metric],
+            1.0e-13);
+    }
+
+    auto constrainedDefinition = makeShellDefinition(false, true);
+    constrainedDefinition.steps[0U].boundaries[0U].value = 1.0;
+    const auto constrainedFixture = makeShellFixture(
+        std::move(constrainedDefinition));
+    auto rejected = fesa::AnalysisState::create(
+        *constrainedFixture.dofs, {"Step-1", 0U});
+    for (std::size_t fullDof = 0U;
+         fullDof < constrainedFixture.dofs->fullDofCount();
+         ++fullDof) {
+        rejected.displacement()[fullDof] = 1.0;
+    }
+    const std::vector<fesa::CooContribution> unbalancedEntry{
+        {0U, 0U, 1.0, 0U, 0U}};
+    auto unbalancedStiffness = fesa::SparseMatrix::fromCoo(
+        constrainedFixture.dofs->fullDofCount(),
+        constrainedFixture.dofs->fullDofCount(),
+        unbalancedEntry,
+        constrainedFixture.dofs->sparsePattern());
+    ASSERT_TRUE(unbalancedStiffness.hasValue());
+    expectStatusCode(
+        fesa::ResultRecovery::recover(
+            *constrainedFixture.model,
+            *constrainedFixture.dofs,
+            unbalancedStiffness.value(),
+            rejected),
+        "global-equilibrium-tolerance-failure");
+}
+
+// MITC4-REC-005
+TEST(ResultRecovery, InvalidLaterShellLeavesEntirePriorStateUnchanged) {
+    const auto validFixture = makeShellFixture(makeShellDefinition(true));
+    auto state = makeShellPhysicalState(validFixture);
+    ASSERT_TRUE(fesa::ResultRecovery::recover(
+                    *validFixture.model,
+                    *validFixture.dofs,
+                    *validFixture.stiffness,
+                    state)
+                    .isOk());
+    ASSERT_EQ(state.shellResults().size(), 8U);
+    const auto priorFirstRow = state.shellResults().front();
+    const double priorEnergy = state.physicalStrainEnergy();
+    const auto priorEquilibrium = state.equilibrium();
+    const auto priorMetrics = state.verificationMetrics();
+    state.internalForce()[0U] = 91.0;
+    state.residual()[0U] = 92.0;
+    state.reaction()[0U] = 93.0;
+    state.endpointResults().push_back({});
+    state.displacement()[2U * 6U] =
+        (std::numeric_limits<double>::max)();
+    state.displacement()[5U * 6U] =
+        (std::numeric_limits<double>::max)();
+    state.externalForce() = fesa::Vector{validFixture.dofs->fullDofCount()};
+    auto zeroStiffness = fesa::SparseMatrix::fromCoo(
+        validFixture.dofs->fullDofCount(),
+        validFixture.dofs->fullDofCount(),
+        {},
+        validFixture.dofs->sparsePattern());
+    ASSERT_TRUE(zeroStiffness.hasValue());
+
+    const auto status = fesa::ResultRecovery::recover(
+        *validFixture.model,
+        *validFixture.dofs,
+        zeroStiffness.value(),
+        state);
+
+    expectStatusCode(status, "invalid-shell-recovery");
+    EXPECT_DOUBLE_EQ(state.internalForce()[0U], 91.0);
+    EXPECT_DOUBLE_EQ(state.residual()[0U], 92.0);
+    EXPECT_DOUBLE_EQ(state.reaction()[0U], 93.0);
+    EXPECT_EQ(state.endpointResults().size(), 1U);
+    ASSERT_EQ(state.shellResults().size(), 8U);
+    EXPECT_EQ(state.shellResults().front().element, priorFirstRow.element);
+    EXPECT_EQ(state.shellResults().front().location, priorFirstRow.location);
+    EXPECT_EQ(
+        state.shellResults().front().generalizedStrain,
+        priorFirstRow.generalizedStrain);
+    EXPECT_EQ(
+        state.shellResults().front().sectionResultant,
+        priorFirstRow.sectionResultant);
+    EXPECT_DOUBLE_EQ(state.physicalStrainEnergy(), priorEnergy);
+    EXPECT_EQ(state.equilibrium(), priorEquilibrium);
+    EXPECT_EQ(state.verificationMetrics(), priorMetrics);
+}