#include "fesa/elements/nonlinear_mitc4_shell.h"

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

#include "fesa/elements/finite_rotation.h"

namespace fesa {
namespace {

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> reciprocal;
  Mitc4LocalFrame local_frame;
  double reference_jacobian;
};

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

/// @brief Builds one deterministic structured nonlinear-shell failure.
template <class T>
Result<T> Failure(const std::string& code, const SourceLocation& location,
                  const std::string& identity, std::string message) {
  return Result<T>::Failure(Status::Failure(
      FailureCategory::kModel,
      {{Severity::kError, code, location, "*ELEMENT", identity,
        std::move(message)}}));
}

/// @brief Builds the stable source label concatenation used by diagnostics.
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;
}

/// @brief Forms a weighted vector sum in source-node order.
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;
}

/// @brief Forms a stable derivative sum after cancelling the constant term.
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[0U];
  }
  return WeightedSum(derivatives, relative);
}

/// @brief Returns a frame matrix whose columns are a, b and the director.
Matrix FrameMatrix(const ShellNodeInitialFrame& frame) {
  Matrix result{3U, 3U};
  for (std::size_t row = 0U; row < 3U; ++row) {
    result(row, 0U) = frame.tangent_a[row];
    result(row, 1U) = frame.tangent_b[row];
    result(row, 2U) = frame.director[row];
  }
  return result;
}

/// @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 Checks every scalar in one fixed-size engineering quantity.
template <std::size_t N>
bool IsFinite(const std::array<double, N>& values) {
  return std::all_of(values.begin(), values.end(),
                     [](const double value) { return std::isfinite(value); });
}

/// @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 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;

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

  geometry.reciprocal[0U] =
      geometry.reference_covariant[1U].Cross(geometry.reference_covariant[2U]) /
      geometry.reference_jacobian;
  geometry.reciprocal[1U] =
      geometry.reference_covariant[2U].Cross(geometry.reference_covariant[0U]) /
      geometry.reference_jacobian;
  geometry.reciprocal[2U] =
      geometry.reference_covariant[0U].Cross(geometry.reference_covariant[1U]) /
      geometry.reference_jacobian;
  if (!geometry.reciprocal[0U].IsFinite() ||
      !geometry.reciprocal[1U].IsFinite() ||
      !geometry.reciprocal[2U].IsFinite()) {
    return false;
  }
  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.current_covariant[row].Dot(
                     geometry.current_covariant[column]) -
                 geometry.reference_covariant[row].Dot(
                     geometry.reference_covariant[column]));
    }
  }
  strain[2U][2U] = 0.0;
  return strain;
}

/// @brief Reconstructs a fixed-reference Cartesian strain from covariant data.
Tensor CartesianStrain(const Tensor& covariant,
                       const std::array<Vector3, 3>& reciprocal) {
  Tensor cartesian{};
  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) {
          cartesian[row][column] += covariant[first][second] *
                                    reciprocal[first][row] *
                                    reciprocal[second][column];
        }
      }
    }
  }
  return cartesian;
}

/// @brief Evaluates a tensor component in a fixed local reference frame.
double FrameComponent(const Vector3& left, const Tensor& 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 Converts fixed-reference Cartesian strain to engineering components.
EngineeringStrain LocalEngineering(const Tensor& 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)};
}

/// @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 Evaluates direct strain at one mandatory edge-midpoint tying point.
bool 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 false;
  }
  strain = DirectStrain(geometry);
  return IsFinite(strain);
}

/// @brief Computes one MITC-projected local strain and material stress point.
bool 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 false;
  }
  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.reciprocal);
  point.strain = LocalEngineering(cartesian, point.geometry.local_frame);
  point.stress = Stress(point.strain, youngs, poisson);
  return IsFinite(covariant) && IsFinite(cartesian) && IsFinite(point.strain) &&
         IsFinite(point.stress);
}

/// @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() : values_{std::make_shared<std::array<double, kHessianSize>>()} {}

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

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

  auto begin() const { return values_->begin(); }

  auto end() const { return values_->end(); }

 private:
  std::shared_ptr<std::array<double, kHessianSize>> values_;
};

/// @brief Carries exact first and second trial-coordinate variations.
struct Jet {
  double value{};
  std::array<double, kGlobalDofCount> first{};
  Hessian second;
};

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

/// @brief Returns one independent global trial coordinate.
Jet Variable(const double value, const std::size_t index) {
  Jet result = Constant(value);
  result.first[index] = 1.0;
  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;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = left.first[index] + right.first[index];
  }
  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;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = left.first[index] - right.first[index];
  }
  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;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] =
        left.first[index] * right.value + left.value * right.first[index];
  }
  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 +
                             left.first[row] * right.first[column] +
                             left.first[column] * right.first[row] +
                             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;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = scalar * value.first[index];
  }
  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) {
  Jet result{};
  result.value = value;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] = first * input.first[index];
  }
  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] +
                             second * input.first[row] * input.first[column];
    }
  }
  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);
}

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

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

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

/// @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;
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    result.first[index] =
        derivative_y * y.first[index] + derivative_x * x.first[index];
  }
  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] +
                             derivative_yy * y.first[row] * y.first[column] +
                             derivative_yx * (y.first[row] * x.first[column] +
                                              x.first[row] * y.first[column]) +
                             derivative_xx * x.first[row] * x.first[column];
    }
  }
  return result;
}

/// @brief Checks every coefficient in an automatic-differentiation quantity.
bool IsFinite(const Jet& value) {
  if (!std::isfinite(value.value)) {
    return false;
  }
  return 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 Builds a Jet identity matrix.
JetMatrix3 Identity() {
  JetMatrix3 result{};
  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{};
  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) {
  JetMatrix3 result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      for (std::size_t index = 0U; index < 3U; ++index) {
        result[row][column] =
            result[row][column] + left[row][index] * right[index][column];
      }
    }
  }
  return result;
}

/// @brief Multiplies a Jet matrix by a Jet vector.
JetVector3 Multiply(const JetMatrix3& matrix, const JetVector3& vector) {
  JetVector3 result{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row] = result[row] + matrix[row][column] * vector[column];
    }
  }
  return result;
}

/// @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 Computes the analytic three-by-three inverse through its adjugate.
JetMatrix3 Inverse(const JetMatrix3& matrix) {
  JetMatrix3 cofactor{};
  cofactor[0U][0U] =
      matrix[1U][1U] * matrix[2U][2U] - matrix[1U][2U] * matrix[2U][1U];
  cofactor[0U][1U] =
      matrix[1U][2U] * matrix[2U][0U] - matrix[1U][0U] * matrix[2U][2U];
  cofactor[0U][2U] =
      matrix[1U][0U] * matrix[2U][1U] - matrix[1U][1U] * matrix[2U][0U];
  cofactor[1U][0U] =
      matrix[0U][2U] * matrix[2U][1U] - matrix[0U][1U] * matrix[2U][2U];
  cofactor[1U][1U] =
      matrix[0U][0U] * matrix[2U][2U] - matrix[0U][2U] * matrix[2U][0U];
  cofactor[1U][2U] =
      matrix[0U][1U] * matrix[2U][0U] - matrix[0U][0U] * matrix[2U][1U];
  cofactor[2U][0U] =
      matrix[0U][1U] * matrix[1U][2U] - matrix[0U][2U] * matrix[1U][1U];
  cofactor[2U][1U] =
      matrix[0U][2U] * matrix[1U][0U] - matrix[0U][0U] * matrix[1U][2U];
  cofactor[2U][2U] =
      matrix[0U][0U] * matrix[1U][1U] - matrix[0U][1U] * matrix[1U][0U];
  const Jet determinant = matrix[0U][0U] * cofactor[0U][0U] +
                          matrix[0U][1U] * cofactor[0U][1U] +
                          matrix[0U][2U] * cofactor[0U][2U];
  JetMatrix3 result{};
  const Jet inverse_determinant = Reciprocal(determinant);
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      result[row][column] = cofactor[column][row] * inverse_determinant;
    }
  }
  return result;
}

/// @brief Evaluates Rodrigues' map with an analytic near-zero series branch.
JetMatrix3 ExpSkew(const JetVector3& rotation) {
  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);
  return Add(Add(Identity(), Scale(sine_over_norm, skew)),
             Scale(one_minus_cosine_over_norm_squared, Multiply(skew, skew)));
}

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

/// @brief Lifts one validated nonlinear trial into independent Jet variables.
JetTrial MakeJetTrial(const NonlinearMitc4TrialDofs& trial,
                      const std::array<Matrix, kNodeCount>& initial_frames) {
  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);
      result.rotation[node][component] = Variable(
          trial.rotation_vector[node][component], 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[node] =
        Multiply(ExpSkew(result.rotation[node]), initial_director);
  }
  return result;
}

/// @brief Forms one current covariant basis from the exact trial variables.
std::array<JetVector3, 3> CurrentCovariant(
    const std::array<Vector3, kNodeCount>& coordinates, const JetTrial& trial,
    const double thickness, const double xi, const double eta,
    const double zeta) {
  const Mitc4ShapeFunctions shape = Mitc4Shell::ShapeFunctions(xi, eta);
  JetVector3 displacement_xi{};
  JetVector3 displacement_eta{};
  JetVector3 director{};
  JetVector3 director_xi{};
  JetVector3 director_eta{};
  Vector3 reference_xi{};
  Vector3 reference_eta{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    reference_xi =
        reference_xi + shape.xi_derivatives[node] * coordinates[node];
    reference_eta =
        reference_eta + shape.eta_derivatives[node] * coordinates[node];
    displacement_xi = Add(displacement_xi, Scale(shape.xi_derivatives[node],
                                                 trial.displacement[node]));
    displacement_eta = Add(displacement_eta, Scale(shape.eta_derivatives[node],
                                                   trial.displacement[node]));
    director = Add(director, Scale(shape.values[node], trial.director[node]));
    director_xi = Add(director_xi,
                      Scale(shape.xi_derivatives[node], trial.director[node]));
    director_eta = Add(
        director_eta, Scale(shape.eta_derivatives[node], trial.director[node]));
  }
  const double half_thickness = 0.5 * thickness;
  return {Add(Add(ConstantVector(reference_xi), displacement_xi),
              Scale(half_thickness * zeta, director_xi)),
          Add(Add(ConstantVector(reference_eta), displacement_eta),
              Scale(half_thickness * zeta, director_eta)),
          Scale(half_thickness, director)};
}

/// @brief Forms the direct covariant Green--Lagrange strain tensor.
JetTensor DirectStrain(const std::array<JetVector3, 3>& current,
                       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(current[row], current[column]) -
                 Constant(reference.reference_covariant[row].Dot(
                     reference.reference_covariant[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) {
  JetTensor result{};
  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) {
          result[row][column] =
              result[row][column] + reference.reciprocal[first][row] *
                                        reference.reciprocal[second][column] *
                                        covariant[first][second];
        }
      }
    }
  }
  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(CurrentCovariant(coordinates, 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(CurrentCovariant(coordinates, 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 Computes a fixed-iteration analytic polar factor for drill
/// corotation.
JetMatrix3 PolarRotation(const JetMatrix3& deformation) {
  JetMatrix3 result = deformation;
  for (std::size_t iteration = 0U; iteration < 12U; ++iteration) {
    result = Scale(0.5, Add(result, Transpose(Inverse(result))));
  }
  return result;
}

/// @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) {
  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 =
      CurrentCovariant(coordinates, trial, thickness, 0.0, 0.0, 0.0);
  JetMatrix3 deformation{};
  for (std::size_t row = 0U; row < 3U; ++row) {
    for (std::size_t column = 0U; column < 3U; ++column) {
      for (std::size_t alpha = 0U; alpha < 3U; ++alpha) {
        deformation[row][column] =
            deformation[row][column] +
            current[alpha][row] * reference.reciprocal[alpha][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);
    energy = energy + 0.5 * drilling_stiffness * gamma * gamma;
  }
  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) {
    if (!EvaluateTyingSample(coordinates_, initial_directors_, directors, trial,
                             thickness_, tying_locations[point][0U],
                             tying_locations[point][1U], tying_samples[point])) {
      return Failure<NonlinearMitc4Recovery>(
          "invalid-nonlinear-shell-geometry", source_location_, identity_,
          "Nonlinear MITC4 tying-point reference or trial geometry 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{};
      if (!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)) {
        return Failure<NonlinearMitc4Recovery>(
            "invalid-nonlinear-shell-geometry", source_location_, identity_,
            "Nonlinear MITC4 volume recovery geometry is 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_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{};
      if (!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)) {
        return Failure<NonlinearMitc4Recovery>(
            "invalid-nonlinear-shell-geometry", source_location_, identity_,
            "Nonlinear MITC4 section-position recovery geometry is 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<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
    const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial) {
  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_);
  const Jet physical =
      PhysicalEnergy(element.coordinates_, element.initial_directors_,
                     jet_trial, element.thickness_, element.youngs_modulus_,
                     element.poisson_ratio_, valid);
  const Jet drilling = DrillingEnergy(
      element.coordinates_, element.initial_directors_, jet_trial,
      element.thickness_, element.drilling_stiffness_, valid);
  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};
  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
    equilibrium.internal_force[index] = equilibrium.internal_residual[index];
  }
  return Result<NonlinearMitc4InternalEquilibrium>::Success(
      std::move(equilibrium));
}

}  // namespace fesa