#include "fesa/results/result_recovery.hpp"

#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"
#include "fesa/model/domain.h"

#include <gtest/gtest.h>

#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

namespace {

constexpr double kYoungsModulus = 100.0;
constexpr double kPoissonRatio = 0.25;
constexpr double kLength = 2.0;

struct RecoveryFixture {
    std::unique_ptr<fesa::Domain> domain;
    std::unique_ptr<fesa::AnalysisModel> model;
    std::unique_ptr<fesa::DofManager> dofs;
    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 = {},
    std::vector<fesa::NodalLoad> loads = {},
    const bool reverseSecond = false,
    const bool sectionJump = false,
    const bool nonzeroPrescription = true) {
    const std::filesystem::path source{"models/result-recovery.inp"};
    fesa::ModelDefinition definition{};
    definition.source_path = source;
    definition.source_content_identity = "fnv1a64:0123456789abcdef";
    definition.nodes = {
        {{"Beam-1", 1, "1"}, {0.0, 0.0, 0.0}, {source, 10U}},
        {{"Beam-1", 2, "2"}, {kLength, 0.0, 0.0}, {source, 11U}}};
    if (twoElements) {
        definition.nodes.push_back(
            {{"Beam-1", 3, "3"}, {2.0 * kLength, 0.0, 0.0}, {source, 12U}});
    }
    definition.materials = {
        {"Material", kYoungsModulus, kPoissonRatio, {source, 20U}}};
    definition.sections = {{
        "Section",
        2.0,
        3.0,
        0.0,
        4.0,
        5.0,
        {0.0, 1.0, 0.0},
        std::move(sectionPoints),
        {source, 30U}}};
    if (sectionJump) {
        auto secondSection = definition.sections.front();
        secondSection.name = "Section-2";
        secondSection.area = 2.5;
        secondSection.location.line = 31U;
        definition.sections.push_back(std::move(secondSection));
    }
    definition.elements = {
        {{"Beam-1", 10, "10"}, {0U, 1U}, 0U, 0U, {source, 40U}}};
    if (twoElements) {
        definition.elements.push_back({
            {"Beam-1", 20, "20"},
            reverseSecond ? std::array<fesa::EntityIndex, 2>{2U, 1U}
                          : std::array<fesa::EntityIndex, 2>{1U, 2U},
            0U,
            sectionJump ? 1U : 0U,
            {source, 41U}});
    }
    definition.steps = {{
        "Step-1",
        {{"1", 1, 1, nonzeroPrescription ? 0.1 : 0.0, {source, 50U}},
         {"1", 2, 6, 0.0, {source, 51U}}},
        std::move(loads),
        0.1,
        1.0,
        0.01,
        1.0,
        {source, 49U}}};
    return definition;
}

RecoveryFixture makeFixture(
    const bool twoElements = false,
    std::vector<std::array<double, 2>> sectionPoints = {},
    std::vector<fesa::NodalLoad> loads = {},
    const bool reverseSecond = false,
    const bool sectionJump = false,
    const bool nonzeroPrescription = true) {
    auto domainResult = fesa::Domain::Create(makeDefinition(
        twoElements,
        std::move(sectionPoints),
        std::move(loads),
        reverseSecond,
        sectionJump,
        nonzeroPrescription));
    if (!domainResult.HasValue()) {
        throw std::runtime_error{"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{"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{"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{"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::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.source_path = source;
    definition.source_content_identity = "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.shell_sections = {
        {"ShellSection", 2.0, 0U, {source, 30U}}};
    definition.shell_elements = {{
        {"Shell-1", 10, "10"},
        fesa::ShellSourceElementType::kS4,
        {0U, 1U, twoElements ? 4U : 3U, twoElements ? 3U : 2U},
        0U,
        0U,
        {source, 40U}}};
    if (twoElements) {
        definition.shell_elements.push_back({
            {"Shell-1", 20, "20"},
            fesa::ShellSourceElementType::kS4r,
            {1U, 2U, 5U, 4U},
            0U,
            0U,
            {source, 41U}});
    }
    for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
        definition.shell_node_initial_frames.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.node_sets.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});
    state.displacement()[0U] = 0.1;
    state.displacement()[6U] = 0.3;
    const fesa::Vector internal = fixture.stiffness->Multiply(state.displacement());
    for (const std::size_t fullDof : fixture.dofs->freeDofs()) {
        state.externalForce()[fullDof] = internal[fullDof];
    }
    return state;
}

fesa::AnalysisState makePatchState(
    const RecoveryFixture& fixture,
    const double epsilon,
    const double twist,
    const double kappaY,
    const double kappaZ) {
    auto state = fesa::AnalysisState::create(
        *fixture.dofs, {"Step-1", 0U});
    state.displacement()[0U] = 0.1;
    state.displacement()[6U] = 0.1 + epsilon * kLength;
    state.displacement()[7U] = 0.5 * kappaZ * kLength * kLength;
    state.displacement()[8U] = -0.5 * kappaY * kLength * kLength;
    state.displacement()[9U] = twist * kLength;
    state.displacement()[10U] = kappaY * kLength;
    state.displacement()[11U] = kappaZ * kLength;
    state.externalForce() = fixture.stiffness->Multiply(state.displacement());
    return state;
}

void expectStatusCode(const fesa::Status& status, const std::string& code) {
    ASSERT_FALSE(status.IsOk());
    EXPECT_EQ(status.Category(), fesa::FailureCategory::kModel);
    ASSERT_EQ(status.Diagnostics().size(), 1U);
    EXPECT_EQ(status.Diagnostics()[0U].code, code);
}

void expectScaledNear(
    const double actual,
    const double expected,
    const double relativeTolerance = 1.0e-12) {
    ASSERT_TRUE(std::isfinite(actual));
    ASSERT_TRUE(std::isfinite(expected));
    EXPECT_LE(
        std::abs(actual - expected),
        relativeTolerance * (std::max)(std::abs(expected), 1.0));
}

std::vector<fesa::EndpointResultRow> makeStationRows(
    const RecoveryFixture& fixture) {
    const auto& nodes = fixture.domain->Nodes();
    const auto& elements = fixture.domain->Elements();
    return {
        {0U, 0, nodes[elements[0U].node_indices[0U]].source_id,
         {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
         {1.0, 2.0, 3.0, 4.0}},
        {0U, 1, nodes[elements[0U].node_indices[1U]].source_id,
         {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
         {5.0, 6.0, 7.0, 8.0}},
        {1U, 0, nodes[elements[1U].node_indices[0U]].source_id,
         {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
         {5.0, 6.0, 7.0, 8.0}},
        {1U, 1, nodes[elements[1U].node_indices[1U]].source_id,
         {0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
         {9.0, 10.0, 11.0, 12.0}}};
}

} // namespace

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);

    ASSERT_TRUE(status.IsOk());
    EXPECT_DOUBLE_EQ(state.internalForce()[0U], -20.0);
    EXPECT_DOUBLE_EQ(state.internalForce()[6U], 20.0);
    EXPECT_DOUBLE_EQ(state.residual()[0U], -20.0);
    EXPECT_DOUBLE_EQ(state.residual()[6U], 0.0);
    EXPECT_DOUBLE_EQ(state.reaction()[0U], -20.0);
    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) {
    const auto fixture = makeFixture();

    auto failed = makeAxialEquilibriumState(fixture);
    failed.internalForce()[0U] = 91.0;
    failed.residual()[0U] = 92.0;
    failed.reaction()[0U] = 93.0;
    failed.reaction()[6U] = 94.0;
    failed.endpointResults().push_back({});
    failed.externalForce()[6U] += 1.0e-7;
    expectStatusCode(
        fesa::ResultRecovery::recover(
            *fixture.model, *fixture.dofs, *fixture.stiffness, failed),
        "free-residual-tolerance-failure");
    EXPECT_DOUBLE_EQ(failed.internalForce()[0U], 91.0);
    EXPECT_DOUBLE_EQ(failed.residual()[0U], 92.0);
    EXPECT_DOUBLE_EQ(failed.reaction()[0U], 93.0);
    EXPECT_DOUBLE_EQ(failed.reaction()[6U], 94.0);
    EXPECT_EQ(failed.endpointResults().size(), 1U);

    auto thresholdPass = makeAxialEquilibriumState(fixture);
    thresholdPass.externalForce()[6U] += 1.0e-10 * 20.0 * 0.5;
    const fesa::Status thresholdStatus = fesa::ResultRecovery::recover(
        *fixture.model,
        *fixture.dofs,
        *fixture.stiffness,
        thresholdPass);
    ASSERT_TRUE(thresholdStatus.IsOk());
    EXPECT_NE(thresholdPass.residual()[6U], 0.0);
    EXPECT_DOUBLE_EQ(
        thresholdPass.reaction()[6U], thresholdPass.residual()[6U]);

    const auto zeroFixture = makeFixture(false, {}, {}, false, false, false);
    auto zeroEquilibrium = fesa::AnalysisState::create(
        *zeroFixture.dofs, {"Step-1", 0U});
    EXPECT_TRUE(fesa::ResultRecovery::recover(
                    *zeroFixture.model,
                    *zeroFixture.dofs,
                    *zeroFixture.stiffness,
                    zeroEquilibrium)
                    .IsOk());

    auto wrongPrescription = makeAxialEquilibriumState(fixture);
    wrongPrescription.displacement()[0U] = 0.0;
    expectStatusCode(
        fesa::ResultRecovery::recover(
            *fixture.model,
            *fixture.dofs,
            *fixture.stiffness,
            wrongPrescription),
        "invalid-recovery-state");

    auto nonfinite = makeAxialEquilibriumState(fixture);
    nonfinite.displacement()[6U] = std::numeric_limits<double>::quiet_NaN();
    expectStatusCode(
        fesa::ResultRecovery::recover(
            *fixture.model, *fixture.dofs, *fixture.stiffness, nonfinite),
        "nonfinite-recovery-value");

    const auto wrongFixture = makeFixture(true);
    auto wrongState = fesa::AnalysisState::create(
        *wrongFixture.dofs, {"Step-1", 0U});
    expectStatusCode(
        fesa::ResultRecovery::recover(
            *fixture.model, *fixture.dofs, *fixture.stiffness, wrongState),
        "invalid-recovery-dimensions");
}

TEST(ResultRecovery, KeepsEndActionSectionAndGaussResultsDistinct) {
    const auto fixture = makeFixture();
    auto state = makePatchState(fixture, 0.02, 0.03, -0.04, 0.05);

    ASSERT_TRUE(fesa::ResultRecovery::recover(
                    *fixture.model, *fixture.dofs, *fixture.stiffness, state)
                    .IsOk());
    ASSERT_EQ(state.endpointResults().size(), 2U);
    ASSERT_EQ(state.gaussResults().size(), 2U);
    EXPECT_EQ(state.endpointResults()[0U].endpoint, 0);
    EXPECT_EQ(state.endpointResults()[1U].endpoint, 1);
    EXPECT_EQ(state.gaussResults()[0U].gaussPoint, 1);
    EXPECT_EQ(state.gaussResults()[1U].gaussPoint, 2);
    EXPECT_DOUBLE_EQ(
        state.endpointResults()[0U].endAction[0U],
        -state.endpointResults()[0U].sectionResultant[0U]);
    EXPECT_DOUBLE_EQ(
        state.endpointResults()[1U].endAction[0U],
        state.endpointResults()[1U].sectionResultant[0U]);
    EXPECT_DOUBLE_EQ(
        state.gaussResults()[0U].generalizedResultant[0U],
        state.endpointResults()[0U].sectionResultant[0U]);
}

TEST(ResultRecovery, MatchesAxialTorsionAndTwoPlaneEndSigns) {
    const auto fixture = makeFixture();
    const double epsilon = 0.02;
    const double twist = -0.03;
    const double kappaY = 0.04;
    const double kappaZ = -0.05;
    auto state = makePatchState(fixture, epsilon, twist, kappaY, kappaZ);

    ASSERT_TRUE(fesa::ResultRecovery::recover(
                    *fixture.model, *fixture.dofs, *fixture.stiffness, state)
                    .IsOk());
    const double shearModulus =
        kYoungsModulus / (2.0 * (1.0 + kPoissonRatio));
    const std::array<double, 4> expected = {
        kYoungsModulus * 2.0 * epsilon,
        shearModulus * 5.0 * twist,
        kYoungsModulus * 3.0 * kappaY,
        kYoungsModulus * 4.0 * kappaZ};
    const std::array<std::size_t, 4> endComponents = {0U, 3U, 4U, 5U};
    for (std::size_t component = 0U; component < expected.size(); ++component) {
        expectScaledNear(
            state.endpointResults()[0U].sectionResultant[component],
            expected[component]);
        expectScaledNear(
            state.endpointResults()[1U].sectionResultant[component],
            expected[component]);
        expectScaledNear(
            state.endpointResults()[0U].endAction[endComponents[component]],
            -expected[component]);
        expectScaledNear(
            state.endpointResults()[1U].endAction[endComponents[component]],
            expected[component]);
    }
}

TEST(ResultRecovery, OrdersStressPointsAndDefaultCentroid) {
    const std::vector<std::array<double, 2>> sectionPoints = {
        {0.25, -0.5}, {-0.4, 0.3}};
    const auto fixture = makeFixture(false, sectionPoints);
    auto state = makePatchState(fixture, 0.01, 0.0, 0.02, -0.03);
    ASSERT_TRUE(fesa::ResultRecovery::recover(
                    *fixture.model, *fixture.dofs, *fixture.stiffness, state)
                    .IsOk());
    ASSERT_EQ(state.stressResults().size(), 4U);
    for (std::size_t gauss = 0U; gauss < 2U; ++gauss) {
        for (std::size_t point = 0U; point < sectionPoints.size(); ++point) {
            const auto& row = state.stressResults()[gauss * 2U + point];
            EXPECT_EQ(row.element, 0U);
            EXPECT_EQ(row.gaussPoint, static_cast<int>(gauss + 1U));
            EXPECT_EQ(row.sectionPoint, point + 1U);
            EXPECT_DOUBLE_EQ(row.x1, sectionPoints[point][0U]);
            EXPECT_DOUBLE_EQ(row.x2, sectionPoints[point][1U]);
            EXPECT_EQ(row.source, "input");
            expectScaledNear(
                row.s11,
                kYoungsModulus *
                    (0.01 + row.x2 * 0.02 - row.x1 * -0.03));
        }
    }

    const auto defaultFixture = makeFixture();
    auto defaultState = makePatchState(defaultFixture, 0.01, 0.0, 0.0, 0.0);
    ASSERT_TRUE(fesa::ResultRecovery::recover(
                    *defaultFixture.model,
                    *defaultFixture.dofs,
                    *defaultFixture.stiffness,
                    defaultState)
                    .IsOk());
    ASSERT_EQ(defaultState.stressResults().size(), 2U);
    for (const auto& row : defaultState.stressResults()) {
        EXPECT_EQ(row.sectionPoint, 0U);
        EXPECT_DOUBLE_EQ(row.x1, 0.0);
        EXPECT_DOUBLE_EQ(row.x2, 0.0);
        EXPECT_EQ(row.source, "fesa-default");
    }
}

TEST(ResultRecovery, RequiresInteriorEndpointConsistencyWithoutAveraging) {
    const auto fixture = makeFixture(true);
    const std::array<double, 4> tolerances = {1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6};
    auto rows = makeStationRows(fixture);
    rows[2U].sectionResultant[0U] += 0.5e-6;

    auto normalized =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *fixture.model, rows, tolerances);
    ASSERT_TRUE(normalized.HasValue());
    ASSERT_EQ(normalized.Value().size(), 3U);
    EXPECT_EQ(normalized.Value()[1U].representativeElement, 0U);
    EXPECT_DOUBLE_EQ(normalized.Value()[1U].sectionResultant[0U], 5.0);

    rows[2U].sectionResultant[0U] = 5.0 + 2.0e-6;
    auto mismatch =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *fixture.model, rows, tolerances);
    ASSERT_FALSE(mismatch.HasValue());
    expectStatusCode(mismatch.GetStatus(), "node-station-tolerance-failure");

    rows = makeStationRows(fixture);
    rows[2U].sectionResultant[1U] =
        std::numeric_limits<double>::infinity();
    auto nonfinite =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *fixture.model, rows, tolerances);
    ASSERT_FALSE(nonfinite.HasValue());
    expectStatusCode(nonfinite.GetStatus(), "nonfinite-node-station-value");

    auto invalidTolerance =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *fixture.model,
            makeStationRows(fixture),
            {1.0e-6, -1.0, 1.0e-6, 1.0e-6});
    ASSERT_FALSE(invalidTolerance.HasValue());
    expectStatusCode(
        invalidTolerance.GetStatus(), "invalid-node-station-tolerance");

    const std::filesystem::path source{"models/result-recovery.inp"};
    const auto loadedFixture = makeFixture(
        true, {}, {{"2", 2, 1.0, {source, 60U}}});
    auto loaded =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *loadedFixture.model,
            makeStationRows(loadedFixture),
            tolerances);
    ASSERT_FALSE(loaded.HasValue());
    expectStatusCode(loaded.GetStatus(), "ineligible-node-station");

    const auto reversedFixture = makeFixture(true, {}, {}, true);
    auto reversed =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *reversedFixture.model,
            makeStationRows(reversedFixture),
            tolerances);
    ASSERT_FALSE(reversed.HasValue());
    expectStatusCode(reversed.GetStatus(), "ineligible-node-station");

    const auto jumpFixture = makeFixture(true, {}, {}, false, true);
    auto jumped =
        fesa::ResultRecovery::normalizeSectionResultantsToNodeStations(
            *jumpFixture.model,
            makeStationRows(jumpFixture),
            tolerances);
    ASSERT_FALSE(jumped.HasValue());
    expectStatusCode(jumped.GetStatus(), "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);
}