diff --git a/include/fesa/elements/nonlinear_mitc4_shell.h b/include/fesa/elements/nonlinear_mitc4_shell.h
index e9fcb73..bab0476 100644
--- a/include/fesa/elements/nonlinear_mitc4_shell.h
+++ b/include/fesa/elements/nonlinear_mitc4_shell.h
@@ -7,6 +7,7 @@
 #include "fesa/core/status.h"
 #include "fesa/elements/mitc4_shell.h"
 #include "fesa/math/matrix.h"
+#include "fesa/math/vector.h"
 #include "fesa/math/vector3.h"
 #include "fesa/model/model_types.h"
 #include "fesa/results/result_records.h"
@@ -35,6 +36,26 @@ struct NonlinearMitc4Recovery {
   double physical_strain_energy;
 };
 
+/// @brief Stores complete internal equilibrium data for one nonlinear trial.
+struct NonlinearMitc4InternalEquilibrium {
+  Vector internal_residual;
+  Matrix internal_tangent;
+  Vector internal_force;
+  NonlinearMitc4Recovery recovery;
+  double physical_strain_energy;
+};
+
+class NonlinearMitc4Shell;
+
+/// @brief Evaluates the internal residual and tangent at one complete trial.
+/// @param element Immutable nonlinear MITC4 reference evaluator.
+/// @param trial Complete global translation and rotation trial coordinates.
+/// @return Internal force, consistent tangent and physical recovery, or
+/// failure.
+/// @note External load vectors and their tangents are intentionally excluded.
+Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
+    const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial);
+
 /// @brief Evaluates fixed-reference Total Lagrangian MITC4 kinematics.
 /// @note This class owns no mutable trial or accepted analysis state.
 class NonlinearMitc4Shell {
@@ -68,8 +89,12 @@ class NonlinearMitc4Shell {
                       std::array<Matrix, 4> initial_frames,
                       std::array<Vector3, 4> initial_directors,
                       double thickness, double youngs_modulus,
-                      double poisson_ratio, SourceLocation source_location,
-                      std::string identity);
+                      double poisson_ratio, double drilling_stiffness,
+                      SourceLocation source_location, std::string identity);
+
+  friend Result<NonlinearMitc4InternalEquilibrium>
+  EvaluateInternalResidualAndTangent(const NonlinearMitc4Shell& element,
+                                     const NonlinearMitc4TrialDofs& trial);
 
   std::array<Vector3, 4> coordinates_;
   std::array<Matrix, 4> initial_frames_;
@@ -77,6 +102,7 @@ class NonlinearMitc4Shell {
   double thickness_;
   double youngs_modulus_;
   double poisson_ratio_;
+  double drilling_stiffness_;
   SourceLocation source_location_;
   std::string identity_;
 };
diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp
index 0fa7548..48d117a 100644
--- a/src/fesa/elements/nonlinear_mitc4_shell.cpp
+++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp
@@ -4,6 +4,8 @@
 #include <array>
 #include <cmath>
 #include <cstddef>
+#include <limits>
+#include <memory>
 #include <string>
 #include <utility>
 
@@ -343,19 +345,675 @@ bool IsFinite(const NonlinearMitc4Recovery& recovery) {
   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,
-    SourceLocation source_location, std::string identity)
+    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)} {}
 
@@ -448,12 +1106,31 @@ Result<NonlinearMitc4Shell> NonlinearMitc4Shell::Create(
         "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{});
@@ -611,4 +1288,67 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::RecoverState(
   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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 69d92e9..d94a2a1 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -21,6 +21,7 @@ add_executable(
     unit/elements/finite_rotation_test.cpp
     unit/elements/mitc4_shell_test.cpp
     unit/elements/nonlinear_mitc4_shell_test.cpp
+    unit/elements/nonlinear_mitc4_shell_tangent_test.cpp
     unit/fem/dof_manager_test.cpp
     unit/math/dense_blas_internal_test.cpp
     unit/math/matrix_test.cpp
diff --git a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp
new file mode 100644
index 0000000..2379f1e
--- /dev/null
+++ b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp
@@ -0,0 +1,420 @@
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+#include <string>
+
+#include "fesa/elements/finite_rotation.h"
+#include "fesa/elements/nonlinear_mitc4_shell.h"
+#include "fesa/model/model_types.h"
+
+namespace {
+
+using Coordinates = std::array<double, 3>;
+
+fesa::Node Node(const std::int64_t label, const Coordinates coordinates) {
+  return {{"Shell-Instance", label, std::to_string(label)},
+          coordinates,
+          {"nonlinear-mitc4-tangent.inp", static_cast<std::size_t>(label)}};
+}
+
+std::array<fesa::Node, 4> Nodes() {
+  return {Node(1, {-1.0, -1.0, 0.0}), Node(2, {1.0, -1.0, 0.0}),
+          Node(3, {1.0, 1.0, 0.0}), Node(4, {-1.0, 1.0, 0.0})};
+}
+
+std::array<const fesa::Node*, 4> NodePointers(
+    const std::array<fesa::Node, 4>& nodes) {
+  return {&nodes[0], &nodes[1], &nodes[2], &nodes[3]};
+}
+
+std::array<fesa::ShellNodeInitialFrame, 4> Frames() {
+  return {{{0U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
+           {1U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
+           {2U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}},
+           {3U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}}}};
+}
+
+fesa::NonlinearMitc4Shell Shell(const std::array<fesa::Node, 4>& nodes) {
+  const fesa::ShellSection section{
+      "Section", 0.2, 0U, {"nonlinear-mitc4-tangent.inp", 20U}};
+  const fesa::LinearElasticMaterial material{
+      "Material", 120.0, 0.25, {"nonlinear-mitc4-tangent.inp", 21U}};
+  const auto shell = fesa::NonlinearMitc4Shell::Create(
+      NodePointers(nodes), Frames(), section, material);
+  if (!shell.HasValue()) {
+    throw std::runtime_error{shell.GetStatus().Diagnostics()[0U].code};
+  }
+  return shell.Value();
+}
+
+fesa::NonlinearMitc4TrialDofs Trial() {
+  fesa::NonlinearMitc4TrialDofs trial{};
+  trial.displacement[0U] = {-0.02, 0.01, 0.03};
+  trial.displacement[1U] = {0.01, -0.03, -0.01};
+  trial.displacement[2U] = {0.03, 0.02, 0.04};
+  trial.displacement[3U] = {-0.01, 0.04, -0.02};
+  trial.rotation_vector[0U] = {0.03, -0.02, 0.01};
+  trial.rotation_vector[1U] = {-0.01, 0.04, -0.02};
+  trial.rotation_vector[2U] = {0.02, 0.01, 0.03};
+  trial.rotation_vector[3U] = {-0.03, 0.02, -0.01};
+  return trial;
+}
+
+fesa::NonlinearMitc4TrialDofs AddScaled(
+    const fesa::NonlinearMitc4TrialDofs& trial,
+    const fesa::NonlinearMitc4TrialDofs& direction, const double scale) {
+  fesa::NonlinearMitc4TrialDofs result = trial;
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    result.displacement[node] =
+        result.displacement[node] + scale * direction.displacement[node];
+    result.rotation_vector[node] =
+        result.rotation_vector[node] + scale * direction.rotation_vector[node];
+  }
+  return result;
+}
+
+fesa::NonlinearMitc4TrialDofs Direction() {
+  fesa::NonlinearMitc4TrialDofs direction{};
+  direction.displacement[0U] = {0.5, -0.2, 0.1};
+  direction.displacement[1U] = {-0.1, 0.3, -0.4};
+  direction.displacement[2U] = {0.2, 0.1, 0.5};
+  direction.displacement[3U] = {-0.3, -0.4, -0.2};
+  direction.rotation_vector[0U] = {0.2, -0.1, 0.3};
+  direction.rotation_vector[1U] = {-0.3, 0.2, -0.1};
+  direction.rotation_vector[2U] = {0.1, 0.3, -0.2};
+  direction.rotation_vector[3U] = {-0.2, -0.3, 0.1};
+  return direction;
+}
+
+fesa::Vector ToVector(const fesa::NonlinearMitc4TrialDofs& trial) {
+  fesa::Vector result{24U};
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    for (std::size_t component = 0U; component < 3U; ++component) {
+      result[6U * node + component] = trial.displacement[node][component];
+      result[6U * node + 3U + component] =
+          trial.rotation_vector[node][component];
+    }
+  }
+  return result;
+}
+
+double ReferenceLength(const std::array<fesa::Node, 4>& nodes) {
+  double length = 0.0;
+  for (std::size_t first = 0U; first < nodes.size(); ++first) {
+    for (std::size_t second = first + 1U; second < nodes.size(); ++second) {
+      length = (std::max)(length, (fesa::Vector3{nodes[first].coordinates} -
+                                   fesa::Vector3{nodes[second].coordinates})
+                                      .Norm());
+    }
+  }
+  return length;
+}
+
+fesa::NonlinearMitc4TrialDofs WorkConjugateDirection(
+    const std::array<fesa::Node, 4>& nodes) {
+  auto direction = Direction();
+  double maximum_component = 0.0;
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    for (std::size_t component = 0U; component < 3U; ++component) {
+      maximum_component =
+          (std::max)(maximum_component,
+                     std::abs(direction.displacement[node][component]));
+      maximum_component =
+          (std::max)(maximum_component,
+                     std::abs(direction.rotation_vector[node][component]));
+    }
+  }
+  const double length = ReferenceLength(nodes);
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    direction.displacement[node] =
+        length * direction.displacement[node] / maximum_component;
+    direction.rotation_vector[node] =
+        direction.rotation_vector[node] / maximum_component;
+  }
+  return direction;
+}
+
+fesa::Vector WorkConjugate(const fesa::Vector& vector, const double length) {
+  fesa::Vector result = vector;
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    for (std::size_t component = 0U; component < 3U; ++component) {
+      result[6U * node + component] *= length;
+    }
+  }
+  return result;
+}
+
+double ChannelNorm(const fesa::Vector& vector, const bool translation) {
+  double sum = 0.0;
+  const std::size_t offset = translation ? 0U : 3U;
+  for (std::size_t node = 0U; node < 4U; ++node) {
+    for (std::size_t component = 0U; component < 3U; ++component) {
+      const double value = vector[6U * node + offset + component];
+      sum += value * value;
+    }
+  }
+  return std::sqrt(sum);
+}
+
+double FrobeniusNorm(const fesa::Matrix& matrix) {
+  double sum = 0.0;
+  for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
+    for (std::size_t column = 0U; column < matrix.Columns(); ++column) {
+      sum += matrix(row, column) * matrix(row, column);
+    }
+  }
+  return std::sqrt(sum);
+}
+
+fesa::Matrix WorkConjugateTangent(const fesa::Matrix& tangent,
+                                  const double length) {
+  fesa::Matrix result = tangent;
+  for (std::size_t row = 0U; row < tangent.Rows(); ++row) {
+    const double row_scale = row % 6U < 3U ? length : 1.0;
+    for (std::size_t column = 0U; column < tangent.Columns(); ++column) {
+      const double column_scale = column % 6U < 3U ? length : 1.0;
+      result(row, column) *= row_scale * column_scale;
+    }
+  }
+  return result;
+}
+
+void ExpectRelative(const double actual, const double expected) {
+  if (expected == 0.0) {
+    EXPECT_LE(std::abs(actual), 1.0e-12);
+    return;
+  }
+  EXPECT_LE(std::abs(actual - expected) / std::abs(expected), 1.0e-9);
+}
+
+// NSNR-ELEM-TAN-001
+TEST(NonlinearMitc4Residual, EvaluatesResidualTangentAndRecoveryAtOneTrial) {
+  const auto nodes = Nodes();
+  const auto shell = Shell(nodes);
+  const auto trial = Trial();
+  const auto recovery = shell.RecoverState(trial);
+  const auto result = fesa::EvaluateInternalResidualAndTangent(shell, trial);
+
+  ASSERT_TRUE(recovery.HasValue());
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_EQ(result.Value().internal_residual.Size(), 24U);
+  EXPECT_EQ(result.Value().internal_tangent.Rows(), 24U);
+  EXPECT_EQ(result.Value().internal_tangent.Columns(), 24U);
+  EXPECT_DOUBLE_EQ(result.Value().physical_strain_energy,
+                   result.Value().recovery.physical_strain_energy);
+  EXPECT_DOUBLE_EQ(result.Value().recovery.physical_strain_energy,
+                   recovery.Value().physical_strain_energy);
+  for (std::size_t index = 0U; index < 24U; ++index) {
+    EXPECT_DOUBLE_EQ(result.Value().internal_residual[index],
+                     result.Value().internal_force[index]);
+  }
+  for (std::size_t point = 0U;
+       point < result.Value().recovery.gauss_points.size(); ++point) {
+    EXPECT_EQ(result.Value().recovery.gauss_points[point].generalized_strain,
+              recovery.Value().gauss_points[point].generalized_strain);
+    EXPECT_EQ(result.Value().recovery.gauss_points[point].section_resultant,
+              recovery.Value().gauss_points[point].section_resultant);
+  }
+}
+
+// NSNR-ELEM-TAN-002 and NSNR-ELEM-TAN-003
+TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) {
+  const auto nodes = Nodes();
+  const auto shell = Shell(nodes);
+  const auto trial = Trial();
+  const auto direction = WorkConjugateDirection(nodes);
+  const double length = ReferenceLength(nodes);
+  const auto equilibrium =
+      fesa::EvaluateInternalResidualAndTangent(shell, trial);
+  ASSERT_TRUE(equilibrium.HasValue());
+
+  const fesa::Vector direction_vector = ToVector(direction);
+  const fesa::Vector analytical = WorkConjugate(
+      equilibrium.Value().internal_tangent.Multiply(direction_vector), length);
+  const std::array<double, 8> steps{1.0e-3, 3.0e-4, 1.0e-4, 3.0e-5,
+                                    1.0e-5, 3.0e-6, 1.0e-6, 3.0e-7};
+  std::array<double, 8> force_errors{};
+  std::array<double, 8> moment_errors{};
+  ASSERT_GT(ChannelNorm(analytical, true), 0.0);
+  ASSERT_GT(ChannelNorm(analytical, false), 0.0);
+  for (std::size_t index = 0U; index < steps.size(); ++index) {
+    const auto plus = fesa::EvaluateInternalResidualAndTangent(
+        shell, AddScaled(trial, direction, steps[index]));
+    const auto minus = fesa::EvaluateInternalResidualAndTangent(
+        shell, AddScaled(trial, direction, -steps[index]));
+    ASSERT_TRUE(plus.HasValue());
+    ASSERT_TRUE(minus.HasValue());
+    fesa::Vector derivative = plus.Value().internal_residual;
+    derivative.Axpy(-1.0, minus.Value().internal_residual);
+    derivative.Scale(0.5 / steps[index]);
+    derivative = WorkConjugate(derivative, length);
+    derivative.Axpy(-1.0, analytical);
+    force_errors[index] =
+        ChannelNorm(derivative, true) / ChannelNorm(analytical, true);
+    moment_errors[index] =
+        ChannelNorm(derivative, false) / ChannelNorm(analytical, false);
+  }
+
+  std::size_t force_order_count = 0U;
+  std::size_t moment_order_count = 0U;
+  for (std::size_t index = 1U; index < force_errors.size(); ++index) {
+    if (force_errors[index - 1U] > 0.0 && force_errors[index] > 0.0) {
+      const double order =
+          std::log(force_errors[index - 1U] / force_errors[index]) /
+          std::log(steps[index - 1U] / steps[index]);
+      if (order >= 1.5) {
+        ++force_order_count;
+      }
+    }
+    if (moment_errors[index - 1U] > 0.0 && moment_errors[index] > 0.0) {
+      const double order =
+          std::log(moment_errors[index - 1U] / moment_errors[index]) /
+          std::log(steps[index - 1U] / steps[index]);
+      if (order >= 1.5) {
+        ++moment_order_count;
+      }
+    }
+  }
+  EXPECT_GE(force_order_count, 2U);
+  EXPECT_GE(moment_order_count, 2U);
+  EXPECT_LE(force_errors[3U], 1.0e-6);
+  EXPECT_LE(moment_errors[3U], 1.0e-6);
+}
+
+// NSNR-ELEM-TAN-004 and NSNR-ELEM-TAN-005
+TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) {
+  const auto nodes = Nodes();
+  const auto equilibrium =
+      fesa::EvaluateInternalResidualAndTangent(Shell(nodes), {});
+  ASSERT_TRUE(equilibrium.HasValue());
+  const double length = ReferenceLength(nodes);
+
+  fesa::NonlinearMitc4TrialDofs rigid_translation{};
+  for (auto& displacement : rigid_translation.displacement) {
+    displacement = {0.3, -0.2, 0.4};
+  }
+  const fesa::Vector null_action =
+      WorkConjugate(equilibrium.Value().internal_tangent.Multiply(
+                        ToVector(rigid_translation)),
+                    length);
+  EXPECT_LE(null_action.Norm(), 1.0e-10);
+
+  const fesa::Matrix scaled_tangent =
+      WorkConjugateTangent(equilibrium.Value().internal_tangent, length);
+  fesa::Matrix antisymmetric = scaled_tangent;
+  for (std::size_t row = 0U; row < antisymmetric.Rows(); ++row) {
+    for (std::size_t column = 0U; column < antisymmetric.Columns(); ++column) {
+      antisymmetric(row, column) -= scaled_tangent(column, row);
+    }
+  }
+  ASSERT_GT(FrobeniusNorm(scaled_tangent), 0.0);
+  EXPECT_LE(FrobeniusNorm(antisymmetric) / FrobeniusNorm(scaled_tangent),
+            1.0e-12);
+}
+
+// NSNR-ELEM-TAN-006
+TEST(NonlinearMitc4Residual, DoesNotAcceptExternalLoadShares) {
+  const auto nodes = Nodes();
+  const auto first =
+      fesa::EvaluateInternalResidualAndTangent(Shell(nodes), Trial());
+  const auto second =
+      fesa::EvaluateInternalResidualAndTangent(Shell(nodes), Trial());
+  ASSERT_TRUE(first.HasValue());
+  ASSERT_TRUE(second.HasValue());
+  for (std::size_t index = 0U; index < 24U; ++index) {
+    EXPECT_DOUBLE_EQ(first.Value().internal_residual[index],
+                     second.Value().internal_residual[index]);
+  }
+}
+
+// NSNR-ELEM-TAN-007
+TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) {
+  const auto nodes = Nodes();
+  const fesa::Vector3 rotation{0.3, -0.2, 0.1};
+  const fesa::Matrix tensor = fesa::ExpSkew(rotation);
+  fesa::NonlinearMitc4TrialDofs trial{};
+  for (std::size_t node = 0U; node < nodes.size(); ++node) {
+    const fesa::Vector3 initial{nodes[node].coordinates};
+    const fesa::Vector3 rotated{
+        tensor(0U, 0U) * initial.X() + tensor(0U, 1U) * initial.Y() +
+            tensor(0U, 2U) * initial.Z(),
+        tensor(1U, 0U) * initial.X() + tensor(1U, 1U) * initial.Y() +
+            tensor(1U, 2U) * initial.Z(),
+        tensor(2U, 0U) * initial.X() + tensor(2U, 1U) * initial.Y() +
+            tensor(2U, 2U) * initial.Z()};
+    trial.displacement[node] =
+        rotated + fesa::Vector3{0.7, -0.4, 1.2} - initial;
+    trial.rotation_vector[node] = rotation;
+  }
+  const auto result =
+      fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial);
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_LE(result.Value().internal_residual.Norm(), 1.0e-10);
+  EXPECT_LE(result.Value().physical_strain_energy, 1.0e-12);
+}
+
+// NSNR-ELEM-TAN-008
+TEST(NonlinearMitc4Patch, RetainsPhysicalRecoveryIdentities) {
+  const auto nodes = Nodes();
+  const auto shell = Shell(nodes);
+  constexpr double kMembraneStrain = 1.0e-2;
+  const double green_strain =
+      kMembraneStrain + 0.5 * kMembraneStrain * kMembraneStrain;
+  constexpr double kYoungsModulus = 120.0;
+  constexpr double kPoissonRatio = 0.25;
+  constexpr double kThickness = 0.2;
+  const double plane_stress =
+      kYoungsModulus / (1.0 - kPoissonRatio * kPoissonRatio);
+  const double s11 = plane_stress * green_strain;
+  const double s22 = plane_stress * kPoissonRatio * green_strain;
+  fesa::NonlinearMitc4TrialDofs trial{};
+  for (std::size_t node = 0U; node < nodes.size(); ++node) {
+    trial.displacement[node] = {kMembraneStrain * nodes[node].coordinates[0U],
+                                0.0, 0.0};
+  }
+  const auto result = fesa::EvaluateInternalResidualAndTangent(shell, trial);
+  ASSERT_TRUE(result.HasValue());
+  for (const auto& point : result.Value().recovery.gauss_points) {
+    ExpectRelative(point.generalized_strain[0U], green_strain);
+    for (std::size_t component = 1U; component < 8U; ++component) {
+      ExpectRelative(point.generalized_strain[component], 0.0);
+    }
+    ExpectRelative(point.section_resultant[0U], kThickness * s11);
+    ExpectRelative(point.section_resultant[1U], kThickness * s22);
+    for (std::size_t component = 2U; component < 8U; ++component) {
+      ExpectRelative(point.section_resultant[component], 0.0);
+    }
+    for (const auto& stress : point.stress) {
+      ExpectRelative(stress.components[0U], s11);
+      ExpectRelative(stress.components[1U], s22);
+      ExpectRelative(stress.components[2U], 0.0);
+    }
+  }
+}
+
+// NSNR-ELEM-TAN-009
+TEST(NonlinearMitc4Residual, RejectsNonfiniteTrialData) {
+  const auto nodes = Nodes();
+  auto trial = Trial();
+  trial.rotation_vector[2U] = {std::numeric_limits<double>::quiet_NaN(), 0.0,
+                               0.0};
+  const auto result =
+      fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial);
+  ASSERT_FALSE(result.HasValue());
+  EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code,
+            "invalid-nonlinear-shell-trial");
+
+  auto overflow = Trial();
+  overflow.displacement[1U] = {1.0e300, 0.0, 0.0};
+  const auto overflow_result =
+      fesa::EvaluateInternalResidualAndTangent(Shell(nodes), overflow);
+  EXPECT_FALSE(overflow_result.HasValue());
+}
+
+}  // namespace