#include "fesa/elements/nonlinear_mitc4_shell.h"

#include <gtest/gtest.h>

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

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

namespace {

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

fesa::Node Node(const std::int64_t label, const Vector3 coordinates) {
  return {{"Shell-Instance", label, std::to_string(label)}, coordinates,
          {"nonlinear-mitc4.inp", static_cast<std::size_t>(label)}};
}

std::array<fesa::Node, 4> Nodes() {
  return {Node(1, {-1.0, -1.0, 0.0}), Node(2, {1.0, -1.0, 0.0}),
          Node(3, {1.0, 1.0, 0.0}), Node(4, {-1.0, 1.0, 0.0})};
}

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

std::array<fesa::ShellNodeInitialFrame, 4> Frames() {
  return {{{0U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
           {1U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
           {2U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
           {3U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}}}};
}

fesa::ShellSection Section() {
  return {"Section", 0.2, 0U, {"nonlinear-mitc4.inp", 20U}};
}

fesa::LinearElasticMaterial Material() {
  return {"Material", 120.0, 0.25, {"nonlinear-mitc4.inp", 21U}};
}

fesa::NonlinearMitc4TrialDofs ZeroTrial() { return {}; }

fesa::Vector3 Multiply(const fesa::Matrix& matrix, const fesa::Vector3& value) {
  return {matrix(0U, 0U) * value.X() + matrix(0U, 1U) * value.Y() +
              matrix(0U, 2U) * value.Z(),
          matrix(1U, 0U) * value.X() + matrix(1U, 1U) * value.Y() +
              matrix(1U, 2U) * value.Z(),
          matrix(2U, 0U) * value.X() + matrix(2U, 1U) * value.Y() +
              matrix(2U, 2U) * value.Z()};
}

void ExpectZeroRecovery(const fesa::NonlinearMitc4Recovery& recovery) {
  EXPECT_DOUBLE_EQ(recovery.physical_strain_energy, 0.0);
  for (const auto& point : recovery.gauss_points) {
    for (const double value : point.generalized_strain) {
      EXPECT_DOUBLE_EQ(value, 0.0);
    }
    for (const double value : point.section_resultant) {
      EXPECT_DOUBLE_EQ(value, 0.0);
    }
    for (const auto& stress : point.stress) {
      for (const double value : stress.components) {
        EXPECT_DOUBLE_EQ(value, 0.0);
      }
    }
  }
}

fesa::NonlinearMitc4Shell Shell(const std::array<fesa::Node, 4>& nodes) {
  const auto result = fesa::NonlinearMitc4Shell::Create(
      NodePointers(nodes), Frames(), Section(), Material());
  EXPECT_TRUE(result.HasValue());
  return result.Value();
}

// NSNR-ELEM-KIN-001
TEST(NonlinearMitc4Shell, UndeformedStateHasZeroPhysicalRecovery) {
  const auto nodes = Nodes();
  const auto recovery = Shell(nodes).RecoverState(ZeroTrial());

  ASSERT_TRUE(recovery.HasValue());
  ExpectZeroRecovery(recovery.Value());
}

TEST(NonlinearMitc4Recovery, BendingResultantsMatchStressThicknessMoments) {
  const auto nodes = Nodes();
  auto trial = ZeroTrial();
  trial.rotation_vector = {{{0.12, -0.05, 0.0},
                            {-0.02, 0.13, 0.0},
                            {-0.12, 0.05, 0.0},
                            {0.02, -0.13, 0.0}}};
  for (const double thickness : {0.2, 0.5, 2.0}) {
    SCOPED_TRACE(thickness);
    auto section = Section();
    section.thickness = thickness;
    const auto shell = fesa::NonlinearMitc4Shell::Create(
        NodePointers(nodes), Frames(), section, Material());
    ASSERT_TRUE(shell.HasValue());
    const auto recovery = shell.Value().RecoverState(trial);
    ASSERT_TRUE(recovery.HasValue());
    for (const auto& point : recovery.Value().gauss_points) {
      for (std::size_t component = 0U; component < 3U; ++component) {
        // Flat reference geometry gives quadratic stress in z. Its odd
        // thickness moment is exact from independently sampled face stresses.
        const double moment = thickness * thickness / 12.0 *
                              (point.stress[2U].components[component] -
                               point.stress[0U].components[component]);
        ASSERT_GT(std::abs(moment), 0.0);
        EXPECT_NEAR(point.section_resultant[3U + component], moment,
                    1.0e-12 * std::abs(moment));
      }
    }
  }
}

// NSNR-ELEM-KIN-001A
TEST(NonlinearMitc4Shell, RejectsInvalidNodeAndInitialFrameContracts) {
  const auto nodes = Nodes();
  auto node_pointers = NodePointers(nodes);
  node_pointers[3U] = nullptr;
  const auto null_node = fesa::NonlinearMitc4Shell::Create(
      node_pointers, Frames(), Section(), Material());
  ASSERT_FALSE(null_node.HasValue());
  EXPECT_EQ(null_node.GetStatus().Diagnostics()[0U].code,
            "invalid-nonlinear-shell-geometry");

  auto frames = Frames();
  frames[2U].node_index = frames[1U].node_index;
  const auto duplicate_frame = fesa::NonlinearMitc4Shell::Create(
      NodePointers(nodes), frames, Section(), Material());
  ASSERT_FALSE(duplicate_frame.HasValue());
  EXPECT_EQ(duplicate_frame.GetStatus().Diagnostics()[0U].code,
            "invalid-nonlinear-shell-frame");

  frames = Frames();
  frames[0U].tangent_b = {1.0, 0.0, 0.0};
  const auto invalid_frame = fesa::NonlinearMitc4Shell::Create(
      NodePointers(nodes), frames, Section(), Material());
  ASSERT_FALSE(invalid_frame.HasValue());
  EXPECT_EQ(invalid_frame.GetStatus().Diagnostics()[0U].code,
            "invalid-nonlinear-shell-frame");
}

// NSNR-ELEM-KIN-002
TEST(NonlinearMitc4Shell, UsesFixedInitialReferenceGeometry) {
  const auto nodes = Nodes();
  auto trial = ZeroTrial();
  for (auto& displacement : trial.displacement) {
    displacement = {0.2, -0.3, 0.4};
  }

  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
  ASSERT_TRUE(recovery.HasValue());
  ExpectZeroRecovery(recovery.Value());
}

// NSNR-ELEM-KIN-003
TEST(NonlinearMitc4Shell, PreservesUnnormalizedInterpolatedDirector) {
  const auto nodes = Nodes();
  auto trial = ZeroTrial();
  trial.rotation_vector[0U] = {0.8, 0.0, 0.0};
  trial.rotation_vector[1U] = {-0.8, 0.0, 0.0};
  trial.rotation_vector[2U] = {-0.8, 0.0, 0.0};
  trial.rotation_vector[3U] = {0.8, 0.0, 0.0};

  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
  ASSERT_TRUE(recovery.HasValue());
  EXPECT_GT(recovery.Value().physical_strain_energy, 0.0);
  const double expected_membrane =
      0.01 * std::sin(0.8) * std::sin(0.8) / 6.0;
  EXPECT_NEAR(recovery.Value().gauss_points[0U].generalized_strain[0U],
              expected_membrane, 1.0e-12);
}

// NSNR-ELEM-KIN-004
TEST(NonlinearMitc4Shell, RejectsOpposedOrDegenerateCurrentDirectors) {
  const auto nodes = Nodes();
  auto trial = ZeroTrial();
  trial.rotation_vector[0U] = {std::acos(-1.0) - 1.0e-5, 0.0, 0.0};
  trial.rotation_vector[1U] = {-(std::acos(-1.0) - 1.0e-5), 0.0, 0.0};

  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
  ASSERT_FALSE(recovery.HasValue());
  EXPECT_EQ(recovery.GetStatus().Diagnostics()[0U].code,
            "invalid-nonlinear-shell-director");
}

// NSNR-ELEM-KIN-005
TEST(NonlinearMitc4Shell, MatchesLinearMitc4RecoveryInSmallRotationLimit) {
  const auto nodes = Nodes();
  const auto linear = fesa::Mitc4Shell::Create(
      NodePointers(nodes),
      std::array<Vector3, 4>{{{0.0, 0.0, 1.0}, {0.0, 0.0, 1.0},
                               {0.0, 0.0, 1.0}, {0.0, 0.0, 1.0}}},
      Section(), Material());
  ASSERT_TRUE(linear.HasValue());

  auto trial = ZeroTrial();
  fesa::Vector linear_dofs{24U};
  for (std::size_t node = 0U; node < 4U; ++node) {
    const double x = nodes[node].coordinates[0U];
    const double y = nodes[node].coordinates[1U];
    trial.displacement[node] = {1.0e-7 * x, -2.0e-7 * y, 3.0e-7 * x * y};
    trial.rotation_vector[node] = {-2.0e-7 * y, 1.0e-7 * x, 0.0};
    for (std::size_t component = 0U; component < 3U; ++component) {
      linear_dofs[6U * node + component] = trial.displacement[node][component];
      linear_dofs[6U * node + 3U + component] =
          trial.rotation_vector[node][component];
    }
  }

  const auto nonlinear_recovery = Shell(nodes).EvaluateKinematics(trial);
  const auto linear_recovery = linear.Value().RecoverPhysical(linear_dofs);
  ASSERT_TRUE(nonlinear_recovery.HasValue());
  ASSERT_TRUE(linear_recovery.HasValue());
  for (std::size_t point = 0U; point < 4U; ++point) {
    for (std::size_t component = 0U; component < 8U; ++component) {
      EXPECT_NEAR(nonlinear_recovery.Value().gauss_points[point]
                      .generalized_strain[component],
                  linear_recovery.Value().points[point]
                      .generalized_strain[component],
                  1.0e-11);
    }
  }
}

// NSNR-ELEM-KIN-006
TEST(NonlinearMitc4Shell, FiniteRigidMotionIsObjective) {
  const auto nodes = Nodes();
  const fesa::Vector3 rotation{0.3, -0.2, 0.1};
  const fesa::Matrix rigid_rotation = fesa::ExpSkew(rotation);
  const fesa::Vector3 translation{0.7, -0.4, 1.2};
  auto trial = ZeroTrial();
  for (std::size_t node = 0U; node < 4U; ++node) {
    const fesa::Vector3 initial{nodes[node].coordinates};
    trial.displacement[node] =
        Multiply(rigid_rotation, initial) + translation - initial;
    trial.rotation_vector[node] = rotation;
  }

  const auto recovery = Shell(nodes).EvaluateKinematics(trial);
  ASSERT_TRUE(recovery.HasValue());
  for (const auto& point : recovery.Value().gauss_points) {
    for (const double value : point.generalized_strain) {
      EXPECT_NEAR(value, 0.0, 1.0e-12);
    }
    for (const auto& stress : point.stress) {
      for (const double value : stress.components) {
        EXPECT_NEAR(value, 0.0, 1.0e-12);
      }
    }
  }
  EXPECT_NEAR(recovery.Value().physical_strain_energy, 0.0, 1.0e-12);
}

// NSNR-ELEM-KIN-007 and NSNR-ELEM-KIN-008
TEST(NonlinearMitc4Recovery, EmitsOnlyFixedPhysicalGaussAndStressRows) {
  const auto nodes = Nodes();
  const auto recovery = Shell(nodes).RecoverState(ZeroTrial());
  ASSERT_TRUE(recovery.HasValue());
  const double gauss = 1.0 / std::sqrt(3.0);
  constexpr std::array<fesa::ShellMidsurfaceLocation, 4> locations{
      fesa::ShellMidsurfaceLocation::kGp1, fesa::ShellMidsurfaceLocation::kGp2,
      fesa::ShellMidsurfaceLocation::kGp3, fesa::ShellMidsurfaceLocation::kGp4};
  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};

  for (std::size_t point = 0U; point < recovery.Value().gauss_points.size();
       ++point) {
    const auto& row = recovery.Value().gauss_points[point];
    EXPECT_EQ(row.location, locations[point]);
    EXPECT_EQ(row.natural_coordinates[0U], point == 0U || point == 3U ? -gauss
                                                                       : gauss);
    EXPECT_EQ(row.natural_coordinates[1U], point < 2U ? -gauss : gauss);
    for (std::size_t stress = 0U; stress < row.stress.size(); ++stress) {
      EXPECT_EQ(row.stress[stress].position, positions[stress]);
      EXPECT_EQ(row.stress[stress].zeta, zeta[stress]);
      EXPECT_EQ(row.stress[stress].components.size(), 3U);
    }
  }
}

}  // namespace