#pragma once

#include "fesa/Math/Math.hpp"
#include "fesa/Util/Util.hpp"

#include <array>
#include <optional>
#include <string>
#include <utility>
#include <vector>

namespace fesa {

struct ShapeData {
  std::array<Real, 4> n{};
  std::array<Real, 4> dr{};
  std::array<Real, 4> ds{};
};

inline ShapeData shapeFunctions(Real r, Real s) {
  return { {
               0.25 * (1.0 - r) * (1.0 - s),
               0.25 * (1.0 + r) * (1.0 - s),
               0.25 * (1.0 + r) * (1.0 + s),
               0.25 * (1.0 - r) * (1.0 + s),
           },
           {
               -0.25 * (1.0 - s),
               0.25 * (1.0 - s),
               0.25 * (1.0 + s),
               -0.25 * (1.0 + s),
           },
           {
               -0.25 * (1.0 - r),
               -0.25 * (1.0 + r),
               0.25 * (1.0 + r),
               0.25 * (1.0 - r),
           } };
}

struct LocalBasis {
  Vec3 e1;
  Vec3 e2;
  Vec3 e3;
};

struct MITC4NaturalPoint {
  Real xi = 0.0;
  Real eta = 0.0;
};

struct MITC4TyingPoint {
  std::string label;
  MITC4NaturalPoint natural;
  std::array<LocalIndex, 2> edge_node_indices{};
};

inline std::array<MITC4NaturalPoint, 4> mitc4NodeNaturalCoordinates() {
  return {{{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}}};
}

inline std::array<MITC4TyingPoint, 4> mitc4TyingPoints() {
  return {{{"A", {0.0, -1.0}, {0, 1}},
           {"B", {-1.0, 0.0}, {0, 3}},
           {"C", {0.0, 1.0}, {3, 2}},
           {"D", {1.0, 0.0}, {1, 2}}}};
}

struct MITC4DirectorFrame {
  Vec3 v1;
  Vec3 v2;
  Vec3 vn;
};

struct MITC4MidsurfaceDerivatives {
  ShapeData shape;
  Vec3 g1;
  Vec3 g2;
};

struct MITC4Geometry {
  std::array<Vec3, 4> coordinates{};
  Real thickness = 0.0;
  ShapeData center_shape;
  Vec3 g1_center;
  Vec3 g2_center;
  Vec3 center_normal;
  std::array<MITC4DirectorFrame, 4> nodal_frames{};
  std::vector<Diagnostic> diagnostics;

  bool ok() const {
    return !hasError(diagnostics);
  }
};

struct MITC4IntegrationBasis {
  ShapeData shape;
  Vec3 g1;
  Vec3 g2;
  Vec3 g3;
  Real jacobian = 0.0;
  LocalBasis local;
  std::vector<Diagnostic> diagnostics;

  bool ok() const {
    return !hasError(diagnostics);
  }
};

inline Vec3 globalEX() {
  return {1.0, 0.0, 0.0};
}

inline Vec3 globalEY() {
  return {0.0, 1.0, 0.0};
}

inline Vec3 globalEZ() {
  return {0.0, 0.0, 1.0};
}

inline Diagnostic mitc4Diagnostic(std::string code, std::string message) {
  return makeDiagnostic(Severity::Error, std::move(code), std::move(message), "mitc4", "<element>", 0);
}

inline void appendDiagnostics(std::vector<Diagnostic>& target, const std::vector<Diagnostic>& source) {
  target.insert(target.end(), source.begin(), source.end());
}

inline MITC4MidsurfaceDerivatives mitc4MidsurfaceDerivatives(const std::array<Vec3, 4>& coordinates,
                                                             Real xi,
                                                             Real eta) {
  MITC4MidsurfaceDerivatives result;
  result.shape = shapeFunctions(xi, eta);
  for (std::size_t i = 0; i < 4; ++i) {
    result.g1 = result.g1 + result.shape.dr[i] * coordinates[i];
    result.g2 = result.g2 + result.shape.ds[i] * coordinates[i];
  }
  return result;
}

inline std::optional<Vec3> firstNormalizedCross(const std::array<Vec3, 3>& axes, const Vec3& vector, Real tolerance) {
  for (const Vec3& axis : axes) {
    auto candidate = normalizedIfValid(cross(axis, vector), tolerance);
    if (candidate) {
      return candidate;
    }
  }
  return std::nullopt;
}

inline std::optional<MITC4DirectorFrame> buildMITC4DirectorFrame(const Vec3& normal, Real tolerance) {
  auto v1 = normalizedIfValid(cross(globalEY(), normal), tolerance);
  if (!v1) {
    v1 = firstNormalizedCross({globalEZ(), globalEX(), globalEY()}, normal, tolerance);
  }
  if (!v1) {
    return std::nullopt;
  }
  auto v2 = normalizedIfValid(cross(normal, *v1), tolerance);
  if (!v2) {
    return std::nullopt;
  }
  return MITC4DirectorFrame{*v1, *v2, normal};
}

inline MITC4Geometry buildMITC4Geometry(const std::array<Vec3, 4>& coordinates,
                                        Real thickness,
                                        Real tolerance = 1.0e-12) {
  MITC4Geometry geometry;
  geometry.coordinates = coordinates;
  geometry.thickness = thickness;
  if (!isFinite(thickness) || thickness <= tolerance) {
    geometry.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-THICKNESS", "MITC4 shell thickness must be positive and finite"));
  }
  for (const Vec3& coordinate : coordinates) {
    if (!isFinite(coordinate)) {
      geometry.diagnostics.push_back(
          mitc4Diagnostic("FESA-MITC4-COORDINATE", "MITC4 element coordinates must be finite"));
      break;
    }
  }

  const auto center = mitc4MidsurfaceDerivatives(coordinates, 0.0, 0.0);
  geometry.center_shape = center.shape;
  geometry.g1_center = center.g1;
  geometry.g2_center = center.g2;
  const auto normal = normalizedIfValid(cross(center.g1, center.g2), tolerance);
  if (!normal) {
    geometry.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-SINGULAR-NORMAL", "MITC4 element center normal is near zero"));
    return geometry;
  }
  geometry.center_normal = *normal;

  const auto frame = buildMITC4DirectorFrame(*normal, tolerance);
  if (!frame) {
    geometry.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 nodal director basis could not be constructed"));
    return geometry;
  }
  geometry.nodal_frames.fill(*frame);
  return geometry;
}

inline MITC4IntegrationBasis computeMITC4IntegrationBasis(const MITC4Geometry& geometry,
                                                          Real xi,
                                                          Real eta,
                                                          Real zeta,
                                                          Real tolerance = 1.0e-12) {
  MITC4IntegrationBasis result;
  result.diagnostics = geometry.diagnostics;
  result.shape = shapeFunctions(xi, eta);
  for (std::size_t i = 0; i < 4; ++i) {
    const Vec3& coordinate = geometry.coordinates[i];
    const Vec3& normal = geometry.nodal_frames[i].vn;
    result.g1 = result.g1 + result.shape.dr[i] * coordinate +
                (0.5 * zeta * geometry.thickness * result.shape.dr[i]) * normal;
    result.g2 = result.g2 + result.shape.ds[i] * coordinate +
                (0.5 * zeta * geometry.thickness * result.shape.ds[i]) * normal;
    result.g3 = result.g3 + (0.5 * geometry.thickness * result.shape.n[i]) * normal;
  }

  result.jacobian = dot(cross(result.g1, result.g2), result.g3);
  if (!isFinite(result.jacobian) || std::fabs(result.jacobian) <= tolerance) {
    result.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-SINGULAR-JACOBIAN", "MITC4 element Jacobian is near zero"));
  }

  const auto e3 = normalizedIfValid(result.g3, tolerance);
  if (!e3) {
    result.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis normal is near zero"));
    return result;
  }
  auto e1 = normalizedIfValid(cross(result.g2, *e3), tolerance);
  if (!e1) {
    e1 = firstNormalizedCross({globalEY(), globalEZ(), globalEX()}, *e3, tolerance);
  }
  if (!e1) {
    result.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis tangent could not be constructed"));
    return result;
  }
  const auto e2 = normalizedIfValid(cross(*e3, *e1), tolerance);
  if (!e2) {
    result.diagnostics.push_back(
        mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis is not right-handed"));
    return result;
  }
  result.local = {*e1, *e2, *e3};
  return result;
}

inline LocalBasis computeLocalBasis(const std::array<Vec3, 4>& coordinates) {
  const MITC4Geometry geometry = buildMITC4Geometry(coordinates, 1.0);
  if (!geometry.ok()) {
    throw std::runtime_error("invalid MITC4 geometry");
  }
  const MITC4DirectorFrame& frame = geometry.nodal_frames[0];
  return {frame.v1, frame.v2, frame.vn};
}

}  // namespace fesa