#ifndef FESA_ELEMENTS_THERMAL_SHELL_H_
#define FESA_ELEMENTS_THERMAL_SHELL_H_

#include <array>

#include "fesa/core/source_identity.h"
#include "fesa/core/status.h"
#include "fesa/math/matrix.h"
#include "fesa/math/vector3.h"

namespace fesa {

/// @brief Holds global heat flux at four CCW surface points times six traces.
/// @note Trace order is -1+, -.5, 0-, 0+, .5, 1-. Frame columns are e1,e2,n.
struct ThermalShellRecovery {
  std::array<Vector3, 24> global_flux;
  std::array<std::array<Vector3, 3>, 4> frames;
};

/// @brief Owns checked local Q4 times two-segment P2 transient-heat operators.
/// @note Twenty DOFs are node-major, five stations bottom to top per node.
class ThermalShell final {
 public:
  /// @brief Builds finite consistent capacity/conduction and flux operators.
  /// @param coordinates Original midsurface Cartesian coordinates, without
  /// projection.
  /// @return A complete kernel or a source-identified geometry/material
  /// failure.
  static Result<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 = {});
  /// @brief Returns the original element identity.
  const SourceEntityId& SourceId() const noexcept;
  /// @brief Returns Domain node indices, not equation numbers.
  const std::array<EntityIndex, 4>& NodeIndices() const noexcept;
  /// @brief Returns the node-major conductivity matrix [power/temperature].
  const Matrix& Conductivity() const noexcept;
  /// @brief Returns the consistent capacity matrix [energy/temperature].
  const Matrix& Capacity() const noexcept;
  /// @brief Returns componentwise outward roundoff bounds, excluding quadrature
  /// error.
  const Matrix& ConductivityError() const noexcept;
  /// @brief Returns componentwise outward roundoff bounds, excluding quadrature
  /// error.
  const Matrix& CapacityError() const noexcept;
  /// @brief Evaluates q=-k*grad(T) without averaging the two midsurface traces.
  /// @return Finite global fluxes and right-handed frames or source-identified
  /// failure.
  Result<ThermalShellRecovery> Recover(
      const std::array<double, 20>& temperature) const;

 private:
  /// @brief Initializes an unpublished candidate whose operators are filled by
  /// Create.
  ThermalShell(SourceEntityId source, std::array<EntityIndex, 4> node_indices,
               SourceLocation location, double conductivity);
  SourceEntityId source_;
  std::array<EntityIndex, 4> node_indices_;
  SourceLocation location_;
  double conductivity_value_;
  Matrix conductivity_{20U, 20U};
  Matrix capacity_{20U, 20U};
  Matrix conductivity_error_{20U, 20U};
  Matrix capacity_error_{20U, 20U};
  std::array<std::array<Vector3, 20>, 24> gradients_{};
  std::array<std::array<Vector3, 3>, 4> frames_{};
};

}  // namespace fesa
#endif  // FESA_ELEMENTS_THERMAL_SHELL_H_