#include <gtest/gtest.h>

#include <array>
#include <cmath>
#include <cstdint>
#include <string>

#include "fesa/elements/mitc4_shell.h"
#include "fesa/model/model_types.h"

namespace {

using Coordinates = std::array<double, 3>;

fesa::Node Node(std::int64_t label, Coordinates coordinates) {
  return {{"I", label, std::to_string(label)}, coordinates, {"mass.inp", 1U}};
}

std::array<const fesa::Node*, 4> Pointers(
    const std::array<fesa::Node, 4>& nodes) {
  return {&nodes[0], &nodes[1], &nodes[2], &nodes[3]};
}

fesa::Mitc4Shell WarpedShell() {
  const std::array<fesa::Node, 4> nodes{
      Node(1, {-1.0, -1.0, 0.1}), Node(2, {1.0, -1.0, -0.1}),
      Node(3, {1.0, 1.0, 0.2}), Node(4, {-1.0, 1.0, -0.1})};
  const std::array<Coordinates, 4> directors{
      {{{0.0, 0.0, 1.0}},
       {{0.6, 0.0, 0.8}},
       {{2.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0}},
       {{0.0, 0.6, 0.8}}}};
  const fesa::ShellSection section{"S", 0.5, 0U, {"mass.inp", 2U}};
  const fesa::LinearElasticMaterial material{"M", 12.0, 0.0, {"mass.inp", 3U}};
  const auto result =
      fesa::Mitc4Shell::Create(Pointers(nodes), directors, section, material);
  EXPECT_TRUE(result.HasValue());
  return result.Value();
}

double Quadratic(const fesa::Matrix& matrix, const fesa::Vector& vector) {
  return vector.Dot(matrix.Multiply(vector));
}

// E2-C02: values derive independently in numerical-review-evidence.md.
TEST(Mitc4ShellMassContract, PreservesWarpedMomentsRowsAndFirstSpatialMoment) {
  const auto mass = WarpedShell().Mass(3.0);
  ASSERT_TRUE(mass.HasValue());
  const auto& value = mass.Value();
  EXPECT_NEAR(value.consistent_moments[0U](0U, 0U), 17609.0 / 28800.0, 1e-12);
  EXPECT_NEAR(value.consistent_moments[1U](0U, 0U), 11003.0 / 1440000.0, 1e-12);
  EXPECT_NEAR(value.consistent_moments[2U](0U, 0U), 88177.0 / 6912000.0, 1e-12);
  double total_m = 0.0, total_h = 0.0, total_j = 0.0;
  for (const auto& row : value.nodal_moments) {
    total_m += row[0U];
    total_h += row[1U];
    total_j += row[2U];
  }
  EXPECT_NEAR(total_m, 7743.0 / 1600.0, 1e-12);
  EXPECT_NEAR(total_h, 739.0 / 12000.0, 1e-12);
  EXPECT_NEAR(total_j, 517.0 / 5120.0, 1e-12);
}

// E2-C03: scalar row-sum lumping retains uniform translation energy and is
// never below the consistent kinetic energy.
TEST(Mitc4ShellMassContract, SatisfiesUniformAndJensenEnergyContracts) {
  const auto mass = WarpedShell().Mass(3.0);
  ASSERT_TRUE(mass.HasValue());
  fesa::Vector translation{24U};
  fesa::Vector field{24U};
  for (std::size_t node = 0U; node < 4U; ++node) {
    translation[6U * node] = 1.0;
    field[6U * node] = static_cast<double>(node + 1U);
    field[6U * node + 1U] = -0.5 * static_cast<double>(node);
    field[6U * node + 3U] = 0.25;
  }
  EXPECT_NEAR(Quadratic(mass.Value().consistent_global24, translation),
              Quadratic(mass.Value().lumped_global24, translation), 1e-12);
  EXPECT_GE(Quadratic(mass.Value().lumped_global24, field) + 1e-12,
            Quadratic(mass.Value().consistent_global24, field));
}

// E2-C04: assemble the physical velocity work over successively finer meshes.
TEST(Mitc4ShellMassContract, ReproducesAffineSquareVelocityWork) {
  const std::array<Coordinates, 4> directors{{{{0.0, 0.0, 1.0}},
                                              {{0.0, 0.0, 1.0}},
                                              {{0.0, 0.0, 1.0}},
                                              {{0.0, 0.0, 1.0}}}};
  const fesa::ShellSection section{"S", 1.0, 0U, {"mass.inp", 2U}};
  const fesa::LinearElasticMaterial material{"M", 12.0, 0.0, {"mass.inp", 3U}};
  for (const int mesh : {1, 2, 4}) {
    double consistent_work = 0.0, lumped_work = 0.0;
    const double edge = 2.0 / mesh;
    for (int row = 0; row < mesh; ++row) {
      for (int column = 0; column < mesh; ++column) {
        const double x = -1.0 + edge * column, y = -1.0 + edge * row;
        const std::array<fesa::Node, 4> nodes{
            Node(1, {x, y, 0.0}), Node(2, {x + edge, y, 0.0}),
            Node(3, {x + edge, y + edge, 0.0}), Node(4, {x, y + edge, 0.0})};
        const auto shell = fesa::Mitc4Shell::Create(Pointers(nodes), directors,
                                                    section, material);
        ASSERT_TRUE(shell.HasValue());
        const auto mass = shell.Value().Mass(1.0);
        ASSERT_TRUE(mass.HasValue());
        fesa::Vector velocity{24U};
        for (std::size_t node = 0U; node < 4U; ++node)
          velocity[6U * node] = nodes[node].coordinates[0U];
        consistent_work +=
            Quadratic(mass.Value().consistent_global24, velocity);
        lumped_work += Quadratic(mass.Value().lumped_global24, velocity);
      }
    }
    EXPECT_NEAR(consistent_work, 4.0 / 3.0, 1e-12) << mesh;
    EXPECT_NEAR(lumped_work - consistent_work, 8.0 / (3.0 * mesh * mesh), 1e-12)
        << mesh;
  }
}

}  // namespace