#include "fesa/elements/mitc4_shell.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>

#include "fesa/model/model_types.h"
#include "math/dense_blas_internal.h"
#include "math/finite_values.h"
#include "math/vector3_operations.h"
#include "elements/shell_kernel_values.h"

namespace fesa {

Mitc4ShellDefinition::Mitc4ShellDefinition(
    SourceEntityId source_id_value,
    const ShellSourceElementType source_type_value,
    std::array<EntityIndex, 4> node_indices_value,
    const EntityIndex material_index_value,
    const EntityIndex section_index_value, SourceLocation location_value)
    : source_id{std::move(source_id_value)},
      source_type{source_type_value},
      node_indices{std::move(node_indices_value)},
      material_index{material_index_value},
      section_index{section_index_value},
      location{std::move(location_value)},
      property_index_{section_index_value},
      node_indices_view_{node_indices.begin(), node_indices.end()} {}

ElementDefinitionKind Mitc4ShellDefinition::Kind() const noexcept {
  return ElementDefinitionKind::kMitc4Shell;
}

const SourceEntityId& Mitc4ShellDefinition::SourceId() const noexcept {
  return source_id;
}

std::string_view Mitc4ShellDefinition::SourceElementType() const noexcept {
  return source_type == ShellSourceElementType::kS4 ? "S4" : "S4R";
}

const std::vector<EntityIndex>& Mitc4ShellDefinition::NodeIndices()
    const noexcept {
  return node_indices_view_;
}

EntityIndex Mitc4ShellDefinition::MaterialIndex() const noexcept {
  return material_index;
}

EntityIndex Mitc4ShellDefinition::PropertyIndex() const noexcept {
  return property_index_;
}

void Mitc4ShellDefinition::SetPropertyIndex(
    const EntityIndex property_index) noexcept {
  property_index_ = property_index;
}

void Mitc4ShellDefinition::SynchronizeNodeIndices() {
  node_indices_view_.assign(node_indices.begin(), node_indices.end());
}

namespace {

using shell_kernel_internal::FrameComponent;
using shell_kernel_internal::LocalEngineeringComponents;
using shell_kernel_internal::ElementIdentity;
using math_internal::DerivativeSum;
using math_internal::WeightedSum;
using math_internal::IsFinite;

constexpr std::size_t kNodeCount = 4U;
constexpr std::size_t kPhysicalDofsPerNode = 5U;
constexpr std::size_t kGlobalDofsPerNode = 6U;
constexpr std::size_t kPhysicalDofCount = 20U;
constexpr std::size_t kGlobalDofCount = 24U;
constexpr double kShearCorrection = 5.0 / 6.0;
constexpr double kFrameTolerance = 1.0e-12;

bool SameCoordinates(const Vector3& left, const Vector3& right) {
  return left == right;
}

/// @brief Selects each director's least-aligned global axis deterministically
/// and projects it into the tangent plane.
std::array<Vector3, kNodeCount> NodalTangentsA(
    const std::array<Vector3, kNodeCount>& directors) {
  const std::array<Vector3, 3> global_axes{
      Vector3{1.0, 0.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{0.0, 0.0, 1.0}};
  std::array<Vector3, kNodeCount> tangents{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    std::size_t selected = 0U;
    double alignment = std::abs(global_axes[0].Dot(directors[node]));
    for (std::size_t axis = 1U; axis < global_axes.size(); ++axis) {
      const double candidate = std::abs(global_axes[axis].Dot(directors[node]));
      if (candidate < alignment) {
        selected = axis;
        alignment = candidate;
      }
    }
    const Vector3 tangent_candidate =
        global_axes[selected] -
        global_axes[selected].Dot(directors[node]) * directors[node];
    tangents[node] = (1.0 / tangent_candidate.Norm()) * tangent_candidate;
  }
  return tangents;
}

/// @brief Completes each right-handed nodal director frame.
std::array<Vector3, kNodeCount> NodalTangentsB(
    const std::array<Vector3, kNodeCount>& directors,
    const std::array<Vector3, kNodeCount>& tangent_a) {
  std::array<Vector3, kNodeCount> tangents{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    tangents[node] = directors[node].Cross(tangent_a[node]);
  }
  return tangents;
}

Result<Mitc4Shell> ModelFailure(std::string code,
                                const SourceLocation& location,
                                const std::string& identity,
                                std::string message) {
  return Result<Mitc4Shell>::Failure(Status::Failure(
      FailureCategory::kModel, {{Severity::kError, std::move(code), location,
                                 "*ELEMENT", identity, std::move(message)}}));
}

/// @brief Forms one physical covariant strain column with thickness stretch
/// excluded.
std::array<std::array<double, 3>, 3> CovariantStrainColumn(
    const std::array<Vector3, 3>& covariant,
    const std::array<Vector3, 3>& derivatives) {
  std::array<std::array<double, 3>, 3> strain{};
  for (std::size_t first = 0U; first < 3U; ++first) {
    for (std::size_t second = 0U; second < 3U; ++second) {
      strain[first][second] = 0.5 * (covariant[first].Dot(derivatives[second]) +
                                     covariant[second].Dot(derivatives[first]));
    }
  }
  // Thickness stretch is excluded from the five-component shell law.
  strain[2U][2U] = 0.0;
  return strain;
}

/// @brief Solves C^T E C = E_cov for the actual strain tensor.
std::array<std::array<double, 3>, 3> ReconstructCartesianStrain(
    const std::array<std::array<double, 3>, 3>& covariant_strain,
    const dense_blas_internal::DenseFactor& factor) {
  Matrix rhs{3U, 3U};
  for (std::size_t i = 0; i < 3U; ++i)
    for (std::size_t j = 0; j < 3U; ++j) rhs(i, j) = covariant_strain[i][j];
  const auto left = factor.Solve(rhs, true);
  if (!left.HasValue())
    throw std::invalid_argument{"MITC4 Cartesian strain is nonfinite."};
  const auto solved = factor.SolveRight(left.Value());
  if (!solved.HasValue())
    throw std::invalid_argument{"MITC4 Cartesian strain is nonfinite."};
  std::array<std::array<double, 3>, 3> tensor{};
  for (std::size_t i = 0; i < 3U; ++i)
    for (std::size_t j = 0; j < 3U; ++j) tensor[i][j] = solved.Value()(i, j);
  return tensor;
}

/// @brief Packs covariant columns without constructing reciprocal vectors.
Matrix CovariantMatrix(const std::array<Vector3, 3>& covariant) {
  Matrix coefficients{3U, 3U};
  for (std::size_t i = 0; i < 3U; ++i)
    for (std::size_t j = 0; j < 3U; ++j) coefficients(i, j) = covariant[j][i];
  return coefficients;
}

Matrix ScaledMatrix(const Matrix& source, double factor) {
  Matrix result{source};
  result.Scale(factor);
  return result;
}

/// @brief Lifts a local matrix through MKL and mirrors the upper triangle.
Matrix Congruence(const Matrix& local, const Matrix& transformation) {
  if (local.Rows() != local.Columns() ||
      local.Rows() != transformation.Rows()) {
    throw std::invalid_argument{
        "MITC4 congruence dimensions are incompatible."};
  }
  Matrix result = dense_blas_internal::MultiplyTransposeLeft(
      transformation, local.Multiply(transformation));
  for (std::size_t row = 0U; row < result.Rows(); ++row) {
    for (std::size_t column = row; column < result.Columns(); ++column) {
      result(column, row) = result(row, column);
    }
  }
  return result;
}

Result<Mitc4Stiffness> StiffnessFailure(const SourceLocation& location,
                                        const std::string& identity,
                                        std::string message) {
  return Result<Mitc4Stiffness>::Failure(
      Status::Failure(FailureCategory::kModel,
                      {{Severity::kError, "invalid-shell-stiffness", location,
                        "*ELEMENT", identity, std::move(message)}}));
}

Result<Mitc4Mass> MassFailure(const SourceLocation& location,
                              const std::string& identity,
                              std::string message) {
  return Result<Mitc4Mass>::Failure(
      Status::Failure(FailureCategory::kModel,
                      {{Severity::kError, "invalid-shell-mass", location,
                        "*ELEMENT", identity, std::move(message)}}));
}

/// @brief Returns S=-[director]_cross for theta cross director.
Matrix DirectorSkew(const Vector3& director) {
  Matrix result{3U, 3U};
  result(0U, 1U) = director.Z();
  result(0U, 2U) = -director.Y();
  result(1U, 0U) = -director.Z();
  result(1U, 2U) = director.X();
  result(2U, 0U) = director.Y();
  result(2U, 1U) = -director.X();
  return result;
}

Result<Mitc4PhysicalRecovery> RecoveryFailure(const SourceLocation& location,
                                              const std::string& identity,
                                              std::string message) {
  return Result<Mitc4PhysicalRecovery>::Failure(
      Status::Failure(FailureCategory::kModel,
                      {{Severity::kError, "invalid-shell-recovery", location,
                        "*ELEMENT", identity, std::move(message)}}));
}

/// @brief Creates a structured failure at the runtime element boundary.
template <class T>
Result<T> RuntimeFailure(const SourceLocation& location,
                         const SourceEntityId& source_id, std::string message) {
  return Result<T>::Failure(Status::Failure(
      FailureCategory::kModel,
      {{Severity::kError, "invalid-runtime-element", location, "*ELEMENT",
        source_id.source_label_text, std::move(message)}}));
}

}  // namespace

Result<Mitc4Shell> Mitc4Shell::Create(
    std::array<const Node*, 4> nodes,
    std::array<std::array<double, 3>, 4> initial_directors,
    const ShellSection& section, const LinearElasticMaterial& material) {
  const std::string identity = ElementIdentity(nodes);
  if (std::any_of(nodes.begin(), nodes.end(),
                  [](const Node* value) { return value == nullptr; })) {
    return ModelFailure("invalid-shell-geometry", section.location, identity,
                        "MITC4 creation requires four valid node references.");
  }

  std::array<Vector3, kNodeCount> coordinates{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    coordinates[node] = Vector3{nodes[node]->coordinates};
    if (!coordinates[node].IsFinite()) {
      return ModelFailure("invalid-shell-geometry", nodes[node]->location,
                          identity, "MITC4 node coordinates must be finite.");
    }
    for (std::size_t previous = 0U; previous < node; ++previous) {
      if (SameCoordinates(coordinates[node], coordinates[previous])) {
        return ModelFailure("invalid-shell-geometry", nodes[node]->location,
                            identity,
                            "MITC4 nodes must have distinct coordinates.");
      }
    }
  }

  const std::array<Vector3, kNodeCount> directors{
      Vector3{initial_directors[0U]}, Vector3{initial_directors[1U]},
      Vector3{initial_directors[2U]}, Vector3{initial_directors[3U]}};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const double director_norm = directors[node].Norm();
    if (!directors[node].IsFinite() || !std::isfinite(director_norm) ||
        std::abs(director_norm - 1.0) > kFrameTolerance) {
      return ModelFailure(
          "invalid-shell-director", nodes[node]->location, identity,
          "MITC4 initial directors must be finite unit vectors.");
    }
  }

  if (!std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
    return ModelFailure("invalid-shell-section", section.location, identity,
                        "MITC4 shell thickness must be finite and positive.");
  }
  if (!std::isfinite(material.youngs_modulus) ||
      !(material.youngs_modulus > 0.0) ||
      !std::isfinite(material.poisson_ratio) ||
      !(material.poisson_ratio > -1.0) || !(material.poisson_ratio < 0.5)) {
    return ModelFailure(
        "invalid-shell-material", material.location, identity,
        "MITC4 isotropic material requires finite E>0 and -1<nu<0.5.");
  }

  const auto center_shape = ShapeFunctions(0.0, 0.0);
  const Vector3 center_xi =
      DerivativeSum(center_shape.xi_derivatives, coordinates);
  const Vector3 center_eta =
      DerivativeSum(center_shape.eta_derivatives, coordinates);
  const Vector3 center_area = center_xi.Cross(center_eta);
  const double center_measure = center_area.Norm();
  if (!center_area.IsFinite() || !std::isfinite(center_measure) ||
      !(center_measure > 0.0)) {
    return ModelFailure(
        "invalid-shell-geometry", nodes[0]->location, identity,
        "MITC4 center surface basis must be finite and nonzero.");
  }
  const Vector3 normal_candidate = (1.0 / center_measure) * center_area;
  if (std::any_of(directors.begin(), directors.end(),
                  [&normal_candidate](const Vector3& director) {
                    return !(normal_candidate.Dot(director) > 0.0);
                  })) {
    return ModelFailure(
        "invalid-shell-director", nodes[0]->location, identity,
        "MITC4 directors must follow the source-order positive face.");
  }

  const auto tangent_a = NodalTangentsA(directors);
  const auto tangent_b = NodalTangentsB(directors, tangent_a);
  Mitc4Shell shell{coordinates,
                   directors,
                   tangent_a,
                   tangent_b,
                   normal_candidate,
                   section.thickness,
                   material.youngs_modulus,
                   material.poisson_ratio,
                   nodes[0]->location,
                   identity};

  const double shear_modulus =
      material.youngs_modulus / (2.0 * (1.0 + material.poisson_ratio));
  const double plane_stress_factor =
      material.youngs_modulus /
      (1.0 - material.poisson_ratio * material.poisson_ratio);
  const double thickness_cubed =
      section.thickness * section.thickness * section.thickness;
  const std::array<double, 6> derived{
      shear_modulus,
      plane_stress_factor,
      kShearCorrection * shear_modulus,
      plane_stress_factor * section.thickness,
      plane_stress_factor * thickness_cubed / 12.0,
      kShearCorrection * shear_modulus * section.thickness};
  if (std::any_of(derived.begin(), derived.end(), [](double value) {
        return !std::isfinite(value) || !(value > 0.0);
      })) {
    return ModelFailure(
        "invalid-shell-material", material.location, identity,
        "Derived MITC4 constitutive coefficients must be finite and positive.");
  }

  GeometryData geometry{};
  if (!shell.EvaluateGeometry(0.0, 0.0, 0.0, geometry)) {
    return ModelFailure("invalid-shell-jacobian", nodes[0]->location, identity,
                        "MITC4 center frame or Jacobian is invalid.");
  }
  for (const auto& point : VolumeQuadrature()) {
    if (!shell.EvaluateGeometry(point.natural_coordinates[0],
                                point.natural_coordinates[1],
                                point.natural_coordinates[2], geometry)) {
      return ModelFailure("invalid-shell-jacobian", nodes[0]->location,
                          identity,
                          "MITC4 quadrature frame or Jacobian is invalid.");
    }
  }
  constexpr std::array<Vector3, 4> tying_points{
      Vector3{0.0, -1.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{-1.0, 0.0, 0.0},
      Vector3{1.0, 0.0, 0.0}};
  for (const auto& point : tying_points) {
    if (!shell.EvaluateGeometry(point[0], point[1], point[2], geometry)) {
      return ModelFailure("invalid-shell-jacobian", nodes[0]->location,
                          identity,
                          "MITC4 tying-point frame or Jacobian is invalid.");
    }
  }

  return Result<Mitc4Shell>::Success(std::move(shell));
}

Mitc4ShapeFunctions Mitc4Shell::ShapeFunctions(double xi, double eta) noexcept {
  constexpr std::array<double, kNodeCount> xi_signs{-1.0, 1.0, 1.0, -1.0};
  constexpr std::array<double, kNodeCount> eta_signs{-1.0, -1.0, 1.0, 1.0};
  Mitc4ShapeFunctions shape{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    shape.values[node] =
        0.25 * (1.0 + xi_signs[node] * xi) * (1.0 + eta_signs[node] * eta);
    shape.xi_derivatives[node] =
        0.25 * xi_signs[node] * (1.0 + eta_signs[node] * eta);
    shape.eta_derivatives[node] =
        0.25 * eta_signs[node] * (1.0 + xi_signs[node] * xi);
  }
  return shape;
}

Mitc4TyingWeights Mitc4Shell::TyingWeights(double xi, double eta) noexcept {
  return {{(1.0 - eta) * 0.5, (1.0 + eta) * 0.5},
          {(1.0 - xi) * 0.5, (1.0 + xi) * 0.5}};
}

const std::array<Mitc4QuadraturePoint, 8>&
Mitc4Shell::VolumeQuadrature() noexcept {
  static const std::array<Mitc4QuadraturePoint, 8> points = [] {
    const double gauss = 1.0 / std::sqrt(3.0);
    return std::array<Mitc4QuadraturePoint, 8>{
        Mitc4QuadraturePoint{{-gauss, -gauss, -gauss}, 1.0},
        {{-gauss, -gauss, gauss}, 1.0},
        {{gauss, -gauss, -gauss}, 1.0},
        {{gauss, -gauss, gauss}, 1.0},
        {{gauss, gauss, -gauss}, 1.0},
        {{gauss, gauss, gauss}, 1.0},
        {{-gauss, gauss, -gauss}, 1.0},
        {{-gauss, gauss, gauss}, 1.0}};
  }();
  return points;
}

Mitc4LocalFrame Mitc4Shell::LocalFrame(double xi, double eta) const {
  GeometryData geometry{};
  if (!EvaluateGeometry(xi, eta, 0.0, geometry)) {
    throw std::invalid_argument{
        "MITC4 local frame is invalid at the requested point."};
  }
  return geometry.frame;
}

Matrix Mitc4Shell::PhysicalTransformation20() const {
  Matrix transformation{kPhysicalDofCount, kGlobalDofCount};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const std::size_t physical_offset = node * kPhysicalDofsPerNode;
    const std::size_t global_offset = node * kGlobalDofsPerNode;
    for (std::size_t component = 0U; component < 3U; ++component) {
      transformation(physical_offset + component, global_offset + component) =
          1.0;
      transformation(physical_offset + 3U, global_offset + 3U + component) =
          tangent_a_[node][component];
      transformation(physical_offset + 4U, global_offset + 3U + component) =
          tangent_b_[node][component];
    }
  }
  return transformation;
}

Matrix Mitc4Shell::DrillingTransformation4() const {
  Matrix transformation{kNodeCount, kGlobalDofCount};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const std::size_t global_offset = node * kGlobalDofsPerNode;
    for (std::size_t component = 0U; component < 3U; ++component) {
      transformation(node, global_offset + 3U + component) =
          directors_[node][component];
    }
  }
  return transformation;
}

Matrix Mitc4Shell::DirectStrainDisplacement20(double xi, double eta,
                                              double zeta) const {
  return StrainDisplacement(xi, eta, zeta, nullptr);
}

Matrix Mitc4Shell::CovariantTyingShearSamples20() const {
  constexpr std::array<Vector3, 4> points{
      Vector3{0.0, -1.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{-1.0, 0.0, 0.0},
      Vector3{1.0, 0.0, 0.0}};
  Matrix samples{4U, kPhysicalDofCount};
  for (std::size_t point = 0U; point < points.size(); ++point) {
    GeometryData geometry{};
    if (!EvaluateGeometry(points[point][0], points[point][1], 0.0, geometry)) {
      throw std::logic_error{"Validated MITC4 tying geometry became invalid."};
    }
    const auto derivatives =
        BasisDerivatives(points[point][0], points[point][1], 0.0);
    const std::size_t first = point < 2U ? 0U : 1U;
    for (std::size_t dof = 0U; dof < kPhysicalDofCount; ++dof) {
      const auto strain =
          CovariantStrainColumn(geometry.covariant, derivatives[dof]);
      samples(point, dof) = strain[first][2U];
    }
  }
  return samples;
}

Matrix Mitc4Shell::StrainDisplacement20(double xi, double eta,
                                        double zeta) const {
  const Matrix samples = CovariantTyingShearSamples20();
  return StrainDisplacement(xi, eta, zeta, &samples);
}

Matrix Mitc4Shell::PlaneStressConstitutive() const {
  const double factor =
      youngs_modulus_ / (1.0 - poisson_ratio_ * poisson_ratio_);
  Matrix constitutive{3U, 3U};
  constitutive(0U, 0U) = factor;
  constitutive(0U, 1U) = factor * poisson_ratio_;
  constitutive(1U, 0U) = factor * poisson_ratio_;
  constitutive(1U, 1U) = factor;
  constitutive(2U, 2U) = factor * (1.0 - poisson_ratio_) * 0.5;
  return constitutive;
}

Matrix Mitc4Shell::MaterialConstitutive5() const {
  Matrix constitutive{5U, 5U};
  const Matrix plane_stress = PlaneStressConstitutive();
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      constitutive(row, column) = plane_stress(row, column);
    }
  }
  const double shear_modulus = youngs_modulus_ / (2.0 * (1.0 + poisson_ratio_));
  constitutive(3U, 3U) = kShearCorrection * shear_modulus;
  constitutive(4U, 4U) = kShearCorrection * shear_modulus;
  return constitutive;
}

Matrix Mitc4Shell::MembraneSectionMatrix() const {
  return ScaledMatrix(PlaneStressConstitutive(), thickness_);
}

Matrix Mitc4Shell::BendingSectionMatrix() const {
  return ScaledMatrix(PlaneStressConstitutive(),
                      thickness_ * thickness_ * thickness_ / 12.0);
}

Matrix Mitc4Shell::TransverseShearSectionMatrix() const {
  const double shear_modulus = youngs_modulus_ / (2.0 * (1.0 + poisson_ratio_));
  Matrix result{2U, 2U};
  result(0U, 0U) = kShearCorrection * shear_modulus * thickness_;
  result(1U, 1U) = result(0U, 0U);
  return result;
}

Result<Mitc4Stiffness> Mitc4Shell::Stiffness() const try {
  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
  const Matrix constitutive = MaterialConstitutive5();
  const Matrix tying_samples = CovariantTyingShearSamples20();
  Matrix physical_local{kPhysicalDofCount, kPhysicalDofCount};
  for (const auto& point : VolumeQuadrature()) {
    GeometryData geometry{};
    if (!EvaluateGeometry(point.natural_coordinates[0],
                          point.natural_coordinates[1],
                          point.natural_coordinates[2], geometry)) {
      return StiffnessFailure(
          source_location_, identity_,
          "Validated MITC4 quadrature geometry became invalid.");
    }
    const Matrix strain = StrainDisplacement(
        point.natural_coordinates[0], point.natural_coordinates[1],
        point.natural_coordinates[2], &tying_samples);
    const Matrix integrand = dense_blas_internal::MultiplyTransposeLeft(
        strain, constitutive.Multiply(strain));
    for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) {
      for (std::size_t column = row; column < kPhysicalDofCount; ++column) {
        const double contribution =
            integrand(row, column) * geometry.jacobian * point.weight;
        physical_local(row, column) += contribution;
        if (row != column) {
          physical_local(column, row) += contribution;
        }
      }
    }
  }

  double maximum_rotary_diagonal = 0.0;
  std::array<double, 8> rotary_diagonals{};
  std::size_t rotary_index = 0U;
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const std::size_t offset = node * kPhysicalDofsPerNode;
    for (std::size_t rotation = 3U; rotation < 5U; ++rotation) {
      const double diagonal =
          physical_local(offset + rotation, offset + rotation);
      if (!std::isfinite(diagonal) || diagonal < 0.0) {
        return StiffnessFailure(
            source_location_, identity_,
            "MITC4 drilling stabilization requires finite nonnegative "
            "physical tangent-rotation diagonals.");
      }
      rotary_diagonals[rotary_index++] = diagonal;
      maximum_rotary_diagonal = (std::max)(maximum_rotary_diagonal, diagonal);
    }
  }
  if (!(maximum_rotary_diagonal > 0.0)) {
    return StiffnessFailure(
        source_location_, identity_,
        "MITC4 drilling stabilization requires a positive physical "
        "tangent-rotation diagonal inventory.");
  }
  if (!IsFinite(physical_local)) {
    return StiffnessFailure(
        source_location_, identity_,
        "MITC4 physical stiffness must contain only finite values.");
  }

  double normalized_diagonal_sum = 0.0;
  for (const double diagonal : rotary_diagonals) {
    normalized_diagonal_sum += diagonal / maximum_rotary_diagonal;
  }
  const double drilling_stiffness =
      1.0e-3 * maximum_rotary_diagonal * (normalized_diagonal_sum / 8.0);
  if (!std::isfinite(drilling_stiffness) || !(drilling_stiffness > 0.0)) {
    return StiffnessFailure(
        source_location_, identity_,
        "MITC4 drilling stiffness must be finite and positive.");
  }

  Matrix physical_global =
      Congruence(physical_local, PhysicalTransformation20());
  Matrix drilling_local{kNodeCount, kNodeCount};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    drilling_local(node, node) = drilling_stiffness;
  }
  Matrix drilling_global =
      Congruence(drilling_local, DrillingTransformation4());
  Matrix stabilized_global{kGlobalDofCount, kGlobalDofCount};
  for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
    for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
      stabilized_global(row, column) =
          physical_global(row, column) + drilling_global(row, column);
    }
  }
  if (!IsFinite(physical_global) || !IsFinite(drilling_global) ||
      !IsFinite(stabilized_global)) {
    return StiffnessFailure(
        source_location_, identity_,
        "MITC4 transformed stiffness must contain only finite values.");
  }

  return Result<Mitc4Stiffness>::Success(
      Mitc4Stiffness{std::move(physical_local), std::move(physical_global),
                     std::move(drilling_global), std::move(stabilized_global),
                     drilling_stiffness});
} catch (const std::exception& error) {
  return StiffnessFailure(source_location_, identity_, error.what());
}

Result<Mitc4Mass> Mitc4Shell::Mass(const double density) const {
  if (!std::isfinite(density) || !(density > 0.0)) {
    return MassFailure(source_location_, identity_,
                       "MITC4 mass requires finite density > 0.");
  }

  std::array<Matrix, 3> moments{Matrix{kNodeCount, kNodeCount},
                                Matrix{kNodeCount, kNodeCount},
                                Matrix{kNodeCount, kNodeCount}};
  constexpr double kGaussCoordinate = 0.77459666924148337704;
  constexpr std::array<double, 3> kGaussPoints{-kGaussCoordinate, 0.0,
                                               kGaussCoordinate};
  constexpr std::array<double, 3> kGaussWeights{5.0 / 9.0, 8.0 / 9.0,
                                                5.0 / 9.0};
  for (std::size_t xi_index = 0U; xi_index < kGaussPoints.size(); ++xi_index) {
    for (std::size_t eta_index = 0U; eta_index < kGaussPoints.size();
         ++eta_index) {
      const Mitc4ShapeFunctions shape =
          ShapeFunctions(kGaussPoints[xi_index], kGaussPoints[eta_index]);
      for (std::size_t zeta_index = 0U; zeta_index < kGaussPoints.size();
           ++zeta_index) {
        GeometryData geometry{};
        const double zeta = kGaussPoints[zeta_index];
        if (!EvaluateGeometry(kGaussPoints[xi_index], kGaussPoints[eta_index],
                              zeta, geometry)) {
          return MassFailure(source_location_, identity_,
                             "MITC4 mass quadrature geometry is invalid.");
        }
        const double z = 0.5 * thickness_ * zeta;
        const double weight =
            density * geometry.jacobian * kGaussWeights[xi_index] *
            kGaussWeights[eta_index] * kGaussWeights[zeta_index];
        if (!std::isfinite(z) || !std::isfinite(weight)) {
          return MassFailure(
              source_location_, identity_,
              "MITC4 mass quadrature contribution is nonfinite.");
        }
        const std::array<double, 3> z_powers{1.0, z, z * z};
        for (std::size_t row = 0U; row < kNodeCount; ++row) {
          for (std::size_t column = row; column < kNodeCount; ++column) {
            const double shape_product =
                shape.values[row] * shape.values[column];
            for (std::size_t moment = 0U; moment < moments.size(); ++moment) {
              const double contribution =
                  weight * shape_product * z_powers[moment];
              moments[moment](row, column) += contribution;
              if (row != column) {
                moments[moment](column, row) += contribution;
              }
            }
          }
        }
      }
    }
  }
  for (const Matrix& moment : moments) {
    if (!IsFinite(moment)) {
      return MassFailure(source_location_, identity_,
                         "MITC4 consistent moments must be finite.");
    }
  }

  Matrix consistent{kGlobalDofCount, kGlobalDofCount};
  for (std::size_t row_node = 0U; row_node < kNodeCount; ++row_node) {
    const Matrix row_skew = DirectorSkew(directors_[row_node]);
    for (std::size_t column_node = 0U; column_node < kNodeCount;
         ++column_node) {
      const Matrix column_skew = DirectorSkew(directors_[column_node]);
      for (std::size_t row = 0U; row < 3U; ++row) {
        for (std::size_t column = 0U; column < 3U; ++column) {
          const std::size_t global_row = kGlobalDofsPerNode * row_node + row;
          const std::size_t global_column =
              kGlobalDofsPerNode * column_node + column;
          consistent(global_row, global_column) =
              moments[0U](row_node, column_node) * (row == column ? 1.0 : 0.0);
          consistent(global_row, global_column + 3U) =
              moments[1U](row_node, column_node) * column_skew(row, column);
          consistent(global_row + 3U, global_column) =
              moments[1U](row_node, column_node) * row_skew(column, row);
          double rotary = 0.0;
          for (std::size_t axis = 0U; axis < 3U; ++axis) {
            rotary += row_skew(axis, row) * column_skew(axis, column);
          }
          consistent(global_row + 3U, global_column + 3U) =
              moments[2U](row_node, column_node) * rotary;
        }
      }
    }
  }

  Matrix lumped{kGlobalDofCount, kGlobalDofCount};
  std::array<std::array<double, 3>, kNodeCount> nodal_moments{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    double mass = 0.0;
    double first_moment = 0.0;
    double second_moment = 0.0;
    for (std::size_t column = 0U; column < kNodeCount; ++column) {
      mass += moments[0U](node, column);
      first_moment += moments[1U](node, column);
      second_moment += moments[2U](node, column);
    }
    const double centered_inertia =
        second_moment - first_moment * (first_moment / mass);
    if (!std::isfinite(mass) || !std::isfinite(first_moment) ||
        !std::isfinite(second_moment) || !std::isfinite(centered_inertia) ||
        !(mass > 0.0) || !(second_moment > 0.0) || !(centered_inertia > 0.0)) {
      return MassFailure(
          source_location_, identity_,
          "MITC4 nodal mass moments must be finite and positive.");
    }
    nodal_moments[node] = {mass, first_moment, second_moment};
    const Matrix skew = DirectorSkew(directors_[node]);
    const std::size_t offset = kGlobalDofsPerNode * node;
    for (std::size_t row = 0U; row < 3U; ++row) {
      for (std::size_t column = 0U; column < 3U; ++column) {
        lumped(offset + row, offset + column) =
            mass * (row == column ? 1.0 : 0.0);
        lumped(offset + row, offset + column + 3U) =
            first_moment * skew(row, column);
        lumped(offset + row + 3U, offset + column) =
            first_moment * skew(column, row);
        double projector = row == column ? 1.0 : 0.0;
        projector -= directors_[node][row] * directors_[node][column];
        lumped(offset + row + 3U, offset + column + 3U) =
            second_moment * projector;
      }
    }
  }
  if (!IsFinite(consistent) || !IsFinite(lumped)) {
    return MassFailure(
        source_location_, identity_,
        "MITC4 transformed mass must contain only finite values.");
  }
  return Result<Mitc4Mass>::Success(
      Mitc4Mass{std::move(moments), std::move(consistent), std::move(lumped),
                std::move(nodal_moments)});
}

Result<Mitc4PhysicalRecovery> Mitc4Shell::RecoverPhysical(
    const Vector& global_element_displacement24) const {
  if (global_element_displacement24.Size() != kGlobalDofCount) {
    return RecoveryFailure(
        source_location_, identity_,
        "MITC4 physical recovery requires exactly 24 global element DOFs.");
  }
  if (!IsFinite(global_element_displacement24)) {
    return RecoveryFailure(
        source_location_, identity_,
        "MITC4 physical recovery displacement must be finite.");
  }

  const Vector physical_displacement =
      PhysicalTransformation20().Multiply(global_element_displacement24);
  const Matrix tying_samples = CovariantTyingShearSamples20();
  const Matrix constitutive = MaterialConstitutive5();
  const Matrix plane_stress = PlaneStressConstitutive();
  const double gauss = 1.0 / std::sqrt(3.0);
  const std::array<std::array<double, 2>, 4> surface_points{
      std::array<double, 2>{-gauss, -gauss},
      std::array<double, 2>{gauss, -gauss}, std::array<double, 2>{gauss, gauss},
      std::array<double, 2>{-gauss, gauss}};
  constexpr std::array<double, 2> thickness_points{-1.0, 1.0};
  constexpr std::array<double, 3> section_positions{-1.0, 0.0, 1.0};

  Mitc4PhysicalRecovery recovery{};
  for (std::size_t surface = 0U; surface < surface_points.size(); ++surface) {
    auto& point = recovery.points[surface];
    point.natural_coordinates = surface_points[surface];
    GeometryData midsurface_geometry{};
    if (!EvaluateGeometry(point.natural_coordinates[0],
                          point.natural_coordinates[1], 0.0,
                          midsurface_geometry)) {
      return RecoveryFailure(source_location_, identity_,
                             "MITC4 midsurface recovery geometry is invalid.");
    }
    point.local_frame = midsurface_geometry.frame;

    for (double thickness_sign : thickness_points) {
      const double zeta = thickness_sign * gauss;
      const Matrix strain_matrix = StrainDisplacement(
          point.natural_coordinates[0], point.natural_coordinates[1], zeta,
          &tying_samples);
      const Vector strain = strain_matrix.Multiply(physical_displacement);
      const Vector stress = constitutive.Multiply(strain);
      GeometryData geometry{};
      if (!EvaluateGeometry(point.natural_coordinates[0],
                            point.natural_coordinates[1], zeta, geometry)) {
        return RecoveryFailure(
            source_location_, identity_,
            "Validated MITC4 recovery geometry became invalid.");
      }

      for (std::size_t component = 0U; component < 3U; ++component) {
        point.generalized_strain[component] += 0.5 * strain[component];
        point.generalized_strain[3U + component] +=
            3.0 * zeta * strain[component] / thickness_;
        point.section_resultant[component] +=
            0.5 * thickness_ * stress[component];
        point.section_resultant[3U + component] +=
            0.25 * thickness_ * thickness_ * zeta * stress[component];
      }
      for (std::size_t component = 0U; component < 2U; ++component) {
        point.generalized_strain[6U + component] +=
            0.5 * strain[3U + component];
        point.section_resultant[6U + component] +=
            0.5 * thickness_ * stress[3U + component];
      }
      recovery.strain_energy += 0.5 * strain.Dot(stress) * geometry.jacobian;
    }

    for (std::size_t position = 0U; position < section_positions.size();
         ++position) {
      GeometryData section_geometry{};
      if (!EvaluateGeometry(point.natural_coordinates[0],
                            point.natural_coordinates[1],
                            section_positions[position], section_geometry)) {
        return RecoveryFailure(
            source_location_, identity_,
            "MITC4 section-position recovery geometry is invalid.");
      }
      const Vector strain =
          StrainDisplacement(point.natural_coordinates[0],
                             point.natural_coordinates[1],
                             section_positions[position], &tying_samples)
              .Multiply(physical_displacement);
      Vector in_plane_strain{3U};
      for (std::size_t component = 0U; component < 3U; ++component) {
        in_plane_strain[component] = strain[component];
      }
      const Vector stress = plane_stress.Multiply(in_plane_strain);
      for (std::size_t component = 0U; component < 3U; ++component) {
        point.in_plane_stress[position][component] = stress[component];
      }
    }
  }

  if (!std::isfinite(recovery.strain_energy)) {
    return RecoveryFailure(source_location_, identity_,
                           "MITC4 physical strain energy must be finite.");
  }
  for (const auto& point : recovery.points) {
    const auto finite = [](const auto& values) {
      return std::all_of(values.begin(), values.end(),
                         [](double value) { return std::isfinite(value); });
    };
    if (!finite(point.generalized_strain) || !finite(point.section_resultant) ||
        std::any_of(
            point.in_plane_stress.begin(), point.in_plane_stress.end(),
            [&finite](const auto& stress) { return !finite(stress); })) {
      return RecoveryFailure(source_location_, identity_,
                             "MITC4 physical recovery values must be finite.");
    }
  }

  return Result<Mitc4PhysicalRecovery>::Success(std::move(recovery));
}

const ElementDofLayout& Mitc4Shell::DofLayout() const noexcept {
  return dof_layout_;
}

Result<ElementStiffnessContribution> Mitc4Shell::ComputeStiffness() const {
  if (dof_layout_.node_indices.size() != kNodeCount ||
      dof_layout_.components_per_node.size() != kGlobalDofsPerNode) {
    return RuntimeFailure<ElementStiffnessContribution>(
        runtime_location_, dof_layout_.source_id,
        "MITC4 runtime stiffness requires a bound four-node six-DOF layout.");
  }
  try {
    auto stiffness = Stiffness();
    if (!stiffness.HasValue()) {
      return Result<ElementStiffnessContribution>::Failure(
          stiffness.GetStatus());
    }
    return Result<ElementStiffnessContribution>::Success(
        {dof_layout_, std::move(stiffness.Value().stabilized_global24)});
  } catch (const std::exception& error) {
    return RuntimeFailure<ElementStiffnessContribution>(
        runtime_location_, dof_layout_.source_id, error.what());
  }
}

Result<ElementResultBundle> Mitc4Shell::Recover(
    const Vector& element_displacement) const {
  if (dof_layout_.node_indices.size() != kNodeCount ||
      dof_layout_.components_per_node.size() != kGlobalDofsPerNode ||
      element_displacement.Size() != kGlobalDofCount) {
    return RuntimeFailure<ElementResultBundle>(
        runtime_location_, dof_layout_.source_id,
        "MITC4 runtime recovery requires a bound four-node layout and 24 "
        "DOFs.");
  }

  Result<Mitc4PhysicalRecovery> recovered = [&]() {
    try {
      return RecoverPhysical(element_displacement);
    } catch (const std::exception& error) {
      return RuntimeFailure<Mitc4PhysicalRecovery>(
          runtime_location_, dof_layout_.source_id, error.what());
    }
  }();
  if (!recovered.HasValue()) {
    return Result<ElementResultBundle>::Failure(recovered.GetStatus());
  }

  constexpr std::array<ShellMidsurfaceLocation, kNodeCount> locations{
      ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
      ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
  constexpr std::array<ShellSectionPosition, 3> positions{
      ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
      ShellSectionPosition::kTop};
  constexpr std::array<double, 3> zeta{-1.0, 0.0, 1.0};

  ShellElementResultRows rows{};
  rows.rows.reserve(kNodeCount);
  rows.physical_strain_energy = recovered.Value().strain_energy;
  for (std::size_t point = 0U; point < recovered.Value().points.size();
       ++point) {
    const auto& physical_point = recovered.Value().points[point];
    ShellResultRow row{};
    row.element = element_index_;
    row.location = locations[point];
    row.natural_coordinates = physical_point.natural_coordinates;
    row.local_frame = {physical_point.local_frame.e1,
                       physical_point.local_frame.e2,
                       physical_point.local_frame.e3};
    row.generalized_strain = physical_point.generalized_strain;
    row.section_resultant = physical_point.section_resultant;
    for (std::size_t position = 0U; position < positions.size(); ++position) {
      row.stress[position] = {positions[position], zeta[position],
                              physical_point.in_plane_stress[position]};
    }
    rows.rows.push_back(std::move(row));
  }
  return Result<ElementResultBundle>::Success(
      {dof_layout_.source_id, std::move(rows)});
}

void Mitc4Shell::BindRuntime(ElementDofLayout layout,
                             const EntityIndex element_index,
                             SourceLocation location) {
  dof_layout_ = std::move(layout);
  element_index_ = element_index;
  runtime_location_ = std::move(location);
}

Mitc4Shell::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)
    : coordinates_{std::move(coordinates)},
      directors_{std::move(directors)},
      tangent_a_{std::move(tangent_a)},
      tangent_b_{std::move(tangent_b)},
      normal_candidate_{std::move(normal_candidate)},
      thickness_{thickness},
      youngs_modulus_{youngs_modulus},
      poisson_ratio_{poisson_ratio},
      source_location_{std::move(source_location)},
      identity_{std::move(identity)} {}

bool Mitc4Shell::EvaluateGeometry(double xi, double eta, double zeta,
                                  GeometryData& result) const noexcept {
  const auto shape = ShapeFunctions(xi, eta);
  const Vector3 midsurface_xi =
      DerivativeSum(shape.xi_derivatives, coordinates_);
  const Vector3 midsurface_eta =
      DerivativeSum(shape.eta_derivatives, coordinates_);
  const Vector3 director_xi = DerivativeSum(shape.xi_derivatives, directors_);
  const Vector3 director_eta = DerivativeSum(shape.eta_derivatives, directors_);
  const Vector3 director_value = WeightedSum(shape.values, directors_);
  const double half_thickness = 0.5 * thickness_;

  result.covariant[0] = midsurface_xi + half_thickness * zeta * director_xi;
  result.covariant[1] = midsurface_eta + half_thickness * zeta * director_eta;
  result.covariant[2] = half_thickness * director_value;
  result.jacobian =
      result.covariant[0].Dot(result.covariant[1].Cross(result.covariant[2]));
  if (!result.covariant[0].IsFinite() || !result.covariant[1].IsFinite() ||
      !result.covariant[2].IsFinite() || !std::isfinite(result.jacobian) ||
      !(result.jacobian > 0.0)) {
    return false;
  }

  if (!std::isfinite(1.0 / result.jacobian) ||
      !result.covariant[1].Cross(result.covariant[2]).IsFinite() ||
      !result.covariant[2].Cross(result.covariant[0]).IsFinite() ||
      !result.covariant[0].Cross(result.covariant[1]).IsFinite()) return false;
  // Preserve the reciprocal finiteness gate, streaming each action directly
  // into validation rather than retaining a reciprocal-basis matrix.
  try {
    const auto factor = dense_blas_internal::DenseFactor::Lu(
        CovariantMatrix(result.covariant));
    if (!factor.HasValue()) return false;
    Matrix rhs{3U, 1U};
    for (std::size_t axis = 0; axis < 3U; ++axis) {
      if (axis != 0U) rhs(axis - 1U, 0U) = 0.0;
      rhs(axis, 0U) = 1.0;
      if (!factor.Value().Solve(rhs, true).HasValue()) return false;
    }
  } catch (const std::exception&) {
    return false;
  }

  const Vector3 area = midsurface_xi.Cross(midsurface_eta);
  const double director_norm = director_value.Norm();
  if (!area.IsFinite() || !director_value.IsFinite() ||
      !std::isfinite(director_norm) || !(director_norm > 0.0) ||
      !(area.Dot(normal_candidate_) > 0.0)) {
    return false;
  }
  const Vector3 e3 = (1.0 / director_norm) * director_value;
  result.frame.e3 = e3.Components();
  if (!(area.Dot(e3) > 0.0)) {
    return false;
  }
  const Vector3 e1_candidate = midsurface_xi - midsurface_xi.Dot(e3) * e3;
  const double e1_norm = e1_candidate.Norm();
  if (!e1_candidate.IsFinite() || !std::isfinite(e1_norm) || !(e1_norm > 0.0)) {
    return false;
  }
  const Vector3 e1 = (1.0 / e1_norm) * e1_candidate;
  const Vector3 e2 = e3.Cross(e1);
  result.frame.e1 = e1.Components();
  result.frame.e2 = e2.Components();
  return e2.IsFinite();
}

std::array<std::array<Vector3, 3>, 20> Mitc4Shell::BasisDerivatives(
    double xi, double eta, double zeta) const noexcept {
  const auto shape = ShapeFunctions(xi, eta);
  const double half_thickness = 0.5 * thickness_;
  std::array<std::array<Vector3, 3>, kPhysicalDofCount> derivatives{};
  constexpr std::array<Vector3, 3> global_axes{
      Vector3{1.0, 0.0, 0.0}, Vector3{0.0, 1.0, 0.0}, Vector3{0.0, 0.0, 1.0}};

  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const std::size_t offset = node * kPhysicalDofsPerNode;
    for (std::size_t component = 0U; component < 3U; ++component) {
      derivatives[offset + component][0] =
          shape.xi_derivatives[node] * global_axes[component];
      derivatives[offset + component][1] =
          shape.eta_derivatives[node] * global_axes[component];
    }

    const Vector3 alpha_direction = -half_thickness * tangent_b_[node];
    derivatives[offset + 3U][0] =
        zeta * shape.xi_derivatives[node] * alpha_direction;
    derivatives[offset + 3U][1] =
        zeta * shape.eta_derivatives[node] * alpha_direction;
    derivatives[offset + 3U][2] = shape.values[node] * alpha_direction;

    const Vector3 beta_direction = half_thickness * tangent_a_[node];
    derivatives[offset + 4U][0] =
        zeta * shape.xi_derivatives[node] * beta_direction;
    derivatives[offset + 4U][1] =
        zeta * shape.eta_derivatives[node] * beta_direction;
    derivatives[offset + 4U][2] = shape.values[node] * beta_direction;
  }
  return derivatives;
}

Matrix Mitc4Shell::StrainDisplacement(double xi, double eta, double zeta,
                                      const Matrix* tying_samples) const {
  GeometryData geometry{};
  if (!EvaluateGeometry(xi, eta, zeta, geometry)) {
    throw std::invalid_argument{
        "MITC4 strain geometry is invalid at the requested point."};
  }
  const auto factor = dense_blas_internal::DenseFactor::Lu(
      CovariantMatrix(geometry.covariant));
  if (!factor.HasValue())
    throw std::invalid_argument{
        "MITC4 strain geometry is invalid at the requested point."};
  const auto derivatives = BasisDerivatives(xi, eta, zeta);
  const Mitc4TyingWeights weights = TyingWeights(xi, eta);
  Matrix result{5U, kPhysicalDofCount};
  for (std::size_t dof = 0U; dof < kPhysicalDofCount; ++dof) {
    auto covariant =
        CovariantStrainColumn(geometry.covariant, derivatives[dof]);
    if (tying_samples != nullptr) {
      covariant[0U][2U] = weights.xi_zeta[0] * (*tying_samples)(0U, dof) +
                          weights.xi_zeta[1] * (*tying_samples)(1U, dof);
      covariant[2U][0U] = covariant[0U][2U];
      covariant[1U][2U] = weights.eta_zeta[0] * (*tying_samples)(2U, dof) +
                          weights.eta_zeta[1] * (*tying_samples)(3U, dof);
      covariant[2U][1U] = covariant[1U][2U];
    }
    const auto engineering = LocalEngineeringComponents(
        ReconstructCartesianStrain(covariant, factor.Value()),
        geometry.frame);
    for (std::size_t component = 0U; component < engineering.size();
         ++component) {
      result(component, dof) = engineering[component];
    }
  }
  return result;
}

}  // namespace fesa