feat(linear-static-mitc4-shell): step 9 - shell-analysis-state

include/fesa/analysis/analysis_state.hpp

@@ -1,9 +1,11 @@
 #pragma once
 
+#include "fesa/core/status.hpp"
 #include "fesa/fem/dof_manager.hpp"
 #include "fesa/math/vector.hpp"
 #include "fesa/results/result_records.hpp"
 
+#include <array>
 #include <cstddef>
 #include <vector>
 
@@ -32,6 +34,13 @@ public:
     const std::vector<GaussResultRow>& gaussResults() const noexcept;
     std::vector<StressS11Row>& stressResults() noexcept;
     const std::vector<StressS11Row>& stressResults() const noexcept;
+    Status commitShellResults(
+        const std::vector<EntityIndex>& expectedElementOrder,
+        ShellStateCandidate candidate);
+    const std::vector<ShellResultRow>& shellResults() const noexcept;
+    double physicalStrainEnergy() const noexcept;
+    const std::array<double, 6>& equilibrium() const noexcept;
+    const std::array<double, 3>& verificationMetrics() const noexcept;
 
 private:
     AnalysisState(std::size_t fullDofCount, StepFrameIdentity identity);
@@ -47,6 +56,11 @@ private:
     std::vector<EndpointResultRow> endpointResults_;
     std::vector<GaussResultRow> gaussResults_;
     std::vector<StressS11Row> stressResults_;
+    // Shell recovery is replaced only through validated candidate commit.
+    std::vector<ShellResultRow> shellResults_;
+    double physicalStrainEnergy_{0.0};
+    std::array<double, 6> equilibrium_{};
+    std::array<double, 3> verificationMetrics_{};
 };
 
 } // namespace fesa

include/fesa/results/result_records.hpp

@@ -5,6 +5,7 @@
 #include <array>
 #include <cstddef>
 #include <string>
+#include <vector>
 
 namespace fesa {
 
@@ -38,4 +39,45 @@ struct StressS11Row {
     std::string source;
 };
 
+enum class ShellMidsurfaceLocation {
+    gp1,
+    gp2,
+    gp3,
+    gp4
+};
+
+enum class ShellSectionPosition {
+    bottom,
+    middle,
+    top
+};
+
+struct ShellSectionStressRow {
+    ShellSectionPosition position;
+    double zeta;
+    std::array<double, 3> components;
+};
+
+struct ShellResultRow {
+    EntityIndex element;
+    ShellMidsurfaceLocation location;
+    std::array<double, 2> naturalCoordinates;
+    // Axis rows [e1,e2,e3], global-component columns.
+    std::array<std::array<double, 3>, 3> localFrame;
+    std::array<double, 8> generalizedStrain;
+    std::array<double, 8> sectionResultant;
+    // Fixed BOTTOM, MIDDLE, TOP order; components are [S11,S22,S12].
+    std::array<ShellSectionStressRow, 3> stress;
+};
+
+struct ShellStateCandidate {
+    std::vector<ShellResultRow> rows;
+    double physicalStrainEnergy{0.0};
+    // [FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3].
+    std::array<double, 6> equilibrium{};
+    // [FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED,
+    //  MOMENT_BALANCE_NORMALIZED].
+    std::array<double, 3> verificationMetrics{};
+};
+
 } // namespace fesa

src/fesa/analysis/analysis_state.cpp

@@ -1,8 +1,57 @@
 #include "fesa/analysis/analysis_state.hpp"
 
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <string>
 #include <utility>
 
 namespace fesa {
+namespace {
+
+constexpr std::size_t kShellLocationsPerElement = 4U;
+
+Status shellCandidateFailure(
+    const std::string& code,
+    const std::string& identity,
+    const std::string& message) {
+    return Status::failure(
+        FailureCategory::model,
+        {{Severity::error,
+          code,
+          {},
+          "ANALYSIS_STATE",
+          identity,
+          message}});
+}
+
+template<std::size_t Size>
+bool finite(const std::array<double, Size>& values) {
+    return std::all_of(
+        values.begin(), values.end(),
+        [](const double value) { return std::isfinite(value); });
+}
+
+bool finite(const ShellResultRow& row) {
+    if (!finite(row.naturalCoordinates) ||
+        !finite(row.generalizedStrain) ||
+        !finite(row.sectionResultant)) {
+        return false;
+    }
+    for (const auto& axis : row.localFrame) {
+        if (!finite(axis)) {
+            return false;
+        }
+    }
+    return std::all_of(
+        row.stress.begin(), row.stress.end(),
+        [](const ShellSectionStressRow& stress) {
+            return std::isfinite(stress.zeta) && finite(stress.components);
+        });
+}
+
+} // namespace
 
 AnalysisState AnalysisState::create(
     const DofManager& dofs, StepFrameIdentity identity) {
@@ -77,6 +126,123 @@ const std::vector<StressS11Row>& AnalysisState::stressResults() const noexcept {
     return stressResults_;
 }
 
+Status AnalysisState::commitShellResults(
+    const std::vector<EntityIndex>& expectedElementOrder,
+    ShellStateCandidate candidate) {
+    if (expectedElementOrder.size() >
+        (std::numeric_limits<std::size_t>::max)() /
+            kShellLocationsPerElement) {
+        return shellCandidateFailure(
+            "invalid-shell-state-inventory",
+            identity_.stepName,
+            "The expected shell result inventory is too large.");
+    }
+    const std::size_t expectedRowCount =
+        expectedElementOrder.size() * kShellLocationsPerElement;
+    if (candidate.rows.size() != expectedRowCount) {
+        return shellCandidateFailure(
+            "invalid-shell-state-inventory",
+            identity_.stepName,
+            "Shell results require exactly four rows per expected element.");
+    }
+    if (std::adjacent_find(
+            expectedElementOrder.begin(), expectedElementOrder.end(),
+            [](const EntityIndex left, const EntityIndex right) {
+                return left >= right;
+            }) != expectedElementOrder.end()) {
+        return shellCandidateFailure(
+            "invalid-shell-state-inventory",
+            identity_.stepName,
+            "Expected shell elements must be unique and in stable index order.");
+    }
+
+    const std::array<ShellMidsurfaceLocation, kShellLocationsPerElement>
+        expectedLocations{
+            ShellMidsurfaceLocation::gp1,
+            ShellMidsurfaceLocation::gp2,
+            ShellMidsurfaceLocation::gp3,
+            ShellMidsurfaceLocation::gp4};
+    const double gauss = 1.0 / std::sqrt(3.0);
+    const std::array<std::array<double, 2>, kShellLocationsPerElement>
+        expectedCoordinates{
+            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<ShellSectionPosition, 3> expectedPositions{
+        ShellSectionPosition::bottom,
+        ShellSectionPosition::middle,
+        ShellSectionPosition::top};
+    constexpr std::array<double, 3> expectedZeta{-1.0, 0.0, 1.0};
+    for (std::size_t elementOrder = 0U;
+         elementOrder < expectedElementOrder.size();
+         ++elementOrder) {
+        for (std::size_t point = 0U;
+             point < kShellLocationsPerElement;
+             ++point) {
+            const auto& row = candidate.rows[
+                elementOrder * kShellLocationsPerElement + point];
+            if (row.element != expectedElementOrder[elementOrder] ||
+                row.location != expectedLocations[point] ||
+                row.naturalCoordinates != expectedCoordinates[point]) {
+                return shellCandidateFailure(
+                    "invalid-shell-state-inventory",
+                    std::to_string(row.element),
+                    "Shell rows must preserve element and GP1 through GP4 identity.");
+            }
+            for (std::size_t position = 0U;
+                 position < expectedPositions.size();
+                 ++position) {
+                if (row.stress[position].position !=
+                        expectedPositions[position] ||
+                    row.stress[position].zeta != expectedZeta[position]) {
+                    return shellCandidateFailure(
+                        "invalid-shell-state-inventory",
+                        std::to_string(row.element),
+                        "Shell stress rows require BOTTOM, MIDDLE, TOP identity.");
+                }
+            }
+            if (!finite(row)) {
+                return shellCandidateFailure(
+                    "nonfinite-shell-state-value",
+                    std::to_string(row.element),
+                    "Shell result rows must contain only finite values.");
+            }
+        }
+    }
+
+    if (!std::isfinite(candidate.physicalStrainEnergy) ||
+        !finite(candidate.equilibrium) ||
+        !finite(candidate.verificationMetrics)) {
+        return shellCandidateFailure(
+            "nonfinite-shell-state-value",
+            identity_.stepName,
+            "Shell energy, equilibrium, and normalized metrics must be finite.");
+    }
+
+    shellResults_ = std::move(candidate.rows);
+    physicalStrainEnergy_ = candidate.physicalStrainEnergy;
+    equilibrium_ = candidate.equilibrium;
+    verificationMetrics_ = candidate.verificationMetrics;
+    return Status::ok();
+}
+
+const std::vector<ShellResultRow>& AnalysisState::shellResults() const noexcept {
+    return shellResults_;
+}
+
+double AnalysisState::physicalStrainEnergy() const noexcept {
+    return physicalStrainEnergy_;
+}
+
+const std::array<double, 6>& AnalysisState::equilibrium() const noexcept {
+    return equilibrium_;
+}
+
+const std::array<double, 3>& AnalysisState::verificationMetrics() const noexcept {
+    return verificationMetrics_;
+}
+
 AnalysisState::AnalysisState(
     std::size_t fullDofCount, StepFrameIdentity identity)
     : identity_{std::move(identity)},

tests/unit/results/result_records_test.cpp

@@ -3,8 +3,11 @@
 #include <gtest/gtest.h>
 
 #include <array>
+#include <cmath>
 #include <filesystem>
+#include <limits>
 #include <utility>
+#include <vector>
 
 namespace {
 
@@ -25,6 +28,130 @@ fesa::DofManager makeEmptyDofs() {
     return std::move(dofs.value());
 }
 
+fesa::ShellResultRow makeShellRow(
+    fesa::EntityIndex element,
+    fesa::ShellMidsurfaceLocation location,
+    std::array<double, 2> naturalCoordinates,
+    double seed) {
+    return {
+        element,
+        location,
+        naturalCoordinates,
+        {{{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}},
+        {seed + 1.0,
+         seed + 2.0,
+         seed + 3.0,
+         seed + 4.0,
+         seed + 5.0,
+         seed + 6.0,
+         seed + 7.0,
+         seed + 8.0},
+        {seed + 11.0,
+         seed + 12.0,
+         seed + 13.0,
+         seed + 14.0,
+         seed + 15.0,
+         seed + 16.0,
+         seed + 17.0,
+         seed + 18.0},
+        {{{fesa::ShellSectionPosition::bottom,
+           -1.0,
+           {seed + 21.0, seed + 22.0, seed + 23.0}},
+          {fesa::ShellSectionPosition::middle,
+           0.0,
+           {seed + 24.0, seed + 25.0, seed + 26.0}},
+          {fesa::ShellSectionPosition::top,
+           1.0,
+           {seed + 27.0, seed + 28.0, seed + 29.0}}}}};
+}
+
+fesa::ShellStateCandidate makeShellCandidate(
+    const std::vector<fesa::EntityIndex>& elements) {
+    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}};
+
+    fesa::ShellStateCandidate candidate{};
+    for (const auto element : elements) {
+        for (std::size_t point = 0U; point < locations.size(); ++point) {
+            candidate.rows.push_back(makeShellRow(
+                element,
+                locations[point],
+                coordinates[point],
+                100.0 * static_cast<double>(element) +
+                    10.0 * static_cast<double>(point)));
+        }
+    }
+    candidate.physicalStrainEnergy = 35.5;
+    candidate.equilibrium = {1.0, -2.0, 3.0, -4.0, 5.0, -6.0};
+    candidate.verificationMetrics = {1.0e-11, 2.0e-11, 3.0e-11};
+    return candidate;
+}
+
+void expectShellStateEquals(
+    const fesa::AnalysisState& state,
+    const std::vector<fesa::ShellResultRow>& rows,
+    double physicalStrainEnergy,
+    const std::array<double, 6>& equilibrium,
+    const std::array<double, 3>& verificationMetrics) {
+    ASSERT_EQ(state.shellResults().size(), rows.size());
+    for (std::size_t row = 0U; row < rows.size(); ++row) {
+        EXPECT_EQ(state.shellResults()[row].element, rows[row].element);
+        EXPECT_EQ(state.shellResults()[row].location, rows[row].location);
+        EXPECT_EQ(
+            state.shellResults()[row].naturalCoordinates,
+            rows[row].naturalCoordinates);
+        EXPECT_EQ(state.shellResults()[row].localFrame, rows[row].localFrame);
+        EXPECT_EQ(
+            state.shellResults()[row].generalizedStrain,
+            rows[row].generalizedStrain);
+        EXPECT_EQ(
+            state.shellResults()[row].sectionResultant,
+            rows[row].sectionResultant);
+        for (std::size_t position = 0U;
+             position < rows[row].stress.size();
+             ++position) {
+            EXPECT_EQ(
+                state.shellResults()[row].stress[position].position,
+                rows[row].stress[position].position);
+            EXPECT_DOUBLE_EQ(
+                state.shellResults()[row].stress[position].zeta,
+                rows[row].stress[position].zeta);
+            EXPECT_EQ(
+                state.shellResults()[row].stress[position].components,
+                rows[row].stress[position].components);
+        }
+    }
+    EXPECT_DOUBLE_EQ(state.physicalStrainEnergy(), physicalStrainEnergy);
+    EXPECT_EQ(state.equilibrium(), equilibrium);
+    EXPECT_EQ(state.verificationMetrics(), verificationMetrics);
+}
+
+void expectShellCandidateRejectedWithoutMutation(
+    fesa::AnalysisState& state,
+    const std::vector<fesa::EntityIndex>& expectedElements,
+    const fesa::ShellStateCandidate& candidate,
+    const fesa::ShellStateCandidate& committed) {
+    const auto status = state.commitShellResults(expectedElements, candidate);
+    EXPECT_FALSE(status.isOk());
+    EXPECT_EQ(status.failureCategory(), fesa::FailureCategory::model);
+    expectShellStateEquals(
+        state,
+        committed.rows,
+        committed.physicalStrainEnergy,
+        committed.equilibrium,
+        committed.verificationMetrics);
+    ASSERT_TRUE(state.commitShellResults(expectedElements, committed).isOk());
+}
+
 } // namespace
 
 TEST(AnalysisState, PreservesStepFrameAndStableRowOrder) {
@@ -100,3 +227,111 @@ TEST(AnalysisState, PreservesStepFrameAndStableRowOrder) {
     EXPECT_EQ(constState.stressResults()[0].source, "input");
     EXPECT_EQ(constState.stressResults()[1].sectionPoint, 2U);
 }
+
+// MITC4-STATE-001
+TEST(AnalysisState, OwnsExactShellRowsInStableElementAndLocationOrder) {
+    const auto dofs = makeEmptyDofs();
+    auto state = fesa::AnalysisState::create(dofs, {"Step-1", 0U});
+    const std::vector<fesa::EntityIndex> expectedElements{3U, 7U};
+    auto candidate = makeShellCandidate(expectedElements);
+
+    const auto status = state.commitShellResults(expectedElements, candidate);
+
+    ASSERT_TRUE(status.isOk());
+    const fesa::AnalysisState& constState = state;
+    ASSERT_EQ(constState.shellResults().size(), 8U);
+    EXPECT_EQ(constState.shellResults()[0].element, 3U);
+    EXPECT_EQ(
+        constState.shellResults()[0].location,
+        fesa::ShellMidsurfaceLocation::gp1);
+    EXPECT_EQ(
+        constState.shellResults()[3].location,
+        fesa::ShellMidsurfaceLocation::gp4);
+    EXPECT_EQ(constState.shellResults()[4].element, 7U);
+    EXPECT_EQ(
+        constState.shellResults()[4].location,
+        fesa::ShellMidsurfaceLocation::gp1);
+    EXPECT_EQ(
+        constState.shellResults()[0].naturalCoordinates,
+        (std::array<double, 2>{
+            -1.0 / std::sqrt(3.0), -1.0 / std::sqrt(3.0)}));
+    EXPECT_EQ(
+        constState.shellResults()[2].generalizedStrain,
+        (std::array<double, 8>{
+            321.0, 322.0, 323.0, 324.0,
+            325.0, 326.0, 327.0, 328.0}));
+    EXPECT_EQ(
+        constState.shellResults()[7].sectionResultant,
+        (std::array<double, 8>{
+            741.0, 742.0, 743.0, 744.0,
+            745.0, 746.0, 747.0, 748.0}));
+    EXPECT_EQ(
+        constState.shellResults()[7].stress[0].position,
+        fesa::ShellSectionPosition::bottom);
+    EXPECT_DOUBLE_EQ(constState.shellResults()[7].stress[0].zeta, -1.0);
+    EXPECT_EQ(
+        constState.shellResults()[7].stress[2].components,
+        (std::array<double, 3>{757.0, 758.0, 759.0}));
+}
+
+// MITC4-STATE-002
+TEST(AnalysisState, CommitsFiniteShellGlobalEvidence) {
+    const auto dofs = makeEmptyDofs();
+    auto state = fesa::AnalysisState::create(dofs, {"Step-1", 0U});
+    const std::vector<fesa::EntityIndex> expectedElements{5U};
+    const auto candidate = makeShellCandidate(expectedElements);
+
+    const auto status = state.commitShellResults(expectedElements, candidate);
+
+    ASSERT_TRUE(status.isOk());
+    EXPECT_DOUBLE_EQ(state.physicalStrainEnergy(), 35.5);
+    EXPECT_EQ(
+        state.equilibrium(),
+        (std::array<double, 6>{1.0, -2.0, 3.0, -4.0, 5.0, -6.0}));
+    EXPECT_EQ(
+        state.verificationMetrics(),
+        (std::array<double, 3>{1.0e-11, 2.0e-11, 3.0e-11}));
+}
+
+// MITC4-STATE-003
+TEST(AnalysisState, InvalidShellCandidatesLeavePriorStateUnchanged) {
+    const auto dofs = makeEmptyDofs();
+    auto state = fesa::AnalysisState::create(dofs, {"Step-1", 0U});
+    const std::vector<fesa::EntityIndex> expectedElements{5U};
+    const auto committed = makeShellCandidate(expectedElements);
+    ASSERT_TRUE(state.commitShellResults(expectedElements, committed).isOk());
+
+    auto invalidLocation = makeShellCandidate(expectedElements);
+    invalidLocation.rows[0].location = fesa::ShellMidsurfaceLocation::gp2;
+    expectShellCandidateRejectedWithoutMutation(
+        state, expectedElements, invalidLocation, committed);
+
+    auto nonfinite = makeShellCandidate(expectedElements);
+    nonfinite.rows[2].generalizedStrain[6] =
+        (std::numeric_limits<double>::quiet_NaN)();
+    expectShellCandidateRejectedWithoutMutation(
+        state, expectedElements, nonfinite, committed);
+
+    auto nonfiniteFrame = makeShellCandidate(expectedElements);
+    nonfiniteFrame.rows[1].localFrame[2][0] =
+        (std::numeric_limits<double>::infinity)();
+    expectShellCandidateRejectedWithoutMutation(
+        state, expectedElements, nonfiniteFrame, committed);
+
+    auto invalidSectionPosition = makeShellCandidate(expectedElements);
+    invalidSectionPosition.rows[3].stress[0].position =
+        fesa::ShellSectionPosition::top;
+    expectShellCandidateRejectedWithoutMutation(
+        state, expectedElements, invalidSectionPosition, committed);
+
+    auto nonfiniteGlobalEvidence = makeShellCandidate(expectedElements);
+    nonfiniteGlobalEvidence.verificationMetrics[1] =
+        (std::numeric_limits<double>::infinity)();
+    expectShellCandidateRejectedWithoutMutation(
+        state, expectedElements, nonfiniteGlobalEvidence, committed);
+
+    auto incompleteInventory = makeShellCandidate(expectedElements);
+    incompleteInventory.rows.pop_back();
+    expectShellCandidateRejectedWithoutMutation(
+        state, expectedElements, incompleteInventory, committed);
+}