#ifndef FESA_ELEMENTS_EULER_BEAM_3D_H_
#define FESA_ELEMENTS_EULER_BEAM_3D_H_

#include <array>
#include <cstddef>
#include <string>
#include <vector>

#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 constant line-load components in the beam local frame.
struct ConstantLocalLineLoad {
  double px;
  double py;
  double pz;
  double mx;
};

/// @brief Stores one axial stress at a Gauss and section-point identity.
struct BeamStressPoint {
  int gauss_point;
  std::size_t section_point;
  double x1;
  double x2;
  double s11;
  std::string source;
};

/// @brief Stores distinct beam end-action, section, Gauss, and stress results.
struct BeamRecovery {
  std::array<std::array<double, 6>, 2> equilibrium_end_actions;
  std::array<std::array<double, 4>, 2> endpoint_section_resultants;
  std::array<std::array<double, 4>, 2> gauss_generalized_strains;
  std::array<std::array<double, 4>, 2> gauss_generalized_resultants;
  std::vector<BeamStressPoint> stress_points;
};

/// @brief Implements the approved two-node prismatic B33 Euler beam kernel.
/// @note Equation numbering and semantic element identity remain external.
class EulerBeam3D {
 public:
  /// @brief Creates a validated beam kernel and right-handed local frame.
  /// @param first_node First source node in the element connectivity.
  /// @param second_node Second source node in the element connectivity.
  /// @param section Supported general beam section and local first axis.
  /// @param material Supported isotropic elastic material.
  /// @return A validated beam or a structured model failure.
  static Result<EulerBeam3D> Create(const Node& first_node,
                                    const Node& second_node,
                                    const GeneralBeamSection& section,
                                    const LinearElasticMaterial& material);

  /// @brief Computes the 12-by-12 stiffness in local DOF order.
  /// @note Uses the approved two-point Gauss operation order.
  Matrix LocalStiffness() const;

  /// @brief Computes the stiffness in stable global element DOF order.
  Matrix GlobalStiffness() const;

  /// @brief Computes the formulation-only constant local line-load vector.
  /// @warning This kernel does not expose distributed loads through parser
  /// input.
  Vector LocalEquivalentLoad(const ConstantLocalLineLoad& load) const;

  /// @brief Recovers signed physical quantities at their distinct locations.
  /// @param global_element_displacement Twelve global element DOF values.
  /// @return Beam recovery rows in deterministic location order.
  BeamRecovery Recover(const Vector& global_element_displacement) const;

 private:
  /// @brief Stores already validated geometry, material, and section state.
  EulerBeam3D(double length, double youngs_modulus, double shear_modulus,
              double area, double iy, double iz, double torsional_constant,
              std::array<double, 9> rotation,
              std::vector<std::array<double, 2>> section_points);

  double length_;
  double youngs_modulus_;
  double shear_modulus_;
  double area_;
  double iy_;
  double iz_;
  double torsional_constant_;
  std::array<double, 9> rotation_;
  std::vector<std::array<double, 2>> section_points_;
};

}  // namespace fesa

#endif  // FESA_ELEMENTS_EULER_BEAM_3D_H_