#ifndef FESA_ELEMENTS_MITC4_SHELL_H_
#define FESA_ELEMENTS_MITC4_SHELL_H_

#include <array>
#include <string>

#include "fesa/core/status.h"
#include "fesa/math/matrix.h"
#include "fesa/math/vector.h"
#include "fesa/model/model_types.h"

namespace fesa {

/// @brief Stores bilinear shape values and natural-coordinate derivatives.
struct Mitc4ShapeFunctions {
  std::array<double, 4> values;
  std::array<double, 4> xi_derivatives;
  std::array<double, 4> eta_derivatives;
};

/// @brief Stores a right-handed local shell frame at one location.
struct Mitc4LocalFrame {
  std::array<double, 3> e1;
  std::array<double, 3> e2;
  std::array<double, 3> e3;
};

/// @brief Stores canonical MITC4 covariant shear interpolation weights.
struct Mitc4TyingWeights {
  std::array<double, 2> xi_zeta;
  std::array<double, 2> eta_zeta;
};

/// @brief Stores one fixed 2-by-2-by-2 integration point and weight.
struct Mitc4QuadraturePoint {
  std::array<double, 3> natural_coordinates;
  double weight;
};

/// @brief Separates physical, drilling, and stabilized stiffness matrices.
struct Mitc4Stiffness {
  Matrix physical_local20;
  Matrix physical_global24;
  Matrix drilling_global24;
  Matrix stabilized_global24;
  double drilling_stiffness;
};

/// @brief Stores physical shell recovery at one midsurface Gauss location.
struct Mitc4PhysicalRecoveryPoint {
  std::array<double, 2> natural_coordinates;
  Mitc4LocalFrame local_frame;
  std::array<double, 8> generalized_strain;
  std::array<double, 8> section_resultant;
  std::array<std::array<double, 3>, 3> in_plane_stress;
};

/// @brief Stores physical-only recovery rows and strain energy.
struct Mitc4PhysicalRecovery {
  std::array<Mitc4PhysicalRecoveryPoint, 4> points;
  double strain_energy;
};

/// @brief Implements the approved small-rotation FESA-MITC4 shell kernel.
/// @note Physical and numerical drilling contributions remain separate.
class Mitc4Shell {
 public:
  /// @brief Creates a validated shell kernel from four non-owning node
  /// pointers.
  /// @param nodes Node pointers valid for the duration of this call.
  /// @param initial_directors Validated unit initial directors in node order.
  /// @param section Centered constant-thickness shell section.
  /// @param material Supported isotropic elastic material.
  /// @return A validated shell or a structured model failure.
  static Result<Mitc4Shell> Create(
      std::array<const Node*, 4> nodes,
      std::array<std::array<double, 3>, 4> initial_directors,
      const ShellSection& section, const LinearElasticMaterial& material);

  /// @brief Evaluates bilinear shape functions and derivatives.
  static Mitc4ShapeFunctions ShapeFunctions(double xi, double eta) noexcept;

  /// @brief Evaluates the canonical edge-midpoint tying weights.
  static Mitc4TyingWeights TyingWeights(double xi, double eta) noexcept;

  /// @brief Returns the fixed 2-by-2-by-2 quadrature inventory.
  static const std::array<Mitc4QuadraturePoint, 8>& VolumeQuadrature() noexcept;

  /// @brief Evaluates the right-handed local frame at a midsurface location.
  [[nodiscard]] Mitc4LocalFrame LocalFrame(double xi, double eta) const;

  /// @brief Returns the physical 24-to-20 transformation.
  [[nodiscard]] Matrix PhysicalTransformation20() const;

  /// @brief Returns the numerical drilling 24-to-4 transformation.
  [[nodiscard]] Matrix DrillingTransformation4() const;

  /// @brief Evaluates the direct five-component physical strain operator.
  [[nodiscard]] Matrix DirectStrainDisplacement20(double xi, double eta,
                                                  double zeta) const;

  /// @brief Evaluates all four canonical covariant tying shear samples.
  [[nodiscard]] Matrix CovariantTyingShearSamples20() const;

  /// @brief Evaluates the MITC-projected five-component strain operator.
  [[nodiscard]] Matrix StrainDisplacement20(double xi, double eta,
                                            double zeta) const;

  /// @brief Returns the isotropic in-plane plane-stress matrix.
  [[nodiscard]] Matrix PlaneStressConstitutive() const;

  /// @brief Returns the five-component plane-stress and shear matrix.
  [[nodiscard]] Matrix MaterialConstitutive5() const;

  /// @brief Returns the centered membrane section matrix.
  [[nodiscard]] Matrix MembraneSectionMatrix() const;

  /// @brief Returns the centered bending section matrix.
  [[nodiscard]] Matrix BendingSectionMatrix() const;

  /// @brief Returns the corrected transverse-shear section matrix.
  [[nodiscard]] Matrix TransverseShearSectionMatrix() const;

  /// @brief Computes physical, drilling, and stabilized stiffness matrices.
  /// @return Finite stiffness matrices or a structured model failure.
  [[nodiscard]] Result<Mitc4Stiffness> Stiffness() const;

  /// @brief Recovers physical shell quantities without drilling results.
  /// @param global_element_displacement24 Global element DOFs in node order.
  /// @return Physical recovery rows or a structured model failure.
  [[nodiscard]] Result<Mitc4PhysicalRecovery> RecoverPhysical(
      const Vector& global_element_displacement24) const;

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

  /// @brief Stores covariant, reciprocal, frame, and Jacobian data at one
  /// point.
  struct GeometryData {
    std::array<Vector3, 3> covariant;
    std::array<Vector3, 3> reciprocal;
    Mitc4LocalFrame frame;
    double jacobian;
  };

  /// @brief Stores validated shell geometry and constitutive state.
  Mitc4Shell(std::array<Vector3, 4> coordinates,
             std::array<Vector3, 4> directors, std::array<Vector3, 4> tangent_a,
             std::array<Vector3, 4> tangent_b, Vector3 normal_candidate,
             double thickness, double youngs_modulus, double poisson_ratio,
             SourceLocation source_location, std::string identity);

  /// @brief Evaluates a finite positive Jacobian and right-handed frame.
  bool EvaluateGeometry(double xi, double eta, double zeta,
                        GeometryData& result) const noexcept;

  /// @brief Evaluates displacement-basis derivatives in covariant directions.
  std::array<std::array<Vector3, 3>, 20> BasisDerivatives(
      double xi, double eta, double zeta) const noexcept;

  /// @brief Builds direct or tied strain without changing projection order.
  Matrix StrainDisplacement(double xi, double eta, double zeta,
                            const Matrix* tying_samples) const;

  std::array<Vector3, 4> coordinates_;
  std::array<Vector3, 4> directors_;
  std::array<Vector3, 4> tangent_a_;
  std::array<Vector3, 4> tangent_b_;
  Vector3 normal_candidate_;
  double thickness_;
  double youngs_modulus_;
  double poisson_ratio_;
  SourceLocation source_location_;
  std::string identity_;
};

}  // namespace fesa

#endif  // FESA_ELEMENTS_MITC4_SHELL_H_