#include "fesa/elements/mitc4_shell.h"

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

namespace fesa {
namespace {

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;

Vector3 WeightedSum(const std::array<double, kNodeCount>& weights,
                    const std::array<Vector3, kNodeCount>& values) {
  Vector3 result{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    result = result + weights[node] * values[node];
  }
  return result;
}

Vector3 DerivativeSum(const std::array<double, kNodeCount>& derivatives,
                      const std::array<Vector3, kNodeCount>& values) {
  std::array<Vector3, kNodeCount> relative{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    relative[node] = values[node] - values[0];
  }
  return WeightedSum(derivatives, relative);
}

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;
}

std::string ElementIdentity(const std::array<const Node*, kNodeCount>& nodes) {
  std::string identity;
  for (const Node* node : nodes) {
    if (node == nullptr) {
      continue;
    }
    if (!identity.empty()) {
      identity += "-";
    }
    identity += node->source_id.source_label_text;
  }
  return identity;
}

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 Reconstructs Cartesian strain in the approved reciprocal-basis
/// reduction order.
std::array<std::array<double, 3>, 3> ReconstructCartesianStrain(
    const std::array<std::array<double, 3>, 3>& covariant_strain,
    const std::array<Vector3, 3>& reciprocal) {
  std::array<std::array<double, 3>, 3> tensor{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      for (std::size_t first = 0U; first < 3U; ++first) {
        for (std::size_t second = 0U; second < 3U; ++second) {
          tensor[row][column] += covariant_strain[first][second] *
                                 reciprocal[first][row] *
                                 reciprocal[second][column];
        }
      }
    }
  }
  return tensor;
}

double FrameComponent(const Vector3& left,
                      const std::array<std::array<double, 3>, 3>& tensor,
                      const Vector3& right) {
  double value = 0.0;
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      value += left[row] * tensor[row][column] * right[column];
    }
  }
  return value;
}

/// @brief Projects tensor strain into signed local engineering components.
std::array<double, 5> LocalEngineeringComponents(
    const std::array<std::array<double, 3>, 3>& tensor,
    const Mitc4LocalFrame& frame) {
  const Vector3 e1{frame.e1};
  const Vector3 e2{frame.e2};
  const Vector3 e3{frame.e3};
  return {FrameComponent(e1, tensor, e1), FrameComponent(e2, tensor, e2),
          2.0 * FrameComponent(e1, tensor, e2),
          2.0 * FrameComponent(e1, tensor, e3),
          2.0 * FrameComponent(e2, tensor, e3)};
}

Matrix ScaledMatrix(const Matrix& source, double factor) {
  Matrix result{source.Rows(), source.Columns()};
  for (std::size_t row = 0U; row < source.Rows(); ++row) {
    for (std::size_t column = 0U; column < source.Columns(); ++column) {
      result(row, column) = factor * source(row, column);
    }
  }
  return result;
}

/// @brief Lifts a local matrix by a deterministic serial congruence product.
/// @note The loop and accumulation order are part of reproducible stiffness.
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{transformation.Columns(), transformation.Columns()};
  for (std::size_t row = 0U; row < result.Rows(); ++row) {
    for (std::size_t column = row; column < result.Columns(); ++column) {
      double value = 0.0;
      for (std::size_t local_row = 0U; local_row < local.Rows(); ++local_row) {
        for (std::size_t local_column = 0U; local_column < local.Columns();
             ++local_column) {
          value += transformation(local_row, row) *
                   local(local_row, local_column) *
                   transformation(local_column, column);
        }
      }
      result(row, column) = value;
      result(column, row) = value;
    }
  }
  return result;
}

bool IsFinite(const Matrix& matrix) {
  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
      if (!std::isfinite(matrix(row, column))) {
        return false;
      }
    }
  }
  return true;
}

bool IsFinite(const Vector& vector) {
  for (std::size_t index = 0U; index < vector.Size(); ++index) {
    if (!std::isfinite(vector[index])) {
      return false;
    }
  }
  return true;
}

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<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)}}));
}

}  // 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 {
  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);
    for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) {
      for (std::size_t column = row; column < kPhysicalDofCount; ++column) {
        double integrand = 0.0;
        for (std::size_t first = 0U; first < 5U; ++first) {
          for (std::size_t second = 0U; second < 5U; ++second) {
            integrand += strain(first, row) * constitutive(first, second) *
                         strain(second, column);
          }
        }
        const double contribution =
            integrand * geometry.jacobian * point.weight;
        physical_local(row, column) += contribution;
        if (row != column) {
          physical_local(column, row) += contribution;
        }
      }
    }
  }

  double drilling_reference = (std::numeric_limits<double>::max)();
  bool has_drilling_reference = false;
  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) {
        drilling_reference = (std::min)(drilling_reference, diagonal);
        has_drilling_reference = true;
      }
    }
  }
  if (!has_drilling_reference) {
    return StiffnessFailure(
        source_location_, identity_,
        "MITC4 drilling stabilization requires a finite positive physical "
        "tangent-rotation diagonal.");
  }
  if (!IsFinite(physical_local)) {
    return StiffnessFailure(
        source_location_, identity_,
        "MITC4 physical stiffness must contain only finite values.");
  }

  const double drilling_stiffness = 1.0e-3 * drilling_reference;
  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});
}

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));
}

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;
  }

  result.reciprocal[0] =
      (1.0 / result.jacobian) * result.covariant[1].Cross(result.covariant[2]);
  result.reciprocal[1] =
      (1.0 / result.jacobian) * result.covariant[2].Cross(result.covariant[0]);
  result.reciprocal[2] =
      (1.0 / result.jacobian) * result.covariant[0].Cross(result.covariant[1]);

  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 result.reciprocal[0].IsFinite() && result.reciprocal[1].IsFinite() &&
         result.reciprocal[2].IsFinite() && 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 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, geometry.reciprocal),
        geometry.frame);
    for (std::size_t component = 0U; component < engineering.size();
         ++component) {
      result(component, dof) = engineering[component];
    }
  }
  return result;
}

}  // namespace fesa