#include "fesa/model/shell_geometry.h"

#include <algorithm>
#include <cmath>
#include <tuple>
#include <utility>

#include "fesa/math/vector3.h"

#include "math/dense_blas_internal.h"
#include "math/vector3_operations.h"

namespace fesa {
namespace {

using math_internal::DerivativeSum;
using math_internal::WeightedSum;

constexpr std::array<double, 4> kXiSigns{-1.0, 1.0, 1.0, -1.0};
constexpr std::array<double, 4> kEtaSigns{-1.0, -1.0, 1.0, 1.0};

struct ShapeData {
  std::array<double, 4> values;
  std::array<double, 4> xi_derivatives;
  std::array<double, 4> eta_derivatives;
};

struct ElementWork {
  std::array<Vector3, 4> coordinates;
  Vector3 normal;
  double area_weight;
};

ShapeData ShapeDataAt(double xi, double eta) {
  ShapeData data{};
  for (std::size_t node = 0U; node < 4U; ++node) {
    data.values[node] =
        0.25 * (1.0 + kXiSigns[node] * xi) * (1.0 + kEtaSigns[node] * eta);
    data.xi_derivatives[node] =
        0.25 * kXiSigns[node] * (1.0 + kEtaSigns[node] * eta);
    data.eta_derivatives[node] =
        0.25 * kEtaSigns[node] * (1.0 + kXiSigns[node] * xi);
  }
  return data;
}

Result<ShellGeometry> GeometryFailure(std::string code,
                                      const SourceLocation& location,
                                      std::string keyword, std::string identity,
                                      std::string message) {
  return Result<ShellGeometry>::Failure(Status::Failure(
      FailureCategory::kModel,
      {{Severity::kError, std::move(code), location, std::move(keyword),
        std::move(identity), std::move(message)}}));
}

bool SameCoordinates(const Vector3& left, const Vector3& right) {
  return left == right;
}

double Orientation(const Vector3& first, const Vector3& second,
                   const Vector3& third, const Vector3& normal) {
  return (second - first).Cross(third - first).Dot(normal);
}

bool HasOppositeSigns(double first, double second) {
  return (first < 0.0 && second > 0.0) || (first > 0.0 && second < 0.0);
}

/// @brief Detects a proper projected edge crossing without a tolerance clamp.
bool SegmentsProperlyIntersect(const Vector3& first_start,
                               const Vector3& first_end,
                               const Vector3& second_start,
                               const Vector3& second_end,
                               const Vector3& normal) {
  const double first_side_start =
      Orientation(first_start, first_end, second_start, normal);
  const double first_side_end =
      Orientation(first_start, first_end, second_end, normal);
  const double second_side_start =
      Orientation(second_start, second_end, first_start, normal);
  const double second_side_end =
      Orientation(second_start, second_end, first_end, normal);
  return HasOppositeSigns(first_side_start, first_side_end) &&
         HasOppositeSigns(second_side_start, second_side_end);
}

/// @brief Orders incident shell contributions by stable source identity.
bool SourceIdentityLess(const Mitc4ShellDefinition& left,
                        std::size_t left_index,
                        const Mitc4ShellDefinition& right,
                        std::size_t right_index) {
  return std::tie(left.source_id.instance_name, left.source_id.source_label,
                  left.source_id.source_label_text, left_index) <
         std::tie(right.source_id.instance_name, right.source_id.source_label,
                  right.source_id.source_label_text, right_index);
}

}  // namespace

const std::array<ShellGeometryValidationPoint, 17>&
ShellGeometryValidationPoints() noexcept {
  static const std::array<ShellGeometryValidationPoint, 17> points = [] {
    const double g = 1.0 / std::sqrt(3.0);
    return std::array<ShellGeometryValidationPoint, 17>{
        ShellGeometryValidationPoint{
            ShellGeometryPointKind::kCenter, 0U, {0.0, 0.0, 0.0}},
        {ShellGeometryPointKind::kStiffness, 0U, {-g, -g, -g}},
        {ShellGeometryPointKind::kStiffness, 1U, {-g, -g, g}},
        {ShellGeometryPointKind::kStiffness, 2U, {g, -g, -g}},
        {ShellGeometryPointKind::kStiffness, 3U, {g, -g, g}},
        {ShellGeometryPointKind::kStiffness, 4U, {g, g, -g}},
        {ShellGeometryPointKind::kStiffness, 5U, {g, g, g}},
        {ShellGeometryPointKind::kStiffness, 6U, {-g, g, -g}},
        {ShellGeometryPointKind::kStiffness, 7U, {-g, g, g}},
        {ShellGeometryPointKind::kTying, 0U, {0.0, -1.0, 0.0}},
        {ShellGeometryPointKind::kTying, 1U, {0.0, 1.0, 0.0}},
        {ShellGeometryPointKind::kTying, 2U, {-1.0, 0.0, 0.0}},
        {ShellGeometryPointKind::kTying, 3U, {1.0, 0.0, 0.0}},
        {ShellGeometryPointKind::kRecovery, 0U, {-g, -g, 0.0}},
        {ShellGeometryPointKind::kRecovery, 1U, {g, -g, 0.0}},
        {ShellGeometryPointKind::kRecovery, 2U, {g, g, 0.0}},
        {ShellGeometryPointKind::kRecovery, 3U, {-g, g, 0.0}}};
  }();
  return points;
}

Result<ShellGeometry> PreprocessShellGeometry(
    const std::vector<Node>& nodes,
    const std::vector<Mitc4ShellDefinition>& elements,
    const std::vector<ShellSection>& sections) {
  ShellGeometry geometry;
  geometry.element_data.reserve(elements.size());
  std::vector<ElementWork> work;
  work.reserve(elements.size());

  const double g = 1.0 / std::sqrt(3.0);
  const std::array<Vector3, 4> surface_gauss_points{
      Vector3{-g, -g, 0.0}, Vector3{g, -g, 0.0}, Vector3{g, g, 0.0},
      Vector3{-g, g, 0.0}};

  for (std::size_t element_index = 0U; element_index < elements.size();
       ++element_index) {
    const auto& element = elements[element_index];
    ElementWork current{};
    for (std::size_t local_node = 0U; local_node < 4U; ++local_node) {
      if (element.node_indices[local_node] >= nodes.size()) {
        return GeometryFailure(
            "invalid-shell-geometry", element.location, "ELEMENT",
            element.source_id.source_label_text,
            "Shell geometry references an unavailable internal node.");
      }
      current.coordinates[local_node] =
          Vector3{nodes[element.node_indices[local_node]].coordinates};
      if (!current.coordinates[local_node].IsFinite()) {
        return GeometryFailure(
            "invalid-shell-geometry", element.location, "ELEMENT",
            element.source_id.source_label_text,
            "Shell geometry contains a nonfinite source coordinate.");
      }
    }
    for (std::size_t first = 0U; first < 4U; ++first) {
      for (std::size_t second = first + 1U; second < 4U; ++second) {
        if (element.node_indices[first] == element.node_indices[second] ||
            SameCoordinates(current.coordinates[first],
                            current.coordinates[second])) {
          return GeometryFailure("invalid-shell-geometry", element.location,
                                 "ELEMENT", element.source_id.source_label_text,
                                 "Shell geometry contains duplicate nodes.");
        }
      }
    }

    const ShapeData center = ShapeDataAt(0.0, 0.0);
    const Vector3 center_xi =
        DerivativeSum(center.xi_derivatives, current.coordinates);
    const Vector3 center_eta =
        DerivativeSum(center.eta_derivatives, current.coordinates);
    const Vector3 center_cross = center_xi.Cross(center_eta);
    const double center_measure = center_cross.Norm();
    if (!center_xi.IsFinite() || !center_eta.IsFinite() ||
        !center_cross.IsFinite() || !std::isfinite(center_measure) ||
        !(center_measure > 0.0)) {
      return GeometryFailure(
          "invalid-shell-geometry", element.location, "ELEMENT",
          element.source_id.source_label_text,
          "Shell center has no finite nonzero normal candidate.");
    }
    current.normal = (1.0 / center_measure) * center_cross;

    if (SegmentsProperlyIntersect(
            current.coordinates[0], current.coordinates[1],
            current.coordinates[2], current.coordinates[3], current.normal) ||
        SegmentsProperlyIntersect(
            current.coordinates[1], current.coordinates[2],
            current.coordinates[3], current.coordinates[0], current.normal)) {
      return GeometryFailure("invalid-shell-geometry", element.location,
                             "ELEMENT", element.source_id.source_label_text,
                             "Shell boundary is self-intersecting in the "
                             "center-normal projection.");
    }

    double area_weight = 0.0;
    for (const auto& point : surface_gauss_points) {
      const ShapeData shape = ShapeDataAt(point[0], point[1]);
      const Vector3 tangent_xi =
          DerivativeSum(shape.xi_derivatives, current.coordinates);
      const Vector3 tangent_eta =
          DerivativeSum(shape.eta_derivatives, current.coordinates);
      const Vector3 area_vector = tangent_xi.Cross(tangent_eta);
      const double measure = area_vector.Norm();
      if (!tangent_xi.IsFinite() || !tangent_eta.IsFinite() ||
          !area_vector.IsFinite() || !std::isfinite(measure) ||
          !(measure > 0.0) || !(area_vector.Dot(current.normal) > 0.0)) {
        return GeometryFailure("invalid-shell-geometry", element.location,
                               "ELEMENT", element.source_id.source_label_text,
                               "Shell surface is zero-area or locally reversed "
                               "at a required point.");
      }
      area_weight += measure;
    }
    if (!std::isfinite(area_weight) || !(area_weight > 0.0)) {
      return GeometryFailure("invalid-shell-geometry", element.location,
                             "ELEMENT", element.source_id.source_label_text,
                             "Shell surface-area weight is nonfinite or zero.");
    }
    current.area_weight = area_weight;
    work.push_back(current);
    geometry.element_data.push_back({static_cast<EntityIndex>(element_index),
                                     current.normal.Components(), area_weight});
  }

  std::vector<std::vector<std::size_t>> incident(nodes.size());
  for (std::size_t element_index = 0U; element_index < elements.size();
       ++element_index) {
    for (const EntityIndex node_index : elements[element_index].node_indices) {
      incident[node_index].push_back(element_index);
    }
  }

  geometry.nodal_frames.reserve(nodes.size());
  std::vector<const ShellNodeInitialFrame*> frame_by_node(nodes.size(),
                                                          nullptr);
  for (std::size_t node_index = 0U; node_index < incident.size();
       ++node_index) {
    auto& node_incident = incident[node_index];
    if (node_incident.empty()) {
      continue;
    }
    std::sort(node_incident.begin(), node_incident.end(),
              [&elements](std::size_t left, std::size_t right) {
                return SourceIdentityLess(elements[left], left, elements[right],
                                          right);
              });
    for (std::size_t first = 0U; first < node_incident.size(); ++first) {
      for (std::size_t second = first + 1U; second < node_incident.size();
           ++second) {
        const double pair_dot = work[node_incident[first]].normal.Dot(
            work[node_incident[second]].normal);
        if (!std::isfinite(pair_dot) || !(pair_dot > 0.0)) {
          return GeometryFailure("opposed-incident-normal",
                                 nodes[node_index].location, "NODE",
                                 nodes[node_index].source_id.source_label_text,
                                 "Incident shell normal candidates do not "
                                 "share a positive orientation hemisphere.");
        }
      }
    }

    double maximum_weight = 0.0;
    for (const std::size_t element_index : node_incident) {
      maximum_weight =
          std::max(maximum_weight, work[element_index].area_weight);
    }
    Vector3 director_sum{};
    for (const std::size_t element_index : node_incident) {
      director_sum =
          director_sum + (work[element_index].area_weight / maximum_weight) *
                             work[element_index].normal;
    }
    const double director_norm = director_sum.Norm();
    if (!director_sum.IsFinite() || !std::isfinite(director_norm) ||
        !(director_norm > 0.0)) {
      return GeometryFailure(
          "invalid-shell-director", nodes[node_index].location, "NODE",
          nodes[node_index].source_id.source_label_text,
          "Area-weighted shell director is nonfinite or zero.");
    }
    const Vector3 director = (1.0 / director_norm) * director_sum;

    const std::array<Vector3, 3> global_axes{
        Vector3{1.0, 0.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{0.0, 0.0, 1.0}};
    std::size_t selected_axis = 0U;
    double selected_alignment = std::abs(global_axes[0].Dot(director));
    for (std::size_t axis = 1U; axis < global_axes.size(); ++axis) {
      const double alignment = std::abs(global_axes[axis].Dot(director));
      if (alignment < selected_alignment) {
        selected_alignment = alignment;
        selected_axis = axis;
      }
    }
    const Vector3 tangent_candidate =
        global_axes[selected_axis] -
        global_axes[selected_axis].Dot(director) * director;
    const double tangent_norm = tangent_candidate.Norm();
    if (!tangent_candidate.IsFinite() || !std::isfinite(tangent_norm) ||
        !(tangent_norm > 0.0)) {
      return GeometryFailure(
          "invalid-shell-director", nodes[node_index].location, "NODE",
          nodes[node_index].source_id.source_label_text,
          "Least-aligned-axis tangent frame construction failed.");
    }
    const Vector3 tangent_a = (1.0 / tangent_norm) * tangent_candidate;
    const Vector3 tangent_b = director.Cross(tangent_a);
    if (!tangent_b.IsFinite() || !(tangent_b.Norm() > 0.0)) {
      return GeometryFailure(
          "invalid-shell-director", nodes[node_index].location, "NODE",
          nodes[node_index].source_id.source_label_text,
          "Right-handed shell tangent frame construction failed.");
    }
    geometry.nodal_frames.push_back(
        {static_cast<EntityIndex>(node_index), director.Components(),
         tangent_a.Components(), tangent_b.Components()});
  }

  // Build lookup only after frame storage is complete so later code never
  // observes a pointer invalidated by vector growth.
  for (const auto& frame : geometry.nodal_frames) {
    frame_by_node[frame.node_index] = &frame;
  }

  for (std::size_t element_index = 0U; element_index < elements.size();
       ++element_index) {
    const auto& element = elements[element_index];
    if (element.section_index >= sections.size()) {
      return GeometryFailure("invalid-shell-jacobian", element.location,
                             "ELEMENT", element.source_id.source_label_text,
                             "Shell geometry cannot resolve its thickness for "
                             "Jacobian validation.");
    }
    const double thickness = sections[element.section_index].thickness;
    std::array<Vector3, 4> directors{};
    for (std::size_t local_node = 0U; local_node < 4U; ++local_node) {
      const auto* frame = frame_by_node[element.node_indices[local_node]];
      if (frame == nullptr) {
        return GeometryFailure("invalid-shell-director", element.location,
                               "ELEMENT", element.source_id.source_label_text,
                               "Shell element is missing a nodal director.");
      }
      directors[local_node] = Vector3{frame->director};
    }

    for (const auto& point : ShellGeometryValidationPoints()) {
      const double xi = point.natural_coordinates[0];
      const double eta = point.natural_coordinates[1];
      const double zeta = point.natural_coordinates[2];
      const ShapeData shape = ShapeDataAt(xi, eta);
      const Vector3 midsurface_xi =
          DerivativeSum(shape.xi_derivatives, work[element_index].coordinates);
      const Vector3 midsurface_eta =
          DerivativeSum(shape.eta_derivatives, work[element_index].coordinates);
      const Vector3 area_vector = midsurface_xi.Cross(midsurface_eta);
      const double surface_measure = area_vector.Norm();
      if (!midsurface_xi.IsFinite() || !midsurface_eta.IsFinite() ||
          !area_vector.IsFinite() || !std::isfinite(surface_measure) ||
          !(surface_measure > 0.0) ||
          !(area_vector.Dot(work[element_index].normal) > 0.0)) {
        return GeometryFailure("invalid-shell-geometry", element.location,
                               "ELEMENT", element.source_id.source_label_text,
                               "Shell surface basis is nonfinite, zero, or "
                               "reversed at a required point.");
      }

      const Vector3 director_xi = WeightedSum(shape.xi_derivatives, directors);
      const Vector3 director_eta =
          WeightedSum(shape.eta_derivatives, directors);
      const Vector3 director_value = WeightedSum(shape.values, directors);
      const Vector3 covariant_xi =
          midsurface_xi + 0.5 * thickness * zeta * director_xi;
      const Vector3 covariant_eta =
          midsurface_eta + 0.5 * thickness * zeta * director_eta;
      const Vector3 covariant_zeta = 0.5 * thickness * director_value;
      const double jacobian =
          covariant_xi.Dot(covariant_eta.Cross(covariant_zeta));
      if (!covariant_xi.IsFinite() || !covariant_eta.IsFinite() ||
          !covariant_zeta.IsFinite() || !std::isfinite(jacobian) ||
          !(jacobian > 0.0)) {
        return GeometryFailure(
            "invalid-shell-jacobian", element.location, "ELEMENT",
            element.source_id.source_label_text,
            "Shell Jacobian is nonfinite or nonpositive at a required point.");
      }
      bool finite_actions = std::isfinite(1.0 / jacobian) &&
          covariant_eta.Cross(covariant_zeta).IsFinite() &&
          covariant_zeta.Cross(covariant_xi).IsFinite() &&
          covariant_xi.Cross(covariant_eta).IsFinite();
      if (finite_actions) {
        Matrix coefficients{3U, 3U};
        for (std::size_t i = 0; i < 3U; ++i) {
          coefficients(i, 0) = covariant_xi[i];
          coefficients(i, 1) = covariant_eta[i];
          coefficients(i, 2) = covariant_zeta[i];
        }
        const auto factor = dense_blas_internal::DenseFactor::Lu(coefficients);
        finite_actions = factor.HasValue();
        Matrix rhs{3U, 1U};
        for (std::size_t axis = 0; finite_actions && axis < 3U; ++axis) {
          if (axis != 0U) rhs(axis - 1U, 0U) = 0.0;
          rhs(axis, 0U) = 1.0;
          finite_actions = factor.Value().Solve(rhs, true).HasValue();
        }
      }
      if (!finite_actions) {
        return GeometryFailure(
            "invalid-shell-jacobian", element.location, "ELEMENT",
            element.source_id.source_label_text,
            "Shell reciprocal basis is nonfinite at a required point.");
      }
    }
  }

  return Result<ShellGeometry>::Success(std::move(geometry));
}

}  // namespace fesa