#include "fesa/results/result_recovery.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>

#include "fesa/analysis/analysis_model.h"
#include "fesa/analysis/analysis_state.h"
#include "fesa/assembly/load_assembler.h"
#include "fesa/assembly/parallel_for.h"
#include "fesa/assembly/sparse_assembler.h"
#include "fesa/fem/dof_manager.h"
#include "fesa/math/vector3.h"
#include "fesa/model/domain.h"

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 two_elements = false,
    std::vector<std::array<double, 2>> section_points = {},
    std::vector<fesa::NodalLoad> loads = {}, const bool reverse_second = false,
    const bool section_jump = false, const bool nonzero_prescription = 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 (two_elements) {
    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(section_points),
                          {source, 30U}}};
  if (section_jump) {
    auto second_section = definition.sections.front();
    second_section.name = "Section-2";
    second_section.area = 2.5;
    second_section.location.line = 31U;
    definition.sections.push_back(std::move(second_section));
  }
  definition.elements = {
      {{"Beam-1", 10, "10"}, {0U, 1U}, 0U, 0U, {source, 40U}}};
  if (two_elements) {
    definition.elements.push_back(
        {{"Beam-1", 20, "20"},
         reverse_second ? std::array<fesa::EntityIndex, 2>{2U, 1U}
                        : std::array<fesa::EntityIndex, 2>{1U, 2U},
         0U,
         section_jump ? 1U : 0U,
         {source, 41U}});
  }
  definition.steps = {
      {"Step-1",
       {{"1", 1, 1, nonzero_prescription ? 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 two_elements = false,
    std::vector<std::array<double, 2>> section_points = {},
    std::vector<fesa::NodalLoad> loads = {}, const bool reverse_second = false,
    const bool section_jump = false, const bool nonzero_prescription = true) {
  auto domain_result = fesa::Domain::Create(
      MakeDefinition(two_elements, std::move(section_points), std::move(loads),
                     reverse_second, section_jump, nonzero_prescription));
  if (!domain_result.HasValue()) {
    throw std::runtime_error{"Recovery fixture Domain construction failed."};
  }
  auto domain =
      std::make_unique<fesa::Domain>(std::move(domain_result.Value()));

  auto model_result = fesa::AnalysisModel::Create(*domain);
  if (!model_result.HasValue()) {
    throw std::runtime_error{
        "Recovery fixture AnalysisModel construction failed."};
  }
  auto model =
      std::make_unique<fesa::AnalysisModel>(std::move(model_result.Value()));

  auto dofs_result = fesa::DofManager::Create(*model);
  if (!dofs_result.HasValue()) {
    throw std::runtime_error{
        "Recovery fixture DofManager construction failed."};
  }
  auto dofs =
      std::make_unique<fesa::DofManager>(std::move(dofs_result.Value()));

  fesa::SerialParallelFor serial;
  auto stiffness_result =
      fesa::SparseAssembler::AssembleStiffness(*model, *dofs, serial);
  if (!stiffness_result.HasValue()) {
    throw std::runtime_error{"Recovery fixture stiffness assembly failed."};
  }
  auto stiffness =
      std::make_unique<fesa::SparseMatrix>(std::move(stiffness_result.Value()));
  return {std::move(domain), std::move(model), std::move(dofs),
          std::move(stiffness)};
}

fesa::ModelDefinition MakeShellDefinition(
    const bool two_elements = false, const bool constrain_all = 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 (!two_elements) {
    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, two_elements ? 4U : 3U, two_elements ? 3U : 2U},
       0U,
       0U,
       {source, 40U}}};
  if (two_elements) {
    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 (constrain_all) {
    std::vector<fesa::EntityIndex> all_nodes;
    all_nodes.reserve(definition.nodes.size());
    for (std::size_t node = 0U; node < definition.nodes.size(); ++node) {
      all_nodes.push_back(static_cast<fesa::EntityIndex>(node));
    }
    definition.node_sets.push_back(
        {"All", {}, std::move(all_nodes), {source, 50U}});
  }
  definition.steps = {
      {"Step-1",
       constrain_all ? 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 domain_result = fesa::Domain::Create(std::move(definition));
  if (!domain_result.HasValue()) {
    throw std::runtime_error{
        "Shell recovery fixture Domain construction failed."};
  }
  auto domain =
      std::make_unique<fesa::Domain>(std::move(domain_result.Value()));

  auto model_result = fesa::AnalysisModel::Create(*domain);
  if (!model_result.HasValue()) {
    throw std::runtime_error{
        "Shell recovery fixture AnalysisModel construction failed."};
  }
  auto model =
      std::make_unique<fesa::AnalysisModel>(std::move(model_result.Value()));

  auto dofs_result = fesa::DofManager::Create(*model);
  if (!dofs_result.HasValue()) {
    throw std::runtime_error{
        "Shell recovery fixture DofManager construction failed."};
  }
  auto dofs =
      std::make_unique<fesa::DofManager>(std::move(dofs_result.Value()));

  fesa::SerialParallelFor serial;
  auto stiffness_result =
      fesa::SparseAssembler::AssembleStiffness(*model, *dofs, serial);
  if (!stiffness_result.HasValue()) {
    throw std::runtime_error{
        "Shell recovery fixture stiffness assembly failed."};
  }
  auto stiffness =
      std::make_unique<fesa::SparseMatrix>(std::move(stiffness_result.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 full_dof : fixture.dofs->FreeDofs()) {
    state.ExternalForce()[full_dof] = internal[full_dof];
  }
  return state;
}

fesa::AnalysisState MakePatchState(const RecoveryFixture& fixture,
                                   const double epsilon, const double twist,
                                   const double kappa_y, const double kappa_z) {
  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 * kappa_z * kLength * kLength;
  state.Displacement()[8U] = -0.5 * kappa_y * kLength * kLength;
  state.Displacement()[9U] = twist * kLength;
  state.Displacement()[10U] = kappa_y * kLength;
  state.Displacement()[11U] = kappa_z * 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 relative_tolerance = 1.0e-12) {
  ASSERT_TRUE(std::isfinite(actual));
  ASSERT_TRUE(std::isfinite(expected));
  EXPECT_LE(std::abs(actual - expected),
            relative_tolerance * (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 stale_shell_evidence{};
  stale_shell_evidence.physical_strain_energy = 123.0;
  stale_shell_evidence.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
  stale_shell_evidence.verification_metrics = {1.0e-11, 2.0e-11, 3.0e-11};
  ASSERT_TRUE(state.CommitShellResults({}, stale_shell_evidence).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 full_dof : fixture.dofs->FreeDofs()) {
    EXPECT_DOUBLE_EQ(state.Reaction()[full_dof], 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 threshold_pass = MakeAxialEquilibriumState(fixture);
  threshold_pass.ExternalForce()[6U] += 1.0e-10 * 20.0 * 0.5;
  const fesa::Status threshold_status = fesa::ResultRecovery::Recover(
      *fixture.model, *fixture.dofs, *fixture.stiffness, threshold_pass);
  ASSERT_TRUE(threshold_status.IsOk());
  EXPECT_NE(threshold_pass.Residual()[6U], 0.0);
  EXPECT_DOUBLE_EQ(threshold_pass.Reaction()[6U],
                   threshold_pass.Residual()[6U]);

  const auto zero_fixture = MakeFixture(false, {}, {}, false, false, false);
  auto zero_equilibrium =
      fesa::AnalysisState::Create(*zero_fixture.dofs, {"Step-1", 0U});
  EXPECT_TRUE(
      fesa::ResultRecovery::Recover(*zero_fixture.model, *zero_fixture.dofs,
                                    *zero_fixture.stiffness, zero_equilibrium)
          .IsOk());

  auto wrong_prescription = MakeAxialEquilibriumState(fixture);
  wrong_prescription.Displacement()[0U] = 0.0;
  ExpectStatusCode(
      fesa::ResultRecovery::Recover(*fixture.model, *fixture.dofs,
                                    *fixture.stiffness, wrong_prescription),
      "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 wrong_fixture = MakeFixture(true);
  auto wrong_state =
      fesa::AnalysisState::Create(*wrong_fixture.dofs, {"Step-1", 0U});
  ExpectStatusCode(
      fesa::ResultRecovery::Recover(*fixture.model, *fixture.dofs,
                                    *fixture.stiffness, wrong_state),
      "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].gauss_point, 1);
  EXPECT_EQ(state.GaussResults()[1U].gauss_point, 2);
  EXPECT_DOUBLE_EQ(state.EndpointResults()[0U].end_action[0U],
                   -state.EndpointResults()[0U].section_resultant[0U]);
  EXPECT_DOUBLE_EQ(state.EndpointResults()[1U].end_action[0U],
                   state.EndpointResults()[1U].section_resultant[0U]);
  EXPECT_DOUBLE_EQ(state.GaussResults()[0U].generalized_resultant[0U],
                   state.EndpointResults()[0U].section_resultant[0U]);
}

TEST(ResultRecovery, MatchesAxialTorsionAndTwoPlaneEndSigns) {
  const auto fixture = MakeFixture();
  const double epsilon = 0.02;
  const double twist = -0.03;
  const double kappa_y = 0.04;
  const double kappa_z = -0.05;
  auto state = MakePatchState(fixture, epsilon, twist, kappa_y, kappa_z);

  ASSERT_TRUE(fesa::ResultRecovery::Recover(*fixture.model, *fixture.dofs,
                                            *fixture.stiffness, state)
                  .IsOk());
  const double shear_modulus = kYoungsModulus / (2.0 * (1.0 + kPoissonRatio));
  const std::array<double, 4> expected = {
      kYoungsModulus * 2.0 * epsilon, shear_modulus * 5.0 * twist,
      kYoungsModulus * 3.0 * kappa_y, kYoungsModulus * 4.0 * kappa_z};
  const std::array<std::size_t, 4> end_components = {0U, 3U, 4U, 5U};
  for (std::size_t component = 0U; component < expected.size(); ++component) {
    ExpectScaledNear(state.EndpointResults()[0U].section_resultant[component],
                     expected[component]);
    ExpectScaledNear(state.EndpointResults()[1U].section_resultant[component],
                     expected[component]);
    ExpectScaledNear(
        state.EndpointResults()[0U].end_action[end_components[component]],
        -expected[component]);
    ExpectScaledNear(
        state.EndpointResults()[1U].end_action[end_components[component]],
        expected[component]);
  }
}

TEST(ResultRecovery, OrdersStressPointsAndDefaultCentroid) {
  const std::vector<std::array<double, 2>> section_points = {{0.25, -0.5},
                                                             {-0.4, 0.3}};
  const auto fixture = MakeFixture(false, section_points);
  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 < section_points.size(); ++point) {
      const auto& row = state.StressResults()[gauss * 2U + point];
      EXPECT_EQ(row.element, 0U);
      EXPECT_EQ(row.gauss_point, static_cast<int>(gauss + 1U));
      EXPECT_EQ(row.section_point, point + 1U);
      EXPECT_DOUBLE_EQ(row.x1, section_points[point][0U]);
      EXPECT_DOUBLE_EQ(row.x2, section_points[point][1U]);
      EXPECT_EQ(row.source, "input");
      ExpectScaledNear(
          row.s11, kYoungsModulus * (0.01 + row.x2 * 0.02 - row.x1 * -0.03));
    }
  }

  const auto default_fixture = MakeFixture();
  auto default_state = MakePatchState(default_fixture, 0.01, 0.0, 0.0, 0.0);
  ASSERT_TRUE(fesa::ResultRecovery::Recover(
                  *default_fixture.model, *default_fixture.dofs,
                  *default_fixture.stiffness, default_state)
                  .IsOk());
  ASSERT_EQ(default_state.StressResults().size(), 2U);
  for (const auto& row : default_state.StressResults()) {
    EXPECT_EQ(row.section_point, 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].section_resultant[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].representative_element, 0U);
  EXPECT_DOUBLE_EQ(normalized.Value()[1U].section_resultant[0U], 5.0);

  rows[2U].section_resultant[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].section_resultant[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 invalid_tolerance =
      fesa::ResultRecovery::NormalizeSectionResultantsToNodeStations(
          *fixture.model, MakeStationRows(fixture),
          {1.0e-6, -1.0, 1.0e-6, 1.0e-6});
  ASSERT_FALSE(invalid_tolerance.HasValue());
  ExpectStatusCode(invalid_tolerance.GetStatus(),
                   "invalid-node-station-tolerance");

  const std::filesystem::path source{"models/result-recovery.inp"};
  const auto loaded_fixture =
      MakeFixture(true, {}, {{"2", 2, 1.0, {source, 60U}}});
  auto loaded = fesa::ResultRecovery::NormalizeSectionResultantsToNodeStations(
      *loaded_fixture.model, MakeStationRows(loaded_fixture), tolerances);
  ASSERT_FALSE(loaded.HasValue());
  ExpectStatusCode(loaded.GetStatus(), "ineligible-node-station");

  const auto reversed_fixture = MakeFixture(true, {}, {}, true);
  auto reversed =
      fesa::ResultRecovery::NormalizeSectionResultantsToNodeStations(
          *reversed_fixture.model, MakeStationRows(reversed_fixture),
          tolerances);
  ASSERT_FALSE(reversed.HasValue());
  ExpectStatusCode(reversed.GetStatus(), "ineligible-node-station");

  const auto jump_fixture = MakeFixture(true, {}, {}, false, true);
  auto jumped = fesa::ResultRecovery::NormalizeSectionResultantsToNodeStations(
      *jump_fixture.model, MakeStationRows(jump_fixture), 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::kGp1, fesa::ShellMidsurfaceLocation::kGp2,
      fesa::ShellMidsurfaceLocation::kGp3, fesa::ShellMidsurfaceLocation::kGp4};
  const std::array<std::array<double, 2>, 4> coordinates{
      std::array<double, 2>{-gauss, -gauss},
      std::array<double, 2>{gauss, -gauss}, std::array<double, 2>{gauss, gauss},
      std::array<double, 2>{-gauss, gauss}};
  constexpr std::array<double, 8> expected_strain{0.1, -0.05, 0.2, 0.0,
                                                  0.0, 0.0,   0.0, 0.0};
  constexpr std::array<double, 8> expected_resultant{22.4, -6.4, 19.2, 0.0,
                                                     0.0,  0.0,  0.0,  0.0};
  const std::array<std::array<double, 3>, 3> expected_frame{
      {{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.natural_coordinates, coordinates[point]);
      EXPECT_EQ(row.local_frame, expected_frame);
      for (std::size_t component = 0U; component < expected_strain.size();
           ++component) {
        EXPECT_NEAR(row.generalized_strain[component],
                    expected_strain[component], 1.0e-12);
        EXPECT_NEAR(row.section_resultant[component],
                    expected_resultant[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::kBottom, fesa::ShellSectionPosition::kMiddle,
      fesa::ShellSectionPosition::kTop};
  constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};
  constexpr std::array<std::array<double, 3>, 3> expected_stress{
      {{-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 < expected_stress[position].size(); ++component) {
        EXPECT_NEAR(row.stress[position].components[component],
                    expected_stress[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 stabilized_energy =
      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(stabilized_energy, 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 full_load =
      fesa::LoadAssembler::AssembleFullNodalLoad(*fixture.model, *fixture.dofs);
  ASSERT_TRUE(full_load.HasValue());
  state.ExternalForce() = std::move(full_load.Value());

  ASSERT_TRUE(fesa::ResultRecovery::Recover(*fixture.model, *fixture.dofs,
                                            *fixture.stiffness, state)
                  .IsOk());
  ASSERT_EQ(state.ShellResults().size(), 4U);
  for (std::size_t full_dof = 0U; full_dof < fixture.dofs->FullDofCount();
       ++full_dof) {
    EXPECT_DOUBLE_EQ(state.InternalForce()[full_dof], 0.0);
    EXPECT_DOUBLE_EQ(state.Residual()[full_dof],
                     -state.ExternalForce()[full_dof]);
    EXPECT_DOUBLE_EQ(state.Reaction()[full_dof], state.Residual()[full_dof]);
  }
  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 free_fixture = MakeShellFixture(MakeShellDefinition());
  auto perturbed = MakeShellPhysicalState(free_fixture);
  perturbed.ExternalForce()[0U] += 1.0e-9;
  ASSERT_TRUE(fesa::ResultRecovery::Recover(*free_fixture.model,
                                            *free_fixture.dofs,
                                            *free_fixture.stiffness, perturbed)
                  .IsOk());
  for (const double metric : perturbed.VerificationMetrics()) {
    EXPECT_GT(metric, 0.0);
    EXPECT_LE(metric, 1.0e-10);
  }
}

// MITC4-REC-004, C-DUP-002
TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) {
  auto centered_definition = MakeShellDefinition();
  auto translated_definition = centered_definition;
  constexpr fesa::Vector3 translation{7.0, 11.0, 0.0};
  for (auto& node : translated_definition.nodes) {
    for (std::size_t component = 0U;
         component < translation.Components().size(); ++component) {
      node.coordinates[component] += translation[component];
    }
  }
  const auto centered_fixture =
      MakeShellFixture(std::move(centered_definition));
  const auto translated_fixture =
      MakeShellFixture(std::move(translated_definition));
  auto centered = MakeShellPhysicalState(centered_fixture);
  auto translated = MakeShellPhysicalState(translated_fixture);
  centered.ExternalForce()[0U] += 1.0e-9;
  translated.ExternalForce()[0U] += 1.0e-9;

  ASSERT_TRUE(fesa::ResultRecovery::Recover(
                  *centered_fixture.model, *centered_fixture.dofs,
                  *centered_fixture.stiffness, centered)
                  .IsOk());
  ASSERT_TRUE(fesa::ResultRecovery::Recover(
                  *translated_fixture.model, *translated_fixture.dofs,
                  *translated_fixture.stiffness, translated)
                  .IsOk());
  std::array<double, 3> centered_force{};
  for (std::size_t component = 0U; component < 3U; ++component) {
    centered_force[component] = centered.Equilibrium()[component];
    EXPECT_NEAR(translated.Equilibrium()[component], centered_force[component],
                1.0e-12);
  }
  const fesa::Vector3 translated_moment_delta =
      translation.Cross(fesa::Vector3{centered_force});
  for (std::size_t component = 0U; component < 3U; ++component) {
    EXPECT_NEAR(translated.Equilibrium()[3U + component] -
                    centered.Equilibrium()[3U + component],
                translated_moment_delta[component], 5.0e-12);
  }
  EXPECT_GT(std::abs(translated_moment_delta[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 subunit_scale = 1.0e-6;
  constexpr double large_scale = 1.0e6;
  subunit.Displacement().Scale(subunit_scale);
  subunit.ExternalForce().Scale(subunit_scale);
  subunit.ExternalForce()[0U] += 1.0e-9 * subunit_scale;
  large.Displacement().Scale(large_scale);
  large.ExternalForce().Scale(large_scale);
  large.ExternalForce()[0U] += 1.0e-9 * large_scale;

  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 constrained_definition = MakeShellDefinition(false, true);
  constrained_definition.steps[0U].boundaries[0U].value = 1.0;
  const auto constrained_fixture =
      MakeShellFixture(std::move(constrained_definition));
  auto rejected =
      fesa::AnalysisState::Create(*constrained_fixture.dofs, {"Step-1", 0U});
  for (std::size_t full_dof = 0U;
       full_dof < constrained_fixture.dofs->FullDofCount(); ++full_dof) {
    rejected.Displacement()[full_dof] = 1.0;
  }
  const std::vector<fesa::CooContribution> unbalanced_entry{
      {0U, 0U, 1.0, 0U, 0U}};
  auto unbalanced_stiffness = fesa::SparseMatrix::FromCoo(
      constrained_fixture.dofs->FullDofCount(),
      constrained_fixture.dofs->FullDofCount(), unbalanced_entry,
      constrained_fixture.dofs->GetSparsePattern());
  ASSERT_TRUE(unbalanced_stiffness.HasValue());
  ExpectStatusCode(fesa::ResultRecovery::Recover(
                       *constrained_fixture.model, *constrained_fixture.dofs,
                       unbalanced_stiffness.Value(), rejected),
                   "global-equilibrium-tolerance-failure");
}

// MITC4-REC-005
TEST(ResultRecovery, InvalidLaterShellLeavesEntirePriorStateUnchanged) {
  const auto valid_fixture = MakeShellFixture(MakeShellDefinition(true));
  auto state = MakeShellPhysicalState(valid_fixture);
  ASSERT_TRUE(fesa::ResultRecovery::Recover(*valid_fixture.model,
                                            *valid_fixture.dofs,
                                            *valid_fixture.stiffness, state)
                  .IsOk());
  ASSERT_EQ(state.ShellResults().size(), 8U);
  const auto prior_first_row = state.ShellResults().front();
  const double prior_energy = state.PhysicalStrainEnergy();
  const auto prior_equilibrium = state.Equilibrium();
  const auto prior_metrics = 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{valid_fixture.dofs->FullDofCount()};
  auto zero_stiffness = fesa::SparseMatrix::FromCoo(
      valid_fixture.dofs->FullDofCount(), valid_fixture.dofs->FullDofCount(),
      {}, valid_fixture.dofs->GetSparsePattern());
  ASSERT_TRUE(zero_stiffness.HasValue());

  const auto status = fesa::ResultRecovery::Recover(
      *valid_fixture.model, *valid_fixture.dofs, zero_stiffness.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, prior_first_row.element);
  EXPECT_EQ(state.ShellResults().front().location, prior_first_row.location);
  EXPECT_EQ(state.ShellResults().front().generalized_strain,
            prior_first_row.generalized_strain);
  EXPECT_EQ(state.ShellResults().front().section_resultant,
            prior_first_row.section_resultant);
  EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), prior_energy);
  EXPECT_EQ(state.Equilibrium(), prior_equilibrium);
  EXPECT_EQ(state.VerificationMetrics(), prior_metrics);
}