#include "fesa/elements/thermal_shell.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <utility>

#include "thermal_shell_geometry.h"

namespace fesa {
namespace {
using thermal_internal::Interval;

/// @brief Measures the full endpoint deviation from a computed coefficient.
double ErrorBound(Interval interval, double value) {
  return thermal_internal::Up(std::max(std::abs(interval.lower - value),
                                       std::abs(interval.upper - value)));
}

}  // namespace

ThermalShell::ThermalShell(SourceEntityId source,
                           std::array<EntityIndex, 4> node_indices,
                           SourceLocation location, double conductivity)
    : source_{std::move(source)},
      node_indices_{node_indices},
      location_{std::move(location)},
      conductivity_value_{conductivity} {}

Result<ThermalShell> ThermalShell::Create(
    SourceEntityId source, std::array<EntityIndex, 4> node_indices,
    const std::array<Vector3, 4>& coordinates, double conductivity,
    double density, double specific_heat, double thickness,
    SourceLocation location) {
  using namespace thermal_internal;
  for (const double property :
       {conductivity, density, specific_heat, thickness}) {
    if (!(property > 0) || !std::isfinite(property)) {
      return Result<ThermalShell>::Failure(
          GeometryFailure(source, location, "invalid-thermal-shell-material",
                          "Thermal properties must be finite and positive."));
    }
  }
  auto sorted = node_indices;
  std::sort(sorted.begin(), sorted.end());
  if (std::adjacent_find(sorted.begin(), sorted.end()) != sorted.end()) {
    return Result<ThermalShell>::Failure(
        GeometryFailure(source, location, "invalid-thermal-shell-geometry",
                        "Patch contains repeated source nodes."));
  }
  const auto geometry = BuildGeometry(coordinates, source, location);
  if (!geometry.HasValue())
    return Result<ThermalShell>::Failure(geometry.GetStatus());
  ThermalShell kernel{source, node_indices, location, conductivity};
  Matrix mass{4, 4};
  Matrix surface{4, 4};
  std::array<std::array<Interval, 4>, 4> mass_interval{};
  std::array<std::array<Interval, 4>, 4> surface_interval{};
  const double gauss = std::sqrt(3. / 5.);
  const Interval gauss_interval = Sqrt(Interval{3} / Interval{5});
  const std::array<double, 3> points{{-gauss, 0., gauss}};
  const std::array<Interval, 3> point_intervals{
      {Interval{0} - gauss_interval, Interval{0}, gauss_interval}};
  const std::array<double, 3> weights{{5. / 9., 8. / 9., 5. / 9.}};
  const std::array<Interval, 3> weight_intervals{{Interval{5} / Interval{9},
                                                  Interval{8} / Interval{9},
                                                  Interval{5} / Interval{9}}};
  for (std::size_t x = 0; x < 3; ++x) {
    for (std::size_t y = 0; y < 3; ++y) {
      const auto result = EvaluateSurface(geometry.Value(), points[x],
                                          points[y], point_intervals[x],
                                          point_intervals[y], source, location);
      if (!result.HasValue())
        return Result<ThermalShell>::Failure(result.GetStatus());
      const auto& point = result.Value();
      const double weight = weights[x] * weights[y] * point.jacobian;
      const Interval weight_interval =
          weight_intervals[x] * weight_intervals[y] * point.jacobian_interval;
      for (std::size_t i = 0; i < 4; ++i) {
        for (std::size_t j = 0; j < 4; ++j) {
          mass(i, j) += point.shape[i] * point.shape[j] * weight;
          surface(i, j) += point.gradient[i].Dot(point.gradient[j]) * weight;
          mass_interval[i][j] = mass_interval[i][j] +
                                point.shape_interval[i] *
                                    point.shape_interval[j] * weight_interval;
          surface_interval[i][j] =
              surface_interval[i][j] +
              Dot(point.gradient_interval[i], point.gradient_interval[j]) *
                  weight_interval;
        }
      }
    }
  }
  Matrix z0{5, 5};
  Matrix z1{5, 5};
  std::array<std::array<Interval, 5>, 5> z0_interval{};
  std::array<std::array<Interval, 5>, 5> z1_interval{};
  constexpr double kZ0[3][3] = {{4, 2, -1}, {2, 16, 2}, {-1, 2, 4}};
  constexpr double kZ1[3][3] = {{7, -8, 1}, {-8, 16, -8}, {1, -8, 7}};
  for (const std::size_t first : {0U, 2U}) {
    for (std::size_t a = 0; a < 3; ++a) {
      for (std::size_t b = 0; b < 3; ++b) {
        z0(first + a, first + b) += kZ0[a][b] / 30.;
        z1(first + a, first + b) += kZ1[a][b] / 3.;
        z0_interval[first + a][first + b] = z0_interval[first + a][first + b] +
                                            Interval{kZ0[a][b]} / Interval{30};
        z1_interval[first + a][first + b] = z1_interval[first + a][first + b] +
                                            Interval{kZ1[a][b]} / Interval{3};
      }
    }
  }
  for (std::size_t i = 0; i < 4; ++i) {
    for (std::size_t j = 0; j < 4; ++j) {
      for (std::size_t a = 0; a < 5; ++a) {
        for (std::size_t b = 0; b < 5; ++b) {
          const std::size_t p = 5 * i + a;
          const std::size_t q = 5 * j + b;
          kernel.capacity_(p, q) =
              density * specific_heat * (thickness / 2) * mass(i, j) * z0(a, b);
          kernel.conductivity_(p, q) =
              conductivity * ((thickness / 2) * surface(i, j) * z0(a, b) +
                              (2 / thickness) * mass(i, j) * z1(a, b));
          const Interval capacity = Interval{density} *
                                    Interval{specific_heat} *
                                    (Interval{thickness} / Interval{2}) *
                                    mass_interval[i][j] * z0_interval[a][b];
          const Interval conduction =
              Interval{conductivity} *
              ((Interval{thickness} / Interval{2}) * surface_interval[i][j] *
                   z0_interval[a][b] +
               (Interval{2} / Interval{thickness}) * mass_interval[i][j] *
                   z1_interval[a][b]);
          kernel.capacity_error_(p, q) =
              ErrorBound(capacity, kernel.capacity_(p, q));
          kernel.conductivity_error_(p, q) =
              ErrorBound(conduction, kernel.conductivity_(p, q));
          if (!capacity.IsFinite() || !conduction.IsFinite() ||
              !std::isfinite(kernel.capacity_(p, q)) ||
              !std::isfinite(kernel.conductivity_(p, q)) ||
              !std::isfinite(kernel.capacity_error_(p, q)) ||
              !std::isfinite(kernel.conductivity_error_(p, q)) ||
              (p == q && (!(kernel.capacity_(p, q) > 0) ||
                          !(kernel.conductivity_(p, q) > 0)))) {
            return Result<ThermalShell>::Failure(GeometryFailure(
                source, location, "nonfinite-thermal-shell-operator",
                "Local matrix/enclosure is not finite or a diagonal lost "
                "positivity."));
          }
        }
      }
    }
  }
  const double recovery = 1 / std::sqrt(3.);
  const Interval recovery_interval = Interval{1} / Sqrt(Interval{3});
  constexpr std::array<double, 4> kXi{{-1, 1, 1, -1}};
  constexpr std::array<double, 4> kEta{{-1, -1, 1, 1}};
  constexpr std::array<double, 6> kLocal{{-1, 0, 1, -1, 0, 1}};
  for (std::size_t ip = 0; ip < 4; ++ip) {
    const auto result = EvaluateSurface(
        geometry.Value(), kXi[ip] * recovery, kEta[ip] * recovery,
        Interval{kXi[ip]} * recovery_interval,
        Interval{kEta[ip]} * recovery_interval, source, location);
    if (!result.HasValue())
      return Result<ThermalShell>::Failure(result.GetStatus());
    const auto& point = result.Value();
    kernel.frames_[ip] = point.axes;
    for (std::size_t trace = 0; trace < 6; ++trace) {
      const double s = kLocal[trace];
      const std::array<double, 3> shape{
          {.5 * s * (s - 1), 1 - s * s, .5 * s * (s + 1)}};
      const std::array<double, 3> derivative{{2 * s - 1, -4 * s, 2 * s + 1}};
      const std::size_t first = trace < 3 ? 0 : 2;
      for (std::size_t i = 0; i < 4; ++i) {
        for (std::size_t a = 0; a < 3; ++a) {
          const Vector3 gradient =
              shape[a] * point.gradient[i] +
              (point.shape[i] * (2 / thickness) * derivative[a]) *
                  point.axes[2];
          if (!gradient.IsFinite()) {
            return Result<ThermalShell>::Failure(GeometryFailure(
                source, location, "nonfinite-thermal-shell-recovery",
                "Flux gradient is not finite."));
          }
          kernel.gradients_[6 * ip + trace][5 * i + first + a] = gradient;
        }
      }
    }
  }
  return Result<ThermalShell>::Success(std::move(kernel));
}

const SourceEntityId& ThermalShell::SourceId() const noexcept {
  return source_;
}
const std::array<EntityIndex, 4>& ThermalShell::NodeIndices() const noexcept {
  return node_indices_;
}
const Matrix& ThermalShell::Conductivity() const noexcept {
  return conductivity_;
}
const Matrix& ThermalShell::Capacity() const noexcept { return capacity_; }
const Matrix& ThermalShell::ConductivityError() const noexcept {
  return conductivity_error_;
}
const Matrix& ThermalShell::CapacityError() const noexcept {
  return capacity_error_;
}

Result<ThermalShellRecovery> ThermalShell::Recover(
    const std::array<double, 20>& temperature) const {
  for (const double value : temperature) {
    if (!std::isfinite(value)) {
      return Result<ThermalShellRecovery>::Failure(
          thermal_internal::GeometryFailure(
              source_, location_, "nonfinite-thermal-shell-temperature",
              "Recovery temperature is not finite."));
    }
  }
  ThermalShellRecovery recovery{};
  recovery.frames = frames_;
  for (std::size_t row = 0; row < 24; ++row) {
    Vector3 gradient;
    for (std::size_t p = 0; p < 20; ++p)
      gradient = gradient + temperature[p] * gradients_[row][p];
    recovery.global_flux[row] = -conductivity_value_ * gradient;
    if (!recovery.global_flux[row].IsFinite()) {
      return Result<ThermalShellRecovery>::Failure(
          thermal_internal::GeometryFailure(
              source_, location_, "nonfinite-thermal-shell-flux",
              "Recovered heat flux is not finite."));
    }
  }
  return Result<ThermalShellRecovery>::Success(recovery);
}

}  // namespace fesa