#include "fesa/elements/nonlinear_mitc4_shell.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <memory>
#include <optional>
#include <string>
#include <utility>

#include "fesa/elements/finite_rotation.h"
#include "elements/nonlinear_mitc4_configuration_internal.h"
#include "math/dense_blas_internal.h"
#include "math/jet_arithmetic.h"
#include "math/jet_linear_solve.h"
#include "math/finite_values.h"
#include "math/vector3_operations.h"
#include "elements/shell_kernel_values.h"

namespace fesa {

nonlinear_mitc4_internal::ConfigurationSnapshot
nonlinear_mitc4_internal::ConfigurationSnapshot::Read(
    const NonlinearMitc4Shell& element) {
  ConfigurationSnapshot result{};
  std::size_t index = 0U;
  for (const auto& coordinates : element.coordinates_)
    for (std::size_t c = 0U; c < 3U; ++c)
      result.values[index++] = coordinates[c];
  for (const auto& frame : element.initial_frames_)
    for (std::size_t i = 0U; i < 3U; ++i)
      for (std::size_t j = 0U; j < 3U; ++j)
        result.values[index++] = frame(i, j);
  for (const auto& director : element.initial_directors_)
    for (std::size_t c = 0U; c < 3U; ++c) result.values[index++] = director[c];
  result.values[index++] = element.thickness_;
  result.values[index++] = element.youngs_modulus_;
  result.values[index++] = element.poisson_ratio_;
  result.values[index] = element.drilling_stiffness_;
  result.source = element.source_location_;
  result.identity = element.identity_;
  return result;
}

namespace {

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

constexpr std::size_t kNodeCount = 4U;
constexpr double kFrameTolerance = 1.0e-12;
constexpr double kDirectorOppositionTolerance = 1.0e-8;
constexpr double kShearCorrection = 5.0 / 6.0;

using Tensor = std::array<std::array<double, 3>, 3>;
using EngineeringStrain = std::array<double, 5>;
using EngineeringStress = std::array<double, 5>;

struct Geometry {
  std::array<Vector3, 3> reference_covariant;
  std::array<Vector3, 3> current_covariant;
  std::array<Vector3, 3> covariant_increment;
  std::optional<dense_blas_internal::DenseFactor> reference_factor;
  Mitc4LocalFrame local_frame;
  double reference_jacobian;
};

struct PointEvaluation {
  Geometry geometry;
  EngineeringStrain strain;
  EngineeringStress stress;
};

/// @brief Checks every scalar in a tensor before recovery commits it.
bool IsFinite(const Tensor& tensor) {
  for (const auto& row : tensor) {
    for (const double value : row) {
      if (!std::isfinite(value)) {
        return false;
      }
    }
  }
  return true;
}

/// @brief Evaluates a determinant from covariant basis columns.
double Determinant(const std::array<Vector3, 3>& covariant) {
  return covariant[0U].Dot(covariant[1U].Cross(covariant[2U]));
}

/// @brief Builds one proper fixed reference local frame at a midsurface point.
bool MakeLocalFrame(const Vector3& midsurface_xi, const Vector3& midsurface_eta,
                    const Vector3& director, Mitc4LocalFrame& frame) {
  const Vector3 area = midsurface_xi.Cross(midsurface_eta);
  const double director_norm = director.Norm();
  if (!area.IsFinite() || !director.IsFinite() ||
      !(director_norm > kDirectorOppositionTolerance) ||
      !std::isfinite(director_norm)) {
    return false;
  }

  const Vector3 e3 = director / director_norm;
  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 = e1_candidate / e1_norm;
  const Vector3 e2 = e3.Cross(e1);
  if (!e2.IsFinite()) {
    return false;
  }
  frame = {e1.Components(), e2.Components(), e3.Components()};
  return true;
}

/// @brief Forms (R(phi)-I)D0 without adding and subtracting the baseline.
Vector3 DirectorIncrement(const Vector3& phi, const Vector3& initial) {
  const double squared_norm = phi.Dot(phi);
  double first = 1.0;
  double second = 0.5;
  if (squared_norm < 0.25) {
    double power = 1.0;
    double first_coefficient = 1.0;
    double second_coefficient = 0.5;
    for (unsigned int k = 1U; k <= 9U; ++k) {
      power *= squared_norm;
      first_coefficient /= -static_cast<double>((2U * k) * (2U * k + 1U));
      second_coefficient /= -static_cast<double>((2U * k + 1U) * (2U * k + 2U));
      first += first_coefficient * power;
      second += second_coefficient * power;
    }
  } else {
    const double norm = std::sqrt(squared_norm);
    first = std::sin(norm) / norm;
    second = (1.0 - std::cos(norm)) / squared_norm;
  }
  const Vector3 crossed = phi.Cross(initial);
  return first * crossed + second * phi.Cross(crossed);
}

/// @brief Evaluates fixed-reference and trial covariant geometry together.
bool EvaluateGeometry(const std::array<Vector3, kNodeCount>& coordinates,
                      const std::array<Vector3, kNodeCount>& initial_directors,
                      const std::array<Vector3, kNodeCount>& directors,
                      const NonlinearMitc4TrialDofs& trial,
                      const double thickness, const double xi, const double eta,
                      const double zeta, Geometry& geometry) {
  const Mitc4ShapeFunctions shape = Mitc4Shell::ShapeFunctions(xi, eta);
  const Vector3 reference_xi = DerivativeSum(shape.xi_derivatives, coordinates);
  const Vector3 reference_eta =
      DerivativeSum(shape.eta_derivatives, coordinates);
  const Vector3 initial_director_xi =
      DerivativeSum(shape.xi_derivatives, initial_directors);
  const Vector3 initial_director_eta =
      DerivativeSum(shape.eta_derivatives, initial_directors);
  const Vector3 initial_director = WeightedSum(shape.values, initial_directors);
  const Vector3 displacement_xi =
      DerivativeSum(shape.xi_derivatives, trial.displacement);
  const Vector3 displacement_eta =
      DerivativeSum(shape.eta_derivatives, trial.displacement);
  const Vector3 director_xi = DerivativeSum(shape.xi_derivatives, directors);
  const Vector3 director_eta = DerivativeSum(shape.eta_derivatives, directors);
  const Vector3 director = WeightedSum(shape.values, directors);
  const double half_thickness = 0.5 * thickness;
  std::array<Vector3, kNodeCount> director_increment{};
  for (std::size_t node = 0U; node < kNodeCount; ++node)
    director_increment[node] =
        DirectorIncrement(trial.rotation_vector[node], initial_directors[node]);
  geometry.covariant_increment[0U] =
      displacement_xi + half_thickness * zeta *
                            DerivativeSum(shape.xi_derivatives, director_increment);
  geometry.covariant_increment[1U] =
      displacement_eta + half_thickness * zeta *
                             DerivativeSum(shape.eta_derivatives, director_increment);
  geometry.covariant_increment[2U] =
      half_thickness * WeightedSum(shape.values, director_increment);

  geometry.reference_covariant[0U] =
      reference_xi + half_thickness * zeta * initial_director_xi;
  geometry.reference_covariant[1U] =
      reference_eta + half_thickness * zeta * initial_director_eta;
  geometry.reference_covariant[2U] = half_thickness * initial_director;
  geometry.current_covariant[0U] =
      reference_xi + displacement_xi + half_thickness * zeta * director_xi;
  geometry.current_covariant[1U] =
      reference_eta + displacement_eta + half_thickness * zeta * director_eta;
  geometry.current_covariant[2U] = half_thickness * director;
  geometry.reference_jacobian = Determinant(geometry.reference_covariant);
  const double current_jacobian = Determinant(geometry.current_covariant);

  if (!reference_xi.IsFinite() || !reference_eta.IsFinite() ||
      !initial_director.IsFinite() || !director.IsFinite() ||
      !(director.Norm() > kDirectorOppositionTolerance) ||
      !std::isfinite(geometry.reference_jacobian) ||
      !(geometry.reference_jacobian > 0.0) ||
      !std::isfinite(current_jacobian) || !(current_jacobian > 0.0)) {
    return false;
  }

  // Retain the old unscaled cross-product finite gate before factor actions.
  for (std::size_t alpha = 0; alpha < 3; ++alpha)
    if (!geometry.reference_covariant[(alpha+1)%3].Cross(
            geometry.reference_covariant[(alpha+2)%3]).IsFinite()) return false;
  Matrix basis{3, 3};
  for (std::size_t row = 0; row < 3; ++row)
    for (std::size_t alpha = 0; alpha < 3; ++alpha)
      basis(row, alpha) = geometry.reference_covariant[alpha][row];
  auto factor = dense_blas_internal::DenseFactor::Lu(basis);
  if (!factor.HasValue()) return false;
  // Stream the former reciprocal finite gate without retaining an inverse.
  Matrix rhs{3, 1};
  for (std::size_t column = 0; column < 3; ++column) {
    if (column != 0) rhs(column-1, 0) = 0;
    rhs(column, 0) = 1;
    if (!factor.Value().Solve(rhs).HasValue()) return false;
  }
  geometry.reference_factor = std::move(factor.Value());
  return MakeLocalFrame(reference_xi, reference_eta, initial_director,
                        geometry.local_frame);
}

/// @brief Forms the direct Green-Lagrange covariant strain with E33 excluded.
Tensor DirectStrain(const Geometry& geometry) {
  Tensor strain{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      strain[row][column] =
          0.5 * (geometry.reference_covariant[row].Dot(
                     geometry.covariant_increment[column]) +
                 geometry.covariant_increment[row].Dot(
                     geometry.reference_covariant[column]) +
                 geometry.covariant_increment[row].Dot(
                     geometry.covariant_increment[column]));
    }
  }
  strain[2U][2U] = 0.0;
  return strain;
}

/// @brief Reconstructs a fixed-reference Cartesian strain from covariant data.
Tensor CartesianStrain(const Tensor& covariant, const Geometry& reference) {
  Matrix rhs{3, 3};
  for (std::size_t row = 0; row < 3; ++row)
    for (std::size_t column = 0; column < 3; ++column)
      rhs(row, column) = covariant[row][column];
  const auto right = reference.reference_factor->SolveRight(rhs);
  Tensor cartesian{};
  if (right.HasValue()) {
    const auto left = reference.reference_factor->Solve(right.Value(), true);
    if (left.HasValue()) {
      for (std::size_t row = 0; row < 3; ++row)
        for (std::size_t column = 0; column < 3; ++column)
          cartesian[row][column] = left.Value()(row, column);
      return cartesian;
    }
  }
  for (auto& row : cartesian)
    row.fill(std::numeric_limits<double>::quiet_NaN());
  return cartesian;
}

/// @brief Applies the approved homogeneous plane-stress plus shear law.
EngineeringStress Stress(const EngineeringStrain& strain, const double youngs,
                         const double poisson) {
  const double factor = youngs / (1.0 - poisson * poisson);
  const double shear = kShearCorrection * youngs / (2.0 * (1.0 + poisson));
  return {factor * (strain[0U] + poisson * strain[1U]),
          factor * (poisson * strain[0U] + strain[1U]),
          factor * (1.0 - poisson) * 0.5 * strain[2U], shear * strain[3U],
          shear * strain[4U]};
}

/// @brief Distinguishes invalid geometry from arithmetic recovery failure.
enum class PointStatus { kValid, kInvalidGeometry, kNonfiniteRecovery };

/// @brief Evaluates direct strain at one mandatory edge-midpoint tying point.
PointStatus EvaluateTyingSample(
    const std::array<Vector3, kNodeCount>& coordinates,
    const std::array<Vector3, kNodeCount>& initial_directors,
    const std::array<Vector3, kNodeCount>& directors,
    const NonlinearMitc4TrialDofs& trial, const double thickness,
    const double xi, const double eta, Tensor& strain) {
  Geometry geometry{};
  if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
                        thickness, xi, eta, 0.0, geometry)) {
    return PointStatus::kInvalidGeometry;
  }
  strain = DirectStrain(geometry);
  return IsFinite(strain) ? PointStatus::kValid : PointStatus::kNonfiniteRecovery;
}

/// @brief Computes one MITC-projected local strain and material stress point.
PointStatus EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
                   const std::array<Vector3, kNodeCount>& initial_directors,
                   const std::array<Vector3, kNodeCount>& directors,
                   const NonlinearMitc4TrialDofs& trial,
                   const std::array<Tensor, 4>& tying_samples,
                   const double thickness, const double youngs,
                   const double poisson, const double xi, const double eta,
                   const double zeta, PointEvaluation& point) {
  if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
                        thickness, xi, eta, zeta, point.geometry)) {
    return PointStatus::kInvalidGeometry;
  }
  Tensor covariant = DirectStrain(point.geometry);
  const Mitc4TyingWeights weights = Mitc4Shell::TyingWeights(xi, eta);
  covariant[0U][2U] = weights.xi_zeta[0U] * tying_samples[0U][0U][2U] +
                      weights.xi_zeta[1U] * tying_samples[1U][0U][2U];
  covariant[2U][0U] = covariant[0U][2U];
  covariant[1U][2U] = weights.eta_zeta[0U] * tying_samples[2U][1U][2U] +
                      weights.eta_zeta[1U] * tying_samples[3U][1U][2U];
  covariant[2U][1U] = covariant[1U][2U];
  const Tensor cartesian =
      CartesianStrain(covariant, point.geometry);
  point.strain = LocalEngineeringComponents(cartesian, point.geometry.local_frame);
  point.stress = Stress(point.strain, youngs, poisson);
  return IsFinite(covariant) && IsFinite(cartesian) && IsFinite(point.strain) &&
                 IsFinite(point.stress)
             ? PointStatus::kValid
             : PointStatus::kNonfiniteRecovery;
}

/// @brief Tests the fixed-order physical recovery inventory for finite values.
bool IsFinite(const NonlinearMitc4Recovery& recovery) {
  if (!std::isfinite(recovery.physical_strain_energy)) {
    return false;
  }
  for (const auto& point : recovery.gauss_points) {
    if (!IsFinite(point.generalized_strain) ||
        !IsFinite(point.section_resultant)) {
      return false;
    }
    for (const auto& stress : point.stress) {
      if (!IsFinite(stress.components)) {
        return false;
      }
    }
  }
  return true;
}

constexpr std::size_t kGlobalDofCount = 24U;
constexpr std::size_t kHessianSize = kGlobalDofCount * kGlobalDofCount;

/// @brief Carries exact first and second trial-coordinate variations.
class Hessian {
 public:
  Hessian() = default;

  double& operator[](const std::size_t index) { return Data()[index]; }

  const double& operator[](const std::size_t index) const {
    return Storage()[index];
  }

  auto begin() const { return Storage().begin(); }

  auto end() const { return Storage().end(); }

  double* Data() {
    if (!values_) {
      values_ = std::make_shared<Values>();
      values_->entries.fill(0.0);
    }
    return values_->entries.data();
  }

  const double* ReadData() const { return Storage().data(); }

  /// @brief Identifies only the immutable zero Hessian, independently of seeds.
  bool IsImplicitZero() const { return !values_; }

  /// @brief Gets fresh output storage without redundant zero initialization.
  /// @note Only fresh operator results use this pointer; all 576 entries must
  /// be written before any read or copy of the result.
  double* OverwriteData() {
    if (!values_) values_ = std::make_shared<Values>();
    return values_->entries.data();
  }

 private:
  struct Values {
    // User-provided construction leaves entries uninitialized for a complete
    // output overwrite while make_shared retains the single allocation.
    Values() noexcept {}
    std::array<double, kHessianSize> entries;
  };

  const std::array<double, kHessianSize>& Storage() const {
    // Constants borrow only this immutable zero. Mutable coefficients retain
    // the existing shared ownership; no trial-dependent cache is introduced.
    static const std::array<double, kHessianSize> kZero{};
    return values_ ? values_->entries : kZero;
  }

  std::shared_ptr<Values> values_;
};

/// @brief Forms an ordinary derivative outer product without a Matrix copy.
/// @note Rejected nonfinite operands propagate to the existing Jet finite gate.
void Outer(const std::array<double, kGlobalDofCount>& left,
           const std::array<double, kGlobalDofCount>& right, double* result) {
  if (!dense_blas_internal::OuterProduct(left.data(), left.size(), right.data(),
                                         right.size(), result))
    std::fill_n(result, kHessianSize, std::numeric_limits<double>::quiet_NaN());
}

/// @brief Carries exact first and second trial-coordinate variations.
struct Jet {
  double value{};
  std::array<double, kGlobalDofCount> first{};
  Hessian second;
  double directional_first{};
  double directional_second{};
  std::array<double, kGlobalDofCount> mixed_first{};
  std::array<double, kGlobalDofCount> mixed_second{};
  // Provenance for numerically zero derivatives; a null Hessian alone is not
  // constant provenance because Variable retains a nonzero gradient seed.
  bool constant_derivatives{false};
};

using JetVector3 = std::array<Jet, 3>;
using JetMatrix3 = std::array<std::array<Jet, 3>, 3>;
using JetTensor = std::array<std::array<Jet, 3>, 3>;

/// @brief Returns a constant second-order trial-coordinate quantity.
Jet Constant(const double value) {
  Jet result{};
  result.value = value;
  result.constant_derivatives = true;
  return result;
}

/// @brief Returns one independent global trial coordinate.
Jet Variable(const double value, const std::size_t index,
             const double direction = 0.0) {
  Jet result = Constant(value);
  result.constant_derivatives = false;
  result.first[index] = 1.0;
  result.directional_first = direction;
  return result;
}

/// @brief Adds two second-order trial-coordinate quantities.
Jet operator+(const Jet& left, const Jet& right) {
  Jet result{};
  result.value = left.value + right.value;
  result.directional_first = left.directional_first + right.directional_first;
  result.directional_second =
      left.directional_second + right.directional_second;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = left.first[index] + right.first[index];
    result.mixed_first[index] =
        left.mixed_first[index] + right.mixed_first[index];
    result.mixed_second[index] =
        left.mixed_second[index] + right.mixed_second[index];
  }
  result.constant_derivatives = left.constant_derivatives &&
                                right.constant_derivatives &&
                                std::isfinite(result.value);
  if (left.second.IsImplicitZero() && right.second.IsImplicitZero())
    return result;
  const double* left_second = left.second.ReadData();
  const double* right_second = right.second.ReadData();
  double* result_second = result.second.OverwriteData();
  for (std::size_t index = 0U; index < kHessianSize; ++index) {
    result_second[index] = left_second[index] + right_second[index];
  }
  return result;
}

/// @brief Subtracts two second-order trial-coordinate quantities.
Jet operator-(const Jet& left, const Jet& right) {
  Jet result{};
  result.value = left.value - right.value;
  result.directional_first = left.directional_first - right.directional_first;
  result.directional_second =
      left.directional_second - right.directional_second;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = left.first[index] - right.first[index];
    result.mixed_first[index] =
        left.mixed_first[index] - right.mixed_first[index];
    result.mixed_second[index] =
        left.mixed_second[index] - right.mixed_second[index];
  }
  result.constant_derivatives = left.constant_derivatives &&
                                right.constant_derivatives &&
                                std::isfinite(result.value);
  if (left.second.IsImplicitZero() && right.second.IsImplicitZero())
    return result;
  const double* left_second = left.second.ReadData();
  const double* right_second = right.second.ReadData();
  double* result_second = result.second.OverwriteData();
  for (std::size_t index = 0U; index < kHessianSize; ++index) {
    result_second[index] = left_second[index] - right_second[index];
  }
  return result;
}

/// @brief Multiplies two second-order trial-coordinate quantities.
Jet operator*(const Jet& left, const Jet& right) {
  Jet result{};
  result.value = left.value * right.value;
  result.directional_first = left.directional_first * right.value +
                             left.value * right.directional_first;
  result.directional_second =
      left.directional_second * right.value +
      2.0 * left.directional_first * right.directional_first +
      left.value * right.directional_second;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] =
        left.first[index] * right.value + left.value * right.first[index];
    result.mixed_first[index] = left.mixed_first[index] * right.value +
                                left.first[index] * right.directional_first +
                                left.directional_first * right.first[index] +
                                left.value * right.mixed_first[index];
    result.mixed_second[index] =
        left.mixed_second[index] * right.value +
        2.0 * left.mixed_first[index] * right.directional_first +
        left.first[index] * right.directional_second +
        left.directional_second * right.first[index] +
        2.0 * left.directional_first * right.mixed_first[index] +
        left.value * right.mixed_second[index];
  }
  const bool finite_values = std::isfinite(left.value) &&
                             std::isfinite(right.value) &&
                             std::isfinite(result.value);
  const auto finite_gradient = [](const auto& gradient) {
    return std::all_of(gradient.begin(), gradient.end(),
                       [](const double value) { return std::isfinite(value); });
  };
  const bool zero_outer = finite_values &&
      ((left.constant_derivatives && finite_gradient(right.first)) ||
       (right.constant_derivatives && finite_gradient(left.first)));
  result.constant_derivatives = finite_values && left.constant_derivatives &&
                                right.constant_derivatives;
  if (zero_outer && left.second.IsImplicitZero() &&
      right.second.IsImplicitZero()) return result;
  const double* left_second = left.second.ReadData();
  const double* right_second = right.second.ReadData();
  double* result_second = result.second.OverwriteData();
  if (zero_outer) {
    for (std::size_t index = 0U; index < kHessianSize; ++index)
      result_second[index] = left_second[index] * right.value + 0.0 + 0.0 +
                             left.value * right_second[index];
    return result;
  }
  Outer(left.first, right.first, result_second);
  std::array<double, kHessianSize> reverse;
  Outer(right.first, left.first, reverse.data());
  for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
    for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
      const std::size_t index = row * kGlobalDofCount + column;
      result_second[index] = left_second[index] * right.value +
                             result_second[index] + reverse[index] +
                             left.value * right_second[index];
    }
  }
  return result;
}

/// @brief Scales a second-order trial-coordinate quantity.
Jet operator*(const double scalar, const Jet& value) {
  Jet result{};
  result.value = scalar * value.value;
  result.directional_first = scalar * value.directional_first;
  result.directional_second = scalar * value.directional_second;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = scalar * value.first[index];
    result.mixed_first[index] = scalar * value.mixed_first[index];
    result.mixed_second[index] = scalar * value.mixed_second[index];
  }
  result.constant_derivatives = value.constant_derivatives &&
                                std::isfinite(scalar) &&
                                std::isfinite(result.value);
  if (std::isfinite(scalar) && value.second.IsImplicitZero()) return result;
  const double* value_second = value.second.ReadData();
  double* result_second = result.second.OverwriteData();
  for (std::size_t index = 0U; index < kHessianSize; ++index) {
    result_second[index] = scalar * value_second[index];
  }
  return result;
}

/// @brief Scales a second-order trial-coordinate quantity.
Jet operator*(const Jet& value, const double scalar) { return scalar * value; }

/// @brief Applies a scalar unary function through its first two derivatives.
Jet Compose(const Jet& input, const double value, const double first,
            const double second, const double third) {
  Jet result{};
  result.value = value;
  result.directional_first = first * input.directional_first;
  result.directional_second =
      second * input.directional_first * input.directional_first +
      first * input.directional_second;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = first * input.first[index];
    result.mixed_first[index] =
        second * input.first[index] * input.directional_first +
        first * input.mixed_first[index];
    result.mixed_second[index] =
        third * input.first[index] * input.directional_first *
            input.directional_first +
        second * (2.0 * input.mixed_first[index] * input.directional_first +
                  input.first[index] * input.directional_second) +
        first * input.mixed_second[index];
  }
  result.constant_derivatives = input.constant_derivatives &&
      std::isfinite(input.value) && std::isfinite(value) &&
      std::isfinite(first) && std::isfinite(second) && std::isfinite(third);
  if (result.constant_derivatives) return result;
  std::array<double, kGlobalDofCount> scaled_first{};
  for (std::size_t row = 0U; row < kGlobalDofCount; ++row)
    scaled_first[row] = second * input.first[row];
  const double* input_second = input.second.ReadData();
  double* result_second = result.second.OverwriteData();
  Outer(scaled_first, input.first, result_second);
  for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
    for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
      const std::size_t index = row * kGlobalDofCount + column;
      result_second[index] = first * input_second[index] + result_second[index];
    }
  }
  return result;
}

/// @brief Forms an analytic reciprocal without a secant approximation.
Jet Reciprocal(const Jet& value) {
  const double inverse = 1.0 / value.value;
  return Compose(value, inverse, -inverse * inverse,
                 2.0 * inverse * inverse * inverse,
                 -6.0 * inverse * inverse * inverse * inverse);
}

/// @brief Divides two second-order trial-coordinate quantities.
Jet operator/(const Jet& left, const Jet& right) {
  return left * Reciprocal(right);
}

/// @brief Computes the analytic square root variation.
Jet Sqrt(const Jet& value) {
  const double root = std::sqrt(value.value);
  return Compose(value, root, 0.5 / root, -0.25 / (root * root * root),
                 0.375 / (root * root * root * root * root));
}

/// @brief Computes the analytic sine variation.
Jet Sin(const Jet& value) {
  return Compose(value, std::sin(value.value), std::cos(value.value),
                 -std::sin(value.value), -std::cos(value.value));
}

/// @brief Computes the analytic cosine variation.
Jet Cos(const Jet& value) {
  return Compose(value, std::cos(value.value), -std::sin(value.value),
                 -std::cos(value.value), std::sin(value.value));
}

/// @brief Computes the two-argument angle and its analytic first/second map.
Jet Atan2(const Jet& y, const Jet& x) {
  const double radius_squared = x.value * x.value + y.value * y.value;
  Jet result{};
  result.value = std::atan2(y.value, x.value);
  const double derivative_y = x.value / radius_squared;
  const double derivative_x = -y.value / radius_squared;
  const double radius_fourth = radius_squared * radius_squared;
  const double derivative_yy = -2.0 * x.value * y.value / radius_fourth;
  const double derivative_yx =
      (y.value * y.value - x.value * x.value) / radius_fourth;
  const double derivative_xx = 2.0 * x.value * y.value / radius_fourth;
  const bool zero_hessian = x.constant_derivatives && y.constant_derivatives &&
      std::isfinite(result.value) && std::isfinite(derivative_y) &&
      std::isfinite(derivative_x) && std::isfinite(derivative_yy) &&
      std::isfinite(derivative_yx) && std::isfinite(derivative_xx);
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] =
        derivative_y * y.first[index] + derivative_x * x.first[index];
  }
  if (!zero_hessian) {
    std::array<double, kGlobalDofCount> scaled_y{};
    std::array<double, kGlobalDofCount> scaled_x{};
    for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
      scaled_y[i] = derivative_yy * y.first[i];
      scaled_x[i] = derivative_xx * x.first[i];
    }
    std::array<std::array<double, kHessianSize>, 4U> products;
    Outer(scaled_y, y.first, products[0U].data());
    Outer(y.first, x.first, products[1U].data());
    Outer(x.first, y.first, products[2U].data());
    Outer(scaled_x, x.first, products[3U].data());
    const double* y_second = y.second.ReadData();
    const double* x_second = x.second.ReadData();
    double* result_second = result.second.OverwriteData();
    for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
      for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
        const std::size_t index = row * kGlobalDofCount + column;
        result_second[index] =
            derivative_y * y_second[index] + derivative_x * x_second[index] +
            products[0U][index] +
            derivative_yx * (products[1U][index] + products[2U][index]) +
            products[3U][index];
      }
    }
  }
  // Use a bounded ratio for the directional third derivative only. Retain
  // the established value/gradient/Hessian reduction above for static runs.
  const bool use_x = std::abs(x.value) >= std::abs(y.value);
  const Jet ratio = use_x ? y / x : x / y;
  const double denominator = 1.0 + ratio.value * ratio.value;
  const Jet angle = Compose(ratio, std::atan(ratio.value), 1.0 / denominator,
                            -2.0 * ratio.value / (denominator * denominator),
                            (6.0 * ratio.value * ratio.value - 2.0) /
                                (denominator * denominator * denominator));
  const double sign = use_x ? 1.0 : -1.0;
  result.directional_first = sign * angle.directional_first;
  result.directional_second = sign * angle.directional_second;
  for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
    result.mixed_first[i] = sign * angle.mixed_first[i];
    result.mixed_second[i] = sign * angle.mixed_second[i];
  }
  result.constant_derivatives = zero_hessian && angle.constant_derivatives;
  return result;
}

/// @brief Checks every coefficient in an automatic-differentiation quantity.
bool IsFinite(const Jet& value) {
  if (!std::isfinite(value.value) || !std::isfinite(value.directional_first) ||
      !std::isfinite(value.directional_second)) {
    return false;
  }
  return std::all_of(value.mixed_first.begin(), value.mixed_first.end(),
                     [](const double coefficient) {
                       return std::isfinite(coefficient);
                     }) &&
         std::all_of(value.mixed_second.begin(), value.mixed_second.end(),
                     [](const double coefficient) {
                       return std::isfinite(coefficient);
                     }) &&
         std::all_of(value.first.begin(), value.first.end(),
                     [](const double coefficient) {
                       return std::isfinite(coefficient);
                     }) &&
         std::all_of(value.second.begin(), value.second.end(),
                     [](const double coefficient) {
                       return std::isfinite(coefficient);
                     });
}

/// @brief Forms a Jet vector from a fixed reference vector.
JetVector3 ConstantVector(const Vector3& value) {
  return {Constant(value.X()), Constant(value.Y()), Constant(value.Z())};
}

/// @brief Adds two Jet vectors.
JetVector3 Add(const JetVector3& left, const JetVector3& right) {
  return {left[0U] + right[0U], left[1U] + right[1U], left[2U] + right[2U]};
}

/// @brief Subtracts two Jet vectors.
JetVector3 Subtract(const JetVector3& left, const JetVector3& right) {
  return {left[0U] - right[0U], left[1U] - right[1U], left[2U] - right[2U]};
}

/// @brief Scales a Jet vector.
JetVector3 Scale(const Jet& scalar, const JetVector3& vector) {
  return {scalar * vector[0U], scalar * vector[1U], scalar * vector[2U]};
}

/// @brief Scales a Jet vector by a fixed value.
JetVector3 Scale(const double scalar, const JetVector3& vector) {
  return {scalar * vector[0U], scalar * vector[1U], scalar * vector[2U]};
}

/// @brief Computes the Jet dot product.
Jet Dot(const JetVector3& left, const JetVector3& right) {
  return left[0U] * right[0U] + left[1U] * right[1U] + left[2U] * right[2U];
}

/// @brief Computes the Jet right-handed cross product.
JetVector3 Cross(const JetVector3& left, const JetVector3& right) {
  return {left[1U] * right[2U] - left[2U] * right[1U],
          left[2U] * right[0U] - left[0U] * right[2U],
          left[0U] * right[1U] - left[1U] * right[0U]};
}

/// @brief Seeds only explicitly zero matrix entries with constant provenance.
JetMatrix3 ConstantZeroMatrix() {
  JetMatrix3 result{};
  for (auto& row : result)
    for (auto& entry : row) entry = Constant(0.0);
  return result;
}

/// @brief Builds a Jet identity matrix.
JetMatrix3 Identity() {
  JetMatrix3 result = ConstantZeroMatrix();
  for (std::size_t index = 0U; index < 3U; ++index) {
    result[index][index] = Constant(1.0);
  }
  return result;
}

/// @brief Forms a skew matrix from a Jet vector.
JetMatrix3 Skew(const JetVector3& value) {
  JetMatrix3 result = ConstantZeroMatrix();
  result[0U][1U] = -1.0 * value[2U];
  result[0U][2U] = value[1U];
  result[1U][0U] = value[2U];
  result[1U][2U] = -1.0 * value[0U];
  result[2U][0U] = -1.0 * value[1U];
  result[2U][1U] = value[0U];
  return result;
}

/// @brief Multiplies two three-by-three Jet matrices.
JetMatrix3 Multiply(const JetMatrix3& left, const JetMatrix3& right) {
  return math_internal::MultiplyJetMatrix(left, right, Constant);
}

/// @brief Multiplies a Jet matrix by a Jet vector.
JetVector3 Multiply(const JetMatrix3& matrix, const JetVector3& vector) {
  return math_internal::MultiplyJetMatrixVector(matrix, vector, Constant);
}

/// @brief Transposes a three-by-three Jet matrix.
JetMatrix3 Transpose(const JetMatrix3& matrix) {
  JetMatrix3 result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row][column] = matrix[column][row];
    }
  }
  return result;
}

/// @brief Returns a scalar multiple of a Jet matrix.
JetMatrix3 Scale(const double scalar, const JetMatrix3& matrix) {
  JetMatrix3 result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row][column] = scalar * matrix[row][column];
    }
  }
  return result;
}

/// @brief Returns a trial-coordinate scalar multiple of a Jet matrix.
JetMatrix3 Scale(const Jet& scalar, const JetMatrix3& matrix) {
  JetMatrix3 result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row][column] = scalar * matrix[row][column];
    }
  }
  return result;
}

/// @brief Adds two three-by-three Jet matrices.
JetMatrix3 Add(const JetMatrix3& left, const JetMatrix3& right) {
  JetMatrix3 result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row][column] = left[row][column] + right[row][column];
    }
  }
  return result;
}

/// @brief Reads shell variations without conflating directional third order.
double ReadChannel(const Jet& jet, std::size_t degree, std::size_t channel) {
  if (degree == 0) return jet.value;
  if (degree == 1) return channel < kGlobalDofCount ? jet.first[channel]
                                                   : jet.directional_first;
  if (degree == 3) return jet.mixed_second[channel];
  if (channel < kHessianSize) return jet.second.ReadData()[channel];
  if (channel == kHessianSize) return jet.directional_second;
  return jet.mixed_first[channel - kHessianSize - 1];
}
void WriteChannel(Jet& jet, std::size_t degree, std::size_t channel,
                   double value) {
  if (degree == 0) { jet.value = value; return; }
  // A null Hessian is not constant provenance, including zero first seeds.
  jet.constant_derivatives = false;
  if (degree == 1) {
    if (channel < kGlobalDofCount) jet.first[channel] = value;
    else jet.directional_first = value;
  } else if (degree == 3) jet.mixed_second[channel] = value;
  else if (channel < kHessianSize) jet.second.OverwriteData()[channel] = value;
  else if (channel == kHessianSize) jet.directional_second = value;
  else jet.mixed_first[channel - kHessianSize - 1] = value;
}

/// @brief Evaluates Rodrigues' map with an analytic near-zero series branch.
JetMatrix3 ExpSkew(const JetVector3& rotation, const bool increment_only = false) {
  const JetMatrix3 identity = increment_only ? JetMatrix3{} : Identity();
  const Jet squared_norm = Dot(rotation, rotation);
  Jet sine_over_norm{};
  Jet one_minus_cosine_over_norm_squared{};
  if (squared_norm.value < 1.0e-8) {
    const Jet squared = squared_norm * squared_norm;
    const Jet cubed = squared * squared_norm;
    sine_over_norm = Constant(1.0) - squared_norm / Constant(6.0) +
                     squared / Constant(120.0) - cubed / Constant(5040.0);
    one_minus_cosine_over_norm_squared =
        Constant(0.5) - squared_norm / Constant(24.0) +
        squared / Constant(720.0) - cubed / Constant(40320.0);
  } else {
    const Jet norm = Sqrt(squared_norm);
    sine_over_norm = Sin(norm) / norm;
    one_minus_cosine_over_norm_squared =
        (Constant(1.0) - Cos(norm)) / squared_norm;
  }
  const JetMatrix3 skew = Skew(rotation);
  JetMatrix3 result =
      Add(Add(identity, Scale(sine_over_norm, skew)),
          Scale(one_minus_cosine_over_norm_squared, Multiply(skew, skew)));
  if (squared_norm.value < 0.25) {
    // Near the original small-angle boundary, cancellation in 1-cos(theta)
    // also corrupts the physical-force gradient and Hessian. Use the existing
    // stable series consistently for the complete Jet and rotation tensor.
    Jet sine_series = Constant(1.0);
    Jet cosine_series = Constant(0.5);
    Jet power = Constant(1.0);
    double sine_coefficient = 1.0;
    double cosine_coefficient = 0.5;
    for (unsigned int k = 1U; k <= 9U; ++k) {
      power = power * squared_norm;
      sine_coefficient /= -static_cast<double>((2U * k) * (2U * k + 1U));
      cosine_coefficient /= -static_cast<double>((2U * k + 1U) * (2U * k + 2U));
      sine_series = sine_series + sine_coefficient * power;
      cosine_series = cosine_series + cosine_coefficient * power;
    }
    const JetMatrix3 stable = Add(Add(identity, Scale(sine_series, skew)),
                                  Scale(cosine_series, Multiply(skew, skew)));
    result = stable;
  }
  return result;
}

/// @brief Stores the complete Jet-valued global element trial coordinates.
struct JetTrial {
  std::array<JetVector3, kNodeCount> displacement;
  std::array<JetVector3, kNodeCount> rotation;
  std::array<JetVector3, kNodeCount> director_increment;
};

/// @brief Lifts one validated nonlinear trial into independent Jet variables.
JetTrial MakeJetTrial(const NonlinearMitc4TrialDofs& trial,
                      const std::array<Matrix, kNodeCount>& initial_frames,
                      const Vector* direction) {
  JetTrial result{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    for (std::size_t component = 0U; component < 3U; ++component) {
      result.displacement[node][component] = Variable(
          trial.displacement[node][component], 6U * node + component,
          direction == nullptr ? 0.0 : (*direction)[6U * node + component]);
      result.rotation[node][component] = Variable(
          trial.rotation_vector[node][component], 6U * node + 3U + component,
          direction == nullptr ? 0.0
                               : (*direction)[6U * node + 3U + component]);
    }
    JetMatrix3 frame{};
    for (std::size_t row = 0U; row < 3U; ++row) {
      for (std::size_t column = 0U; column < 3U; ++column) {
        frame[row][column] = Constant(initial_frames[node](row, column));
      }
    }
    const JetVector3 initial_director{frame[0U][2U], frame[1U][2U],
                                      frame[2U][2U]};
    result.director_increment[node] =
        Multiply(ExpSkew(result.rotation[node], true), initial_director);
  }
  return result;
}

/// @brief Forms covariant increments directly, preserving exact zero motion.
std::array<JetVector3, 3> CovariantIncrement(
    const std::array<Vector3, kNodeCount>& initial_directors,
    const JetTrial& trial, const double thickness, const double xi,
    const double eta, const double zeta) {
  static_cast<void>(initial_directors);
  const auto shape = Mitc4Shell::ShapeFunctions(xi, eta);
  std::array<JetVector3, 3> increment{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const auto displacement =
        Subtract(trial.displacement[node], trial.displacement[0U]);
    const auto director = trial.director_increment[node];
    const auto layer =
        Add(displacement, Scale(0.5 * thickness * zeta, director));
    increment[0U] =
        Add(increment[0U], Scale(shape.xi_derivatives[node], layer));
    increment[1U] =
        Add(increment[1U], Scale(shape.eta_derivatives[node], layer));
    increment[2U] = Add(increment[2U],
                        Scale(0.5 * thickness * shape.values[node], director));
  }
  return increment;
}

/// @brief Forms the direct covariant Green--Lagrange strain tensor.
JetTensor DirectStrain(const std::array<JetVector3, 3>& increment,
                       const Geometry& reference) {
  JetTensor result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row][column] =
          0.5 * (Dot(ConstantVector(reference.reference_covariant[row]),
                     increment[column]) +
                 Dot(increment[row],
                     ConstantVector(reference.reference_covariant[column])) +
                 Dot(increment[row], increment[column]));
    }
  }
  result[2U][2U] = Constant(0.0);
  return result;
}

/// @brief Reconstructs a fixed-reference Cartesian Jet strain tensor.
JetTensor CartesianStrain(const JetTensor& covariant,
                          const Geometry& reference) {
  // E_cart = B^-T E_cov B^-1. First solve each row against B, then each
  // column against B^T, carrying every independent Jet channel unchanged.
  const std::array<std::size_t, 4> channels{
      1, kGlobalDofCount+1, kHessianSize+1+kGlobalDofCount, kGlobalDofCount};
  JetTensor right{}, result{};
  for (std::size_t row = 0; row < 3; ++row) {
    if (!math_internal::SolveConstantJetAction(
        covariant[row], *reference.reference_factor, true, channels,
        Constant, ReadChannel, WriteChannel, right[row])) {
      for (auto& entry : result)
        entry.fill(Constant(std::numeric_limits<double>::quiet_NaN()));
      return result;
    }
  }
  for (std::size_t column = 0; column < 3; ++column) {
    const JetVector3 rhs{right[0][column],right[1][column],right[2][column]};
    JetVector3 action{};
    if (!math_internal::SolveConstantJetAction(
        rhs, *reference.reference_factor, true, channels,
        Constant, ReadChannel, WriteChannel, action)) {
      for (auto& entry : result)
        entry.fill(Constant(std::numeric_limits<double>::quiet_NaN()));
      return result;
    }
    for (std::size_t row = 0; row < 3; ++row) result[row][column] = action[row];
  }
  return result;
}

/// @brief Forms one fixed-frame engineering strain component.
Jet FrameComponent(const JetVector3& left, const JetTensor& tensor,
                   const JetVector3& right) {
  Jet result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result = result + left[row] * tensor[row][column] * right[column];
    }
  }
  return result;
}

/// @brief Evaluates the physical strain energy through exact Jet variations.
Jet PhysicalEnergy(const std::array<Vector3, kNodeCount>& coordinates,
                   const std::array<Vector3, kNodeCount>& initial_directors,
                   const JetTrial& trial, const double thickness,
                   const double youngs, const double poisson, bool& valid) {
  const NonlinearMitc4TrialDofs reference_trial{};
  constexpr std::array<std::array<double, 2>, 4> tying_locations{
      std::array<double, 2>{0.0, -1.0}, std::array<double, 2>{0.0, 1.0},
      std::array<double, 2>{-1.0, 0.0}, std::array<double, 2>{1.0, 0.0}};
  std::array<JetTensor, 4> tying{};
  for (std::size_t point = 0U; point < tying_locations.size(); ++point) {
    Geometry reference{};
    if (!EvaluateGeometry(coordinates, initial_directors, initial_directors,
                          reference_trial, thickness,
                          tying_locations[point][0U],
                          tying_locations[point][1U], 0.0, reference)) {
      valid = false;
      return {};
    }
    tying[point] =
        DirectStrain(CovariantIncrement(initial_directors, trial, thickness,
                                        tying_locations[point][0U],
                                        tying_locations[point][1U], 0.0),
                     reference);
  }

  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}};
  const std::array<double, 2> thickness_points{-gauss, gauss};
  const double plane_stress = youngs / (1.0 - poisson * poisson);
  const double shear = kShearCorrection * youngs / (2.0 * (1.0 + poisson));
  Jet energy{};
  for (const auto& surface : surface_points) {
    for (const double zeta : thickness_points) {
      Geometry reference{};
      if (!EvaluateGeometry(coordinates, initial_directors, initial_directors,
                            reference_trial, thickness, surface[0U],
                            surface[1U], zeta, reference)) {
        valid = false;
        return {};
      }
      JetTensor covariant =
          DirectStrain(CovariantIncrement(initial_directors, trial, thickness,
                                          surface[0U], surface[1U], zeta),
                       reference);
      const Mitc4TyingWeights weights =
          Mitc4Shell::TyingWeights(surface[0U], surface[1U]);
      covariant[0U][2U] = weights.xi_zeta[0U] * tying[0U][0U][2U] +
                          weights.xi_zeta[1U] * tying[1U][0U][2U];
      covariant[2U][0U] = covariant[0U][2U];
      covariant[1U][2U] = weights.eta_zeta[0U] * tying[2U][1U][2U] +
                          weights.eta_zeta[1U] * tying[3U][1U][2U];
      covariant[2U][1U] = covariant[1U][2U];
      const JetTensor cartesian = CartesianStrain(covariant, reference);
      const JetVector3 e1 = ConstantVector(Vector3{reference.local_frame.e1});
      const JetVector3 e2 = ConstantVector(Vector3{reference.local_frame.e2});
      const JetVector3 e3 = ConstantVector(Vector3{reference.local_frame.e3});
      const std::array<Jet, 5> strain{FrameComponent(e1, cartesian, e1),
                                      FrameComponent(e2, cartesian, e2),
                                      2.0 * FrameComponent(e1, cartesian, e2),
                                      2.0 * FrameComponent(e1, cartesian, e3),
                                      2.0 * FrameComponent(e2, cartesian, e3)};
      const std::array<Jet, 5> stress{
          plane_stress * (strain[0U] + poisson * strain[1U]),
          plane_stress * (poisson * strain[0U] + strain[1U]),
          plane_stress * (1.0 - poisson) * 0.5 * strain[2U], shear * strain[3U],
          shear * strain[4U]};
      Jet density{};
      for (std::size_t component = 0U; component < strain.size(); ++component) {
        density = density + strain[component] * stress[component];
      }
      energy = energy + 0.5 * reference.reference_jacobian * density;
    }
  }
  if (!IsFinite(energy)) {
    valid = false;
  }
  return energy;
}

/// @brief Applies one polar iteration for drill corotation.
JetMatrix3 PolarRotation(const JetMatrix3& deformation) {
  JetMatrix3 result = deformation;
  for (std::size_t iteration = 0U; iteration < 1U; ++iteration) {
    const JetMatrix3 coefficients = Transpose(result);
    Matrix primal{3, 3};
    for (std::size_t row = 0; row < 3; ++row)
      for (std::size_t column = 0; column < 3; ++column)
        primal(row, column) = coefficients[row][column].value;
    const auto factor = dense_blas_internal::DenseFactor::Lu(primal);
    JetMatrix3 candidate{};
    for (std::size_t column = 0; column < 3; ++column) {
      JetVector3 rhs{Constant(0), Constant(0), Constant(0)}, action{};
      rhs[column] = Constant(1);
      if (!factor.HasValue() || !math_internal::SolveJetAction(
          coefficients, rhs, factor.Value(),
          {1, kGlobalDofCount + 1, kHessianSize + 1 + kGlobalDofCount,
           kGlobalDofCount}, Constant, ReadChannel, WriteChannel, action)) {
        for (auto& row : candidate)
          row.fill(Constant(std::numeric_limits<double>::quiet_NaN()));
        return candidate;
      }
      for (std::size_t row = 0; row < 3; ++row)
        candidate[row][column] = 0.5 * (result[row][column] + action[row]);
    }
    // All actions used the previous snapshot; never a partial next iterate.
    result = std::move(candidate);
  }
  return result;
}

/// @brief Returns the unnormalized physical drill-null chart vector.
JetVector3 DrillNullPreimage(const JetVector3& phi, const Vector3& director) {
  const Jet square = Dot(phi, phi);
  Jet coefficient{};
  if (square.value < 0.25) {
    constexpr std::array<double, 9> series{1.0 / 12.0,
                                           1.0 / 720.0,
                                           1.0 / 30240.0,
                                           1.0 / 1209600.0,
                                           1.0 / 47900160.0,
                                           691.0 / 1307674368000.0,
                                           1.0 / 74724249600.0,
                                           3617.0 / 10670622842880000.0,
                                           43867.0 / 5109094217170944000.0};
    coefficient = Constant(series.back());
    for (std::size_t i = series.size() - 1U; i > 0U; --i)
      coefficient = Constant(series[i - 1U]) + square * coefficient;
  } else {
    const Jet half = 0.5 * Sqrt(square);
    coefficient = (Constant(1.0) - half * Cos(half) * Reciprocal(Sin(half))) *
                  Reciprocal(square);
  }
  const auto d0 = ConstantVector(director);
  const auto cross = Cross(phi, d0);
  const auto inverse_action =
      Add(Add(d0, Scale(0.5, cross)), Scale(coefficient, Cross(phi, cross)));
  return inverse_action;
}

/// @brief Returns the normalized physical drill-null chart vector.
JetVector3 DrillNull(const JetVector3& phi, const Vector3& director) {
  const JetVector3 preimage = DrillNullPreimage(phi, director);
  return Scale(Reciprocal(Sqrt(Dot(preimage, preimage))), preimage);
}

/// @brief Returns the speed of physical spin for a unit chart null vector.
Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) {
  const JetVector3 preimage = DrillNullPreimage(phi, director);
  return Reciprocal(Sqrt(Dot(preimage, preimage)));
}

/// @brief Evaluates the objective numerical drilling potential and variations.
Jet DrillingEnergy(const std::array<Vector3, kNodeCount>& coordinates,
                   const std::array<Vector3, kNodeCount>& initial_directors,
                   const JetTrial& trial, const double thickness,
                   const double drilling_stiffness, bool& valid, Matrix* metric,
                   std::array<Jet, 4>* constraints) {
  const NonlinearMitc4TrialDofs reference_trial{};
  Geometry reference{};
  if (!EvaluateGeometry(coordinates, initial_directors, initial_directors,
                        reference_trial, thickness, 0.0, 0.0, 0.0, reference)) {
    valid = false;
    return {};
  }
  const auto current =
      CovariantIncrement(initial_directors, trial, thickness, 0.0, 0.0, 0.0);
  JetMatrix3 deformation = Identity();
  const std::array<std::size_t, 4> channels{
      1, kGlobalDofCount+1, kHessianSize+1+kGlobalDofCount, kGlobalDofCount};
  for (std::size_t row = 0; row < 3; ++row) {
    const JetVector3 rhs{current[0][row],current[1][row],current[2][row]};
    JetVector3 action{};
    if (!math_internal::SolveConstantJetAction(
        rhs, *reference.reference_factor, true, channels,
        Constant, ReadChannel, WriteChannel, action)) {
      valid = false;
      return {};
    }
    for (std::size_t column = 0; column < 3; ++column)
      deformation[row][column] = deformation[row][column] + action[column];
  }
  const JetMatrix3 element_rotation = PolarRotation(deformation);
  Jet energy{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const JetMatrix3 nodal_rotation = ExpSkew(trial.rotation[node]);
    const JetMatrix3 corotated =
        Multiply(Transpose(element_rotation), nodal_rotation);
    const JetVector3 initial_director = ConstantVector(initial_directors[node]);
    const JetVector3 swung_director = Multiply(corotated, initial_director);
    const Jet alignment = Dot(initial_director, swung_director);
    if (!std::isfinite(alignment.value) ||
        !(alignment.value > -1.0 + kDirectorOppositionTolerance)) {
      valid = false;
      return {};
    }
    const JetVector3 swing_axis = Cross(initial_director, swung_director);
    const JetMatrix3 swing_skew = Skew(swing_axis);
    const JetMatrix3 swing = Add(Add(Identity(), swing_skew),
                                 Scale(Reciprocal(Constant(1.0) + alignment),
                                       Multiply(swing_skew, swing_skew)));
    const JetMatrix3 twist = Multiply(Transpose(swing), corotated);
    const JetVector3 twist_vex{twist[2U][1U] - twist[1U][2U],
                               twist[0U][2U] - twist[2U][0U],
                               twist[1U][0U] - twist[0U][1U]};
    const Jet sine = 0.5 * Dot(initial_director, twist_vex);
    const Jet cosine =
        0.5 * (twist[0U][0U] + twist[1U][1U] + twist[2U][2U] - Constant(1.0));
    const Jet gamma = Atan2(sine, cosine);
    if (constraints != nullptr) {
      // n^T grad(Uphys)=0 and Dgamma_j[n_i]=delta_ij/||Jl^-1 d||.
      (*constraints)[node] =
          drilling_stiffness *
          DrillSpinRate(trial.rotation[node], initial_directors[node]) * gamma;
    }
    energy = energy + 0.5 * drilling_stiffness * gamma * gamma;
    std::array<double, kGlobalDofCount> scaled_gradient{};
    for (std::size_t row = 0U; row < kGlobalDofCount; ++row)
      scaled_gradient[row] = drilling_stiffness * gamma.first[row];
    std::array<double, kHessianSize> contribution;
    Outer(scaled_gradient, gamma.first, contribution.data());
    for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
      for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
        (*metric)(row, column) += contribution[row * kGlobalDofCount + column];
      }
    }
  }
  if (!IsFinite(energy)) {
    valid = false;
  }
  return energy;
}

}  // namespace

NonlinearMitc4Shell::NonlinearMitc4Shell(
    std::array<Vector3, 4> coordinates, std::array<Matrix, 4> initial_frames,
    std::array<Vector3, 4> initial_directors, const double thickness,
    const double youngs_modulus, const double poisson_ratio,
    const double drilling_stiffness, SourceLocation source_location,
    std::string identity)
    : coordinates_{std::move(coordinates)},
      initial_frames_{std::move(initial_frames)},
      initial_directors_{std::move(initial_directors)},
      thickness_{thickness},
      youngs_modulus_{youngs_modulus},
      poisson_ratio_{poisson_ratio},
      drilling_stiffness_{drilling_stiffness},
      source_location_{std::move(source_location)},
      identity_{std::move(identity)} {}

Result<NonlinearMitc4Shell> NonlinearMitc4Shell::Create(
    std::array<const Node*, 4> nodes,
    std::array<ShellNodeInitialFrame, 4> initial_frames,
    const ShellSection& section, const LinearElasticMaterial& material) {
  const std::string identity = ElementIdentity(nodes);
  if (std::any_of(nodes.begin(), nodes.end(),
                  [](const Node* node) { return node == nullptr; })) {
    return Failure<NonlinearMitc4Shell>(
        "invalid-nonlinear-shell-geometry", section.location, identity,
        "Nonlinear 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 Failure<NonlinearMitc4Shell>(
          "invalid-nonlinear-shell-geometry", nodes[node]->location, identity,
          "Nonlinear MITC4 node coordinates must be finite.");
    }
    for (std::size_t previous = 0U; previous < node; ++previous) {
      if (coordinates[node] == coordinates[previous]) {
        return Failure<NonlinearMitc4Shell>(
            "invalid-nonlinear-shell-geometry", nodes[node]->location, identity,
            "Nonlinear MITC4 nodes must have distinct coordinates.");
      }
      if (initial_frames[node].node_index ==
          initial_frames[previous].node_index) {
        return Failure<NonlinearMitc4Shell>(
            "invalid-nonlinear-shell-frame", nodes[node]->location, identity,
            "Nonlinear MITC4 initial frames must identify four distinct "
            "nodes.");
      }
    }
  }

  std::array<Matrix, kNodeCount> frames{Matrix{3U, 3U}, Matrix{3U, 3U},
                                        Matrix{3U, 3U}, Matrix{3U, 3U}};
  std::array<Vector3, kNodeCount> directors{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    frames[node] = FrameMatrix(initial_frames[node]);
    const auto state = MakeFiniteRotationState(Vector3{}, frames[node]);
    if (!state.HasValue()) {
      return Failure<NonlinearMitc4Shell>(
          "invalid-nonlinear-shell-frame", nodes[node]->location, identity,
          "Nonlinear MITC4 initial frames must be finite proper triads.");
    }
    directors[node] = state.Value().director;
    if ((directors[node] - Vector3{initial_frames[node].director}).Norm() >
        kFrameTolerance) {
      return Failure<NonlinearMitc4Shell>(
          "invalid-nonlinear-shell-frame", nodes[node]->location, identity,
          "Nonlinear MITC4 initial frame director must be its third axis.");
    }
  }

  const Mitc4ShapeFunctions center = Mitc4Shell::ShapeFunctions(0.0, 0.0);
  const Vector3 center_xi = DerivativeSum(center.xi_derivatives, coordinates);
  const Vector3 center_eta = DerivativeSum(center.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 Failure<NonlinearMitc4Shell>(
        "invalid-nonlinear-shell-geometry", nodes[0U]->location, identity,
        "Nonlinear MITC4 center surface basis must be finite and nonzero.");
  }
  const Vector3 positive_face = center_area / center_measure;
  if (std::any_of(directors.begin(), directors.end(),
                  [&positive_face](const Vector3& director) {
                    return !(positive_face.Dot(director) > 0.0);
                  })) {
    return Failure<NonlinearMitc4Shell>(
        "invalid-nonlinear-shell-frame", nodes[0U]->location, identity,
        "Nonlinear MITC4 initial directors must follow the source positive "
        "face.");
  }

  if (!std::isfinite(section.thickness) || !(section.thickness > 0.0)) {
    return Failure<NonlinearMitc4Shell>(
        "invalid-nonlinear-shell-section", section.location, identity,
        "Nonlinear 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 Failure<NonlinearMitc4Shell>(
        "invalid-nonlinear-shell-material", material.location, identity,
        "Nonlinear MITC4 requires finite E>0 and -1<nu<0.5.");
  }

  std::array<std::array<double, 3>, kNodeCount> director_components{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    director_components[node] = directors[node].Components();
  }
  const auto linear_shell =
      Mitc4Shell::Create(nodes, director_components, section, material);
  if (!linear_shell.HasValue()) {
    return Result<NonlinearMitc4Shell>::Failure(linear_shell.GetStatus());
  }
  const auto linear_stiffness = linear_shell.Value().Stiffness();
  if (!linear_stiffness.HasValue() ||
      !std::isfinite(linear_stiffness.Value().drilling_stiffness) ||
      !(linear_stiffness.Value().drilling_stiffness > 0.0)) {
    return Failure<NonlinearMitc4Shell>(
        "invalid-nonlinear-shell-drilling", section.location, identity,
        "Nonlinear MITC4 requires a finite positive fixed drilling stiffness.");
  }

  NonlinearMitc4Shell shell{coordinates,
                            std::move(frames),
                            directors,
                            section.thickness,
                            material.youngs_modulus,
                            material.poisson_ratio,
                            linear_stiffness.Value().drilling_stiffness,
                            nodes[0U]->location,
                            identity};
  const auto reference = shell.EvaluateKinematics(NonlinearMitc4TrialDofs{});
  if (!reference.HasValue()) {
    return Result<NonlinearMitc4Shell>::Failure(reference.GetStatus());
  }
  return Result<NonlinearMitc4Shell>::Success(std::move(shell));
}

Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(
    const NonlinearMitc4TrialDofs& trial) const {
  std::array<Vector3, kNodeCount> directors{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    if (!trial.displacement[node].IsFinite() ||
        !trial.rotation_vector[node].IsFinite()) {
      return Failure<NonlinearMitc4Recovery>(
          "invalid-nonlinear-shell-trial", source_location_, identity_,
          "Nonlinear MITC4 trial displacements and rotations must be finite.");
    }
    const auto state = MakeFiniteRotationState(trial.rotation_vector[node],
                                               initial_frames_[node]);
    if (!state.HasValue()) {
      return Failure<NonlinearMitc4Recovery>(
          "invalid-nonlinear-shell-rotation", source_location_, identity_,
          "Nonlinear MITC4 trial rotations must remain in the principal "
          "chart.");
    }
    directors[node] = state.Value().director;
  }
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    for (std::size_t other = 0U; other < node; ++other) {
      if (!(directors[node].Dot(directors[other]) >
            -1.0 + kDirectorOppositionTolerance)) {
        return Failure<NonlinearMitc4Recovery>(
            "invalid-nonlinear-shell-director", source_location_, identity_,
            "Nonlinear MITC4 trial directors must not be opposed.");
      }
    }
  }

  Geometry center{};
  if (!EvaluateGeometry(coordinates_, initial_directors_, directors, trial,
                        thickness_, 0.0, 0.0, 0.0, center)) {
    return Failure<NonlinearMitc4Recovery>(
        "invalid-nonlinear-shell-geometry", source_location_, identity_,
        "Nonlinear MITC4 center reference or trial geometry is invalid.");
  }

  constexpr std::array<std::array<double, 2>, 4> tying_locations{
      std::array<double, 2>{0.0, -1.0}, std::array<double, 2>{0.0, 1.0},
      std::array<double, 2>{-1.0, 0.0}, std::array<double, 2>{1.0, 0.0}};
  std::array<Tensor, 4> tying_samples{};
  for (std::size_t point = 0U; point < tying_locations.size(); ++point) {
    const auto point_status = EvaluateTyingSample(
        coordinates_, initial_directors_, directors, trial, thickness_,
        tying_locations[point][0U], tying_locations[point][1U],
        tying_samples[point]);
    if (point_status != PointStatus::kValid) {
      return Failure<NonlinearMitc4Recovery>(
          point_status == PointStatus::kInvalidGeometry
              ? "invalid-nonlinear-shell-geometry"
              : "invalid-nonlinear-shell-recovery",
          source_location_, identity_,
          "Nonlinear MITC4 tying-point geometry or strain is invalid.");
    }
  }

  const double gauss = 1.0 / std::sqrt(3.0);
  constexpr std::array<ShellMidsurfaceLocation, 4> locations{
      ShellMidsurfaceLocation::kGp1, ShellMidsurfaceLocation::kGp2,
      ShellMidsurfaceLocation::kGp3, ShellMidsurfaceLocation::kGp4};
  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_gauss{-1.0, 1.0};
  constexpr std::array<double, 3> stress_zeta{-1.0, 0.0, 1.0};
  constexpr std::array<ShellSectionPosition, 3> stress_position{
      ShellSectionPosition::kBottom, ShellSectionPosition::kMiddle,
      ShellSectionPosition::kTop};

  NonlinearMitc4Recovery recovery{};
  for (std::size_t surface = 0U; surface < surface_points.size(); ++surface) {
    auto& row = recovery.gauss_points[surface];
    row.location = locations[surface];
    row.natural_coordinates = surface_points[surface];
    Geometry midpoint{};
    if (!EvaluateGeometry(coordinates_, initial_directors_, directors, trial,
                          thickness_, row.natural_coordinates[0U],
                          row.natural_coordinates[1U], 0.0, midpoint)) {
      return Failure<NonlinearMitc4Recovery>(
          "invalid-nonlinear-shell-geometry", source_location_, identity_,
          "Nonlinear MITC4 midsurface recovery geometry is invalid.");
    }
    row.local_frame = midpoint.local_frame;

    for (const double thickness_sign : thickness_gauss) {
      PointEvaluation point{};
      const auto point_status = EvaluatePoint(
          coordinates_, initial_directors_, directors, trial, tying_samples,
          thickness_, youngs_modulus_, poisson_ratio_,
          row.natural_coordinates[0U], row.natural_coordinates[1U],
          thickness_sign * gauss, point);
      if (point_status != PointStatus::kValid) {
        return Failure<NonlinearMitc4Recovery>(
            point_status == PointStatus::kInvalidGeometry
                ? "invalid-nonlinear-shell-geometry"
                : "invalid-nonlinear-shell-recovery",
            source_location_, identity_,
            "Nonlinear MITC4 volume recovery geometry or values are invalid.");
      }
      for (std::size_t component = 0U; component < 3U; ++component) {
        row.generalized_strain[component] += 0.5 * point.strain[component];
        row.generalized_strain[3U + component] +=
            3.0 * thickness_sign * gauss * point.strain[component] / thickness_;
        row.section_resultant[component] +=
            0.5 * thickness_ * point.stress[component];
        row.section_resultant[3U + component] +=
            0.25 * thickness_ * thickness_ * thickness_sign * gauss *
            point.stress[component];
      }
      for (std::size_t component = 0U; component < 2U; ++component) {
        row.generalized_strain[6U + component] +=
            0.5 * point.strain[3U + component];
        row.section_resultant[6U + component] +=
            0.5 * thickness_ * point.stress[3U + component];
      }
      double density = 0.0;
      for (std::size_t component = 0U; component < point.strain.size();
           ++component) {
        density += point.strain[component] * point.stress[component];
      }
      recovery.physical_strain_energy +=
          0.5 * density * point.geometry.reference_jacobian;
    }

    for (std::size_t position = 0U; position < stress_zeta.size(); ++position) {
      PointEvaluation point{};
      const auto point_status = EvaluatePoint(
          coordinates_, initial_directors_, directors, trial, tying_samples,
          thickness_, youngs_modulus_, poisson_ratio_,
          row.natural_coordinates[0U], row.natural_coordinates[1U],
          stress_zeta[position], point);
      if (point_status != PointStatus::kValid) {
        return Failure<NonlinearMitc4Recovery>(
            point_status == PointStatus::kInvalidGeometry
                ? "invalid-nonlinear-shell-geometry"
                : "invalid-nonlinear-shell-recovery",
            source_location_, identity_,
            "Nonlinear MITC4 section recovery geometry or values are invalid.");
      }
      row.stress[position] = {
          stress_position[position],
          stress_zeta[position],
          {point.stress[0U], point.stress[1U], point.stress[2U]}};
    }
  }

  if (!IsFinite(recovery)) {
    return Failure<NonlinearMitc4Recovery>(
        "invalid-nonlinear-shell-recovery", source_location_, identity_,
        "Nonlinear MITC4 physical recovery must remain finite.");
  }
  return Result<NonlinearMitc4Recovery>::Success(std::move(recovery));
}

Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::RecoverState(
    const NonlinearMitc4TrialDofs& trial) const {
  return EvaluateKinematics(trial);
}

Result<NonlinearMitc4ExternalDrillingConstraint>
NonlinearMitc4Shell::EvaluateExternalDrillingConstraint(
    const std::size_t node, const Vector3& phi, const Vector3& phi_dot,
    const Vector3& force, const Vector3& force_rate,
    const Vector3& force_acceleration) const {
  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
  if (node >= kNodeCount || !phi_dot.IsFinite() || !force.IsFinite() ||
      !force_rate.IsFinite() || !force_acceleration.IsFinite()) {
    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
        "invalid-nonlinear-shell-external-constraint", source_location_,
        identity_, "External drilling constraint inputs must be finite and "
                   "identify one source-local node.");
  }
  const auto rotation = MakeFiniteRotationState(phi, initial_frames_[node]);
  if (!rotation.HasValue()) {
    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
        "invalid-nonlinear-shell-rotation", source_location_, identity_,
        "External drilling constraint rotations must remain in the principal "
        "chart.");
  }

  JetVector3 jet_phi{};
  for (std::size_t component = 0U; component < 3U; ++component) {
    jet_phi[component] = Variable(phi[component], component,
                                  phi_dot[component]);
  }
  const JetVector3 null_direction = DrillNull(jet_phi, initial_directors_[node]);
  const Jet scalar = -1.0 * Dot(null_direction, ConstantVector(force));
  const Jet scalar_time =
      -1.0 * Dot(null_direction, ConstantVector(force_rate));
  const Jet scalar_acceleration =
      -1.0 * Dot(null_direction, ConstantVector(force_acceleration));
  if (!IsFinite(scalar) || !IsFinite(scalar_time) ||
      !IsFinite(scalar_acceleration)) {
    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
        "invalid-nonlinear-shell-external-constraint", source_location_,
        identity_, "External drilling constraint derivatives must remain "
                   "finite.");
  }

  NonlinearMitc4ExternalDrillingConstraint result{};
  result.value = scalar.value;
  result.gradient = Vector{3U};
  result.hessian = Matrix{3U, 3U};
  result.gradient_directional_first = Vector{3U};
  result.gradient_directional_second = Vector{3U};
  result.directional_first =
      scalar.directional_first + scalar_time.value;
  result.directional_second = scalar.directional_second +
                              2.0 * scalar_time.directional_first +
                              scalar_acceleration.value;
  for (std::size_t row = 0U; row < 3U; ++row) {
    result.gradient[row] = scalar.first[row];
    result.gradient_directional_first[row] =
        scalar.mixed_first[row] + scalar_time.first[row];
    result.gradient_directional_second[row] =
        scalar.mixed_second[row] +
        2.0 * scalar_time.mixed_first[row] + scalar_acceleration.first[row];
    for (std::size_t column = 0U; column < 3U; ++column) {
      result.hessian(row, column) =
          scalar.second[row * kGlobalDofCount + column];
    }
  }
  if (!std::isfinite(result.value) || !std::isfinite(result.directional_first) ||
      !std::isfinite(result.directional_second)) {
    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
        "invalid-nonlinear-shell-external-constraint", source_location_,
        identity_, "External drilling constraint derivatives must remain "
                   "finite.");
  }
  for (std::size_t row = 0U; row < 3U; ++row) {
    if (!std::isfinite(result.gradient[row]) ||
        !std::isfinite(result.gradient_directional_first[row]) ||
        !std::isfinite(result.gradient_directional_second[row])) {
      return Failure<NonlinearMitc4ExternalDrillingConstraint>(
          "invalid-nonlinear-shell-external-constraint", source_location_,
          identity_, "External drilling constraint derivatives must remain "
                     "finite.");
    }
    for (std::size_t column = 0U; column < 3U; ++column) {
      if (!std::isfinite(result.hessian(row, column))) {
        return Failure<NonlinearMitc4ExternalDrillingConstraint>(
            "invalid-nonlinear-shell-external-constraint", source_location_,
            identity_, "External drilling constraint derivatives must remain "
                       "finite.");
      }
    }
  }
  return Result<NonlinearMitc4ExternalDrillingConstraint>::Success(
      std::move(result));
}

Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
    const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial,
    const Vector* coordinate_direction) {
  // Element-local MKL work uses one worker only when the caller has no TLS
  // override. The guard restores that request before returning to global solve.
  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
  if (coordinate_direction != nullptr) {
    if (coordinate_direction->Size() != kGlobalDofCount) {
      return Failure<NonlinearMitc4InternalEquilibrium>(
          "invalid-nonlinear-shell-equilibrium", element.source_location_,
          element.identity_, "The directional derivative requires 24 DOFs.");
    }
    for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
      if (!std::isfinite((*coordinate_direction)[i])) {
        return Failure<NonlinearMitc4InternalEquilibrium>(
            "invalid-nonlinear-shell-equilibrium", element.source_location_,
            element.identity_, "The directional derivative must be finite.");
      }
    }
  }
  const auto recovery = element.EvaluateKinematics(trial);
  if (!recovery.HasValue()) {
    return Result<NonlinearMitc4InternalEquilibrium>::Failure(
        recovery.GetStatus());
  }

  bool valid = true;
  const JetTrial jet_trial =
      MakeJetTrial(trial, element.initial_frames_, coordinate_direction);
  const Jet physical =
      PhysicalEnergy(element.coordinates_, element.initial_directors_,
                     jet_trial, element.thickness_, element.youngs_modulus_,
                     element.poisson_ratio_, valid);
  Matrix drilling_metric{kGlobalDofCount, kGlobalDofCount};
  std::array<Jet, 4> constraints{};
  const Jet drilling = DrillingEnergy(
      element.coordinates_, element.initial_directors_, jet_trial,
      element.thickness_, element.drilling_stiffness_, valid, &drilling_metric,
      coordinate_direction != nullptr ? &constraints : nullptr);
  const Jet internal_energy = physical + drilling;
  if (!valid || !IsFinite(internal_energy)) {
    return Failure<NonlinearMitc4InternalEquilibrium>(
        "invalid-nonlinear-shell-equilibrium", element.source_location_,
        element.identity_,
        "Nonlinear MITC4 internal residual and tangent must remain finite.");
  }

  Vector internal_residual{kGlobalDofCount};
  Matrix internal_tangent{kGlobalDofCount, kGlobalDofCount};
  for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
    internal_residual[row] = internal_energy.first[row];
    for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
      internal_tangent(row, column) =
          internal_energy.second[row * kGlobalDofCount + column];
    }
  }
  for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
    if (!std::isfinite(internal_residual[row])) {
      return Failure<NonlinearMitc4InternalEquilibrium>(
          "invalid-nonlinear-shell-equilibrium", element.source_location_,
          element.identity_,
          "Nonlinear MITC4 internal residual and tangent must remain finite.");
    }
    for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
      if (!std::isfinite(internal_tangent(row, column))) {
        return Failure<NonlinearMitc4InternalEquilibrium>(
            "invalid-nonlinear-shell-equilibrium", element.source_location_,
            element.identity_,
            "Nonlinear MITC4 internal residual and tangent must remain "
            "finite.");
      }
    }
  }

  NonlinearMitc4InternalEquilibrium equilibrium{
      std::move(internal_residual), std::move(internal_tangent),
      Vector{kGlobalDofCount}, recovery.Value(),
      recovery.Value().physical_strain_energy};
  if (coordinate_direction != nullptr) {
    equilibrium.drilling_constraint = Vector{4U};
    equilibrium.drilling_constraint_gradient = Matrix{4U, kGlobalDofCount};
    for (std::size_t node = 0U; node < 4U; ++node) {
      if (!IsFinite(constraints[node]))
        return Failure<NonlinearMitc4InternalEquilibrium>(
            "invalid-nonlinear-shell-equilibrium", element.source_location_,
            element.identity_,
            "Drilling constraint derivatives must remain finite.");
      equilibrium.drilling_constraint[node] = constraints[node].value;
      equilibrium.drilling_constraint_hessian[node] =
          Matrix{kGlobalDofCount, kGlobalDofCount};
      equilibrium.drilling_constraint_third_directional_gradient[node] =
          Vector{kGlobalDofCount};
      for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
        equilibrium.drilling_constraint_gradient(node, i) =
            constraints[node].first[i];
        equilibrium.drilling_constraint_third_directional_gradient[node][i] =
            constraints[node].mixed_second[i];
        for (std::size_t j = 0U; j < kGlobalDofCount; ++j)
          equilibrium.drilling_constraint_hessian[node](i, j) =
              constraints[node].second[i * kGlobalDofCount + j];
      }
    }
  }
  equilibrium.physical_force = Vector{kGlobalDofCount};
  equilibrium.drilling_force = Vector{kGlobalDofCount};
  equilibrium.drilling_tangent = Matrix{kGlobalDofCount, kGlobalDofCount};
  equilibrium.drilling_metric = std::move(drilling_metric);
  equilibrium.drilling_energy = drilling.value;
  equilibrium.force_directional_second = Vector{kGlobalDofCount};
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    equilibrium.internal_force[index] = equilibrium.internal_residual[index];
    equilibrium.physical_force[index] = physical.first[index];
    equilibrium.drilling_force[index] = drilling.first[index];
    equilibrium.force_directional_second[index] =
        internal_energy.mixed_second[index];
    for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
      equilibrium.drilling_tangent(index, column) =
          drilling.second[index * kGlobalDofCount + column];
    }
  }
  return Result<NonlinearMitc4InternalEquilibrium>::Success(
      std::move(equilibrium));
}

}  // namespace fesa