feat: adopt lumped mass and reduced-coordinate shell HHT

Integrate Hex8 mass verification and fixed-null-basis shell dynamics with numerical and physics handoff evidence. Preserve interrupted reference and full-validation status; full integration validation is deferred until after performance optimization.
2026-10-13 16:57:19 +09:00
parent ede137c6b4
commit eaa5cbeb21
110 changed files with 15274 additions and 971 deletions
@@ -219,7 +219,7 @@ Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
  if (!density_)
    return Result<ElementStiffnessContribution>::Failure(
        Failure("Hex8 mass requires reference density."));
-  Matrix mass{24U, 24U};
+  Matrix consistent{24U, 24U};
  for (const auto& point : points_)
    for (std::size_t a = 0; a < 8; ++a)
      for (std::size_t b = 0; b < 8; ++b) {
@@ -229,15 +229,18 @@ Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
          return Result<ElementStiffnessContribution>::Failure(
              Failure("Hex8 reference mass is nonfinite."));
        for (std::size_t i = 0; i < 3; ++i)
-          mass(3U * a + i, 3U * b + i) += value;
+          consistent(3U * a + i, 3U * b + i) += value;
      }
+  Matrix mass{24U, 24U};
  for (std::size_t i = 0; i < 24; ++i)
-    for (std::size_t j = 0; j < 24; ++j)
-      if (!std::isfinite(mass(i, j)))
+    for (std::size_t j = 0; j < 24; ++j) {
+      if (!std::isfinite(consistent(i, j)))
        return Result<ElementStiffnessContribution>::Failure(
            Failure("Hex8 accumulated reference mass is nonfinite."));
+      mass(i, i) += consistent(i, j);
+    }
  for (std::size_t i = 0; i < 24; ++i)
-    if (mass(i, i) <= 0.0)
+    if (!std::isfinite(mass(i, i)) || mass(i, i) <= 0.0)
      return Result<ElementStiffnessContribution>::Failure(Failure(
          "Hex8 reference mass has an unrepresentable positive diagonal."));
  return Result<ElementStiffnessContribution>::Success(
@@ -19,6 +19,49 @@ const std::array<double, 3> kGauss{-std::sqrt(3.0 / 5.0), 0.0,
                                   std::sqrt(3.0 / 5.0)};
 constexpr std::array<double, 3> kWeights{5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0};

+class InertiaSum {
+ public:
+  void Add(const double term) {
+    const double next = sum_ + term;
+    correction_ += std::abs(sum_) >= std::abs(term)
+                       ? (sum_ - next) + term
+                       : (term - next) + sum_;
+    sum_ = next;
+  }
+  double Value() const { return sum_ + correction_; }
+
+ private:
+  double sum_ = 0.0;
+  double correction_ = 0.0;
+};
+
+Vector3 WeightedSum(const std::array<double, kNodes>& weights,
+                    const std::array<Vector3, kNodes>& values) {
+  Vector3 result{};
+  for (std::size_t node = 0U; node < kNodes; ++node)
+    result = result + weights[node] * values[node];
+  return result;
+}
+
+// A derivative partition has no absolute-position contribution. Anchoring
+// before reduction preserves that identity at large global offsets and keeps
+// the reference/current covariants in the same canonical form.
+Vector3 DerivativeSum(const std::array<double, kNodes>& derivatives,
+                      const std::array<Vector3, kNodes>& values) {
+  std::array<Vector3, kNodes> relative{};
+  for (std::size_t node = 0U; node < kNodes; ++node)
+    relative[node] = values[node] - values[0U];
+  return WeightedSum(derivatives, relative);
+}
+
+Vector3 InterpolatedSum(const std::array<double, kNodes>& weights,
+                        const std::array<Vector3, kNodes>& values) {
+  std::array<Vector3, kNodes> relative{};
+  for (std::size_t node = 0U; node < kNodes; ++node)
+    relative[node] = values[node] - values[0U];
+  return values[0U] + WeightedSum(weights, relative);
+}
+
 template <class T>
 Result<T> Failure(const std::string& code, const SourceLocation& location,
                  const std::string& identity, std::string message) {
@@ -723,131 +766,232 @@ Result<NonlinearMitc4InertiaEvaluation> NonlinearMitc4Inertia::Evaluate(
                                         {},
                                         {},
                                         {}};
+  using LocalDerivative = std::array<std::array<std::array<double, 6U>, 6U>, 3U>;
+  std::array<LocalDerivative, kNodes> mass_derivative{};
+  // Integrate acceleration/velocity coefficients before their contraction.
+  // This preserves rotary contributions until odd-in-z translation cancels.
+  std::array<Vector3, kDofs> momentum_coefficients{};
+  std::array<Vector3, kDofs> angular_coefficients{};
+  Vector3 convective_momentum{};
+  Vector3 convective_angular{};
  for (std::size_t xi = 0U; xi < kGauss.size(); ++xi) {
    for (std::size_t eta = 0U; eta < kGauss.size(); ++eta) {
      const Mitc4ShapeFunctions shape =
          Mitc4Shell::ShapeFunctions(kGauss[xi], kGauss[eta]);
+      std::array<Vector3, kNodes> displacements{};
+      std::array<Vector3, kNodes> current_coordinates{};
+      std::array<Vector3, kNodes> current_directors{};
+      for (std::size_t node = 0U; node < kNodes; ++node) {
+        displacements[node] = trial.displacement[node];
+        current_coordinates[node] = coordinates_[node] + displacements[node];
+        current_directors[node] = rotations[node].director;
+      }
+      const Vector3 reference_xi_anchor =
+          DerivativeSum(shape.xi_derivatives, coordinates_);
+      const Vector3 reference_eta_anchor =
+          DerivativeSum(shape.eta_derivatives, coordinates_);
+      const Vector3 reference_xi_director =
+          DerivativeSum(shape.xi_derivatives, initial_directors_);
+      const Vector3 reference_eta_director =
+          DerivativeSum(shape.eta_derivatives, initial_directors_);
+      const Vector3 reference_zeta =
+          0.5 * thickness_ * InterpolatedSum(shape.values, initial_directors_);
+      const Vector3 current_xi_anchor =
+          reference_xi_anchor + DerivativeSum(shape.xi_derivatives, displacements);
+      const Vector3 current_eta_anchor =
+          reference_eta_anchor + DerivativeSum(shape.eta_derivatives, displacements);
+      const Vector3 current_xi_director =
+          DerivativeSum(shape.xi_derivatives, current_directors);
+      const Vector3 current_eta_director =
+          DerivativeSum(shape.eta_derivatives, current_directors);
+      const Vector3 current_zeta =
+          0.5 * thickness_ * InterpolatedSum(shape.values, current_directors);
      for (std::size_t zeta = 0U; zeta < kGauss.size(); ++zeta) {
-        const double z = 0.5 * thickness_ * kGauss[zeta];
-        const double weight = kWeights[xi] * kWeights[eta] * kWeights[zeta];
-        Vector3 ref_xi{};
-        Vector3 ref_eta{};
-        Vector3 ref_zeta{};
-        Vector3 current_xi{};
-        Vector3 current_eta{};
-        Vector3 current_zeta{};
-        Vector3 position{};
-        Vector3 velocity_point{};
-        Vector3 acceleration_point{};
-        Vector3 convective_point{};
-        std::array<Vector3, kDofs> h{};
-        std::array<Vector3, kDofs> velocity_derivative{};
-        std::array<Vector3, kDofs> position_derivative{};
-        for (std::size_t node = 0U; node < kNodes; ++node) {
-          const std::size_t base = 6U * node;
-          const auto& rotation = rotations[node];
-          const Vector3 reference =
-              coordinates_[node] + z * initial_directors_[node];
-          const Vector3 current = coordinates_[node] +
-                                  trial.displacement[node] +
-                                  z * rotation.director;
-          ref_xi = ref_xi + shape.xi_derivatives[node] * reference;
-          ref_eta = ref_eta + shape.eta_derivatives[node] * reference;
-          ref_zeta = ref_zeta + 0.5 * thickness_ * shape.values[node] *
-                                    initial_directors_[node];
-          current_xi = current_xi + shape.xi_derivatives[node] * current;
-          current_eta = current_eta + shape.eta_derivatives[node] * current;
-          current_zeta = current_zeta + 0.5 * thickness_ * shape.values[node] *
-                                            rotation.director;
-          position = position + shape.values[node] * current;
-          for (std::size_t component = 0U; component < 3U; ++component) {
-            std::array<double, 3> translation{};
-            translation[component] = shape.values[node];
-            h[base + component] = Vector3{translation};
-            h[base + 3U + component] =
-                z * shape.values[node] * rotation.director_first[component];
-            for (std::size_t other = 0U; other < 3U; ++other) {
-              const Vector3 second = z * shape.values[node] *
-                                     rotation.director_second[component][other];
-              convective_point =
-                  convective_point + second * velocity[base + 3U + component] *
-                                         velocity[base + 3U + other];
-              velocity_derivative[base + 3U + component] =
-                  velocity_derivative[base + 3U + component] +
-                  2.0 * second * velocity[base + 3U + other];
-              position_derivative[base + 3U + component] =
-                  position_derivative[base + 3U + component] +
-                  second * acceleration[base + 3U + other];
-              for (std::size_t third = 0U; third < 3U; ++third) {
-                position_derivative[base + 3U + component] =
-                    position_derivative[base + 3U + component] +
-                    z * shape.values[node] *
-                        rotation.director_third[component][other][third] *
-                        velocity[base + 3U + other] *
-                        velocity[base + 3U + third];
+        // Use the nodal kinetic measure rho * N_i * dV0.
+        for (std::size_t mass_node = 0U; mass_node < kNodes; ++mass_node) {
+          std::array<double, kNodes> kinetic_shape{};
+          kinetic_shape[mass_node] = 1.0;
+          const double z = 0.5 * thickness_ * kGauss[zeta];
+          const double weight = kWeights[xi] * kWeights[eta] * kWeights[zeta];
+          const Vector3 ref_xi = reference_xi_anchor + z * reference_xi_director;
+          const Vector3 ref_eta = reference_eta_anchor + z * reference_eta_director;
+          const Vector3 current_xi = current_xi_anchor + z * current_xi_director;
+          const Vector3 current_eta = current_eta_anchor + z * current_eta_director;
+          Vector3 position{};
+          Vector3 velocity_point{};
+          Vector3 acceleration_point{};
+          Vector3 convective_point{};
+          std::array<Vector3, kDofs> h{};
+          std::array<Vector3, kDofs> velocity_derivative{};
+          std::array<Vector3, kDofs> position_derivative{};
+          for (std::size_t node = 0U; node < kNodes; ++node) {
+            const std::size_t base = 6U * node;
+            const auto& rotation = rotations[node];
+            const Vector3 current =
+                current_coordinates[node] + z * current_directors[node];
+            position = position + kinetic_shape[node] * current;
+            for (std::size_t component = 0U; component < 3U; ++component) {
+              std::array<double, 3> translation{};
+              translation[component] = kinetic_shape[node];
+              h[base + component] = Vector3{translation};
+              h[base + 3U + component] =
+                  z * kinetic_shape[node] * rotation.director_first[component];
+              for (std::size_t other = 0U; other < 3U; ++other) {
+                const Vector3 second =
+                    z * kinetic_shape[node] *
+                    rotation.director_second[component][other];
+                convective_point = convective_point +
+                                   second * velocity[base + 3U + component] *
+                                       velocity[base + 3U + other];
+                velocity_derivative[base + 3U + component] =
+                    velocity_derivative[base + 3U + component] +
+                    2.0 * second * velocity[base + 3U + other];
+                for (std::size_t third = 0U; third < 3U; ++third) {
+                  position_derivative[base + 3U + component] =
+                      position_derivative[base + 3U + component] +
+                      z * kinetic_shape[node] *
+                          rotation.director_third[component][other][third] *
+                          velocity[base + 3U + other] *
+                          velocity[base + 3U + third];
+                }
              }
            }
          }
-        }
-        const double jacobian = Determinant(ref_xi, ref_eta, ref_zeta);
-        const double current_jacobian =
-            Determinant(current_xi, current_eta, current_zeta);
-        if (!std::isfinite(jacobian) || !(jacobian > 0.0) ||
-            !std::isfinite(current_jacobian) || !(current_jacobian > 0.0)) {
-          return Failure<NonlinearMitc4InertiaEvaluation>(
-              "invalid-nonlinear-inertia-geometry", source_location_, identity_,
-              "Inertia requires positive reference and current mass-point "
-              "Jacobians.");
-        }
-        for (std::size_t row = 0U; row < kDofs; ++row) {
-          velocity_point = velocity_point + h[row] * velocity[row];
-          acceleration_point = acceleration_point + h[row] * acceleration[row];
-        }
-        acceleration_point = acceleration_point + convective_point;
-        const double scale = density_ * weight * jacobian;
-        if (!std::isfinite(scale) || !(scale > 0.0)) {
-          return Failure<NonlinearMitc4InertiaEvaluation>(
-              "invalid-nonlinear-inertia-positive-mass", source_location_,
-              identity_,
-              "Mass-point integration weight must be finite and positive.");
-        }
-        result.kinetic_energy +=
-            0.5 * scale * velocity_point.Dot(velocity_point);
-        result.momentum = result.momentum + scale * velocity_point;
-        result.momentum_rate =
-            result.momentum_rate + scale * acceleration_point;
-        result.angular_momentum =
-            result.angular_momentum + scale * position.Cross(velocity_point);
-        result.angular_momentum_rate =
-            result.angular_momentum_rate +
-            scale * position.Cross(acceleration_point);
-        result.momentum_rate_scale += scale * acceleration_point.Norm();
-        result.angular_momentum_rate_scale +=
-            scale * position.Cross(acceleration_point).Norm();
-        for (std::size_t row = 0U; row < kDofs; ++row) {
-          result.force[row] += scale * h[row].Dot(acceleration_point);
-          // Keep c independent of a; subtracting M*a destroys small c values.
-          result.convective_force[row] += scale * h[row].Dot(convective_point);
-          for (std::size_t column = 0U; column < kDofs; ++column) {
-            result.mass(row, column) += scale * h[row].Dot(h[column]);
-            result.velocity_tangent(row, column) +=
-                scale * h[row].Dot(velocity_derivative[column]);
-            double position_entry = h[row].Dot(position_derivative[column]);
-            if (row / 6U == column / 6U && row % 6U >= 3U &&
-                column % 6U >= 3U) {
-              const std::size_t node = row / 6U;
-              position_entry +=
-                  z * shape.values[node] *
-                  rotations[node]
-                      .director_second[row % 6U - 3U][column % 6U - 3U]
-                      .Dot(acceleration_point);
+          const double jacobian =
+              Determinant(ref_xi, ref_eta, reference_zeta);
+          const double current_jacobian =
+              Determinant(current_xi, current_eta, current_zeta);
+          if (!std::isfinite(jacobian) || !(jacobian > 0.0) ||
+              !std::isfinite(current_jacobian) || !(current_jacobian > 0.0)) {
+            return Failure<NonlinearMitc4InertiaEvaluation>(
+                "invalid-nonlinear-inertia-geometry", source_location_,
+                identity_,
+                "Inertia requires positive reference and current mass-point "
+                "Jacobians.");
+          }
+          for (std::size_t row = 6U * mass_node; row < 6U * (mass_node + 1U);
+               ++row) {
+            velocity_point = velocity_point + h[row] * velocity[row];
+            acceleration_point =
+                acceleration_point + h[row] * acceleration[row];
+          }
+          acceleration_point = acceleration_point + convective_point;
+          const double scale =
+              density_ * weight * jacobian * shape.values[mass_node];
+          if (!std::isfinite(scale) || !(scale > 0.0)) {
+            return Failure<NonlinearMitc4InertiaEvaluation>(
+                "invalid-nonlinear-inertia-positive-mass", source_location_,
+                identity_,
+                "Mass-point integration weight must be finite and positive.");
+          }
+          result.kinetic_energy +=
+              0.5 * scale * velocity_point.Dot(velocity_point);
+          for (std::size_t field = 6U * mass_node;
+               field < 6U * (mass_node + 1U); ++field) {
+            momentum_coefficients[field] = momentum_coefficients[field] + scale * h[field];
+            angular_coefficients[field] =
+                angular_coefficients[field] + scale * position.Cross(h[field]);
+          }
+          convective_momentum =
+              convective_momentum + scale * convective_point;
+          convective_angular =
+              convective_angular + scale * position.Cross(convective_point);
+          result.momentum_rate_scale += scale * acceleration_point.Norm();
+          result.angular_momentum_rate_scale +=
+              scale * position.Cross(acceleration_point).Norm();
+          for (std::size_t row = 6U * mass_node; row < 6U * (mass_node + 1U);
+               ++row) {
+            // Keep c independent of a; subtracting M*a destroys small c values.
+            result.convective_force[row] +=
+                scale * h[row].Dot(convective_point);
+            for (std::size_t column = 6U * mass_node;
+                 column < 6U * (mass_node + 1U); ++column) {
+              result.mass(row, column) += scale * h[row].Dot(h[column]);
+              result.velocity_tangent(row, column) +=
+                  scale * h[row].Dot(velocity_derivative[column]);
+              double position_entry = h[row].Dot(position_derivative[column]);
+              if (row / 6U == column / 6U && row % 6U >= 3U &&
+                  column % 6U >= 3U) {
+                const std::size_t node = row / 6U;
+                position_entry +=
+                    z * kinetic_shape[node] *
+                    rotations[node]
+                        .director_second[row % 6U - 3U][column % 6U - 3U]
+                        .Dot(convective_point);
+              }
+              result.position_tangent(row, column) += scale * position_entry;
+              if (column % 6U >= 3U) {
+                const auto direction = column % 6U - 3U;
+                const Vector3 dh_row =
+                    row % 6U < 3U
+                        ? Vector3{}
+                        : z * rotations[mass_node]
+                                  .director_second[row % 6U - 3U][direction];
+                for (std::size_t local = 0U; local < 6U; ++local) {
+                  const auto field = 6U * mass_node + local;
+                  const Vector3 dh_field =
+                      local < 3U
+                          ? Vector3{}
+                          : z * rotations[mass_node]
+                                    .director_second[local - 3U][direction];
+                  mass_derivative[mass_node][direction][row % 6U][local] +=
+                      scale * (dh_row.Dot(h[field]) + h[row].Dot(dh_field));
+                }
+              }
            }
-            result.position_tangent(row, column) += scale * position_entry;
          }
        }
      }
    }
  }
+  // Integrate M and c before contracting acceleration. Mixing a large
+  // translation and tiny rotary acceleration at a mass point can lose the
+  // latter before the odd-in-z translational moment cancels in integration.
+  for (std::size_t row = 0U; row < kDofs; ++row) {
+    InertiaSum sum;
+    for (std::size_t column = 0U; column < kDofs; ++column)
+      sum.Add(result.mass(row, column) * acceleration[column]);
+    sum.Add(result.convective_force[row]);
+    result.force[row] = sum.Value();
+  }
+  // DqI = (DqM)*a + Dqc.  The integrated tangent above contains only Dqc;
+  // its acceleration contribution uses the independently integrated DqM.
+  for (std::size_t node = 0U; node < kNodes; ++node)
+    for (std::size_t direction = 0U; direction < 3U; ++direction)
+      for (std::size_t row = 0U; row < 6U; ++row) {
+        InertiaSum sum;
+        for (std::size_t field = 0U; field < 6U; ++field)
+          sum.Add(mass_derivative[node][direction][row][field] *
+                  acceleration[6U * node + field]);
+        sum.Add(result.position_tangent(6U * node + row,
+                                         6U * node + 3U + direction));
+        result.position_tangent(6U * node + row,
+                                  6U * node + 3U + direction) = sum.Value();
+      }
+  std::array<double, 3U> momentum{};
+  std::array<double, 3U> angular{};
+  std::array<double, 3U> momentum_rate{};
+  std::array<double, 3U> angular_rate{};
+  for (std::size_t component = 0U; component < 3U; ++component) {
+    const auto contract = [&](const std::array<Vector3, kDofs>& coefficients,
+                              const Vector& channel, const double extra) {
+      InertiaSum sum;
+      for (std::size_t field = 0U; field < kDofs; ++field)
+        sum.Add(coefficients[field][component] * channel[field]);
+      sum.Add(extra);
+      return sum.Value();
+    };
+    momentum[component] = contract(momentum_coefficients, velocity, 0.0);
+    angular[component] = contract(angular_coefficients, velocity, 0.0);
+    momentum_rate[component] = contract(momentum_coefficients, acceleration,
+                                        convective_momentum[component]);
+    angular_rate[component] = contract(angular_coefficients, acceleration,
+                                       convective_angular[component]);
+  }
+  result.momentum = Vector3{momentum};
+  result.angular_momentum = Vector3{angular};
+  result.momentum_rate = Vector3{momentum_rate};
+  result.angular_momentum_rate = Vector3{angular_rate};
  if (!IsFinite(result.mass) || !IsFinite(result.position_tangent) ||
      !IsFinite(result.velocity_tangent) || !IsFinite(result.force) ||
      !IsFinite(result.convective_force) ||
@@ -26,6 +26,7 @@ using EngineeringStress = std::array<double, 5>;
 struct Geometry {
  std::array<Vector3, 3> reference_covariant;
  std::array<Vector3, 3> current_covariant;
+  std::array<Vector3, 3> covariant_increment;
  std::array<Vector3, 3> reciprocal;
  Mitc4LocalFrame local_frame;
  double reference_jacobian;
@@ -145,6 +146,31 @@ bool MakeLocalFrame(const Vector3& midsurface_xi, const Vector3& midsurface_eta,
  return true;
 }

+/// @brief Forms (R(phi)-I)D0 without adding and subtracting the baseline.
+Vector3 DirectorIncrement(const Vector3& phi, const Vector3& initial) {
+  const double squared_norm = phi.Dot(phi);
+  double first = 1.0;
+  double second = 0.5;
+  if (squared_norm < 0.25) {
+    double power = 1.0;
+    double first_coefficient = 1.0;
+    double second_coefficient = 0.5;
+    for (unsigned int k = 1U; k <= 9U; ++k) {
+      power *= squared_norm;
+      first_coefficient /= -static_cast<double>((2U * k) * (2U * k + 1U));
+      second_coefficient /= -static_cast<double>((2U * k + 1U) * (2U * k + 2U));
+      first += first_coefficient * power;
+      second += second_coefficient * power;
+    }
+  } else {
+    const double norm = std::sqrt(squared_norm);
+    first = std::sin(norm) / norm;
+    second = (1.0 - std::cos(norm)) / squared_norm;
+  }
+  const Vector3 crossed = phi.Cross(initial);
+  return first * crossed + second * phi.Cross(crossed);
+}
+
 /// @brief Evaluates fixed-reference and trial covariant geometry together.
 bool EvaluateGeometry(const std::array<Vector3, kNodeCount>& coordinates,
                      const std::array<Vector3, kNodeCount>& initial_directors,
@@ -169,6 +195,18 @@ bool EvaluateGeometry(const std::array<Vector3, kNodeCount>& coordinates,
  const Vector3 director_eta = DerivativeSum(shape.eta_derivatives, directors);
  const Vector3 director = WeightedSum(shape.values, directors);
  const double half_thickness = 0.5 * thickness;
+  std::array<Vector3, kNodeCount> director_increment{};
+  for (std::size_t node = 0U; node < kNodeCount; ++node)
+    director_increment[node] =
+        DirectorIncrement(trial.rotation_vector[node], initial_directors[node]);
+  geometry.covariant_increment[0U] =
+      displacement_xi + half_thickness * zeta *
+                            DerivativeSum(shape.xi_derivatives, director_increment);
+  geometry.covariant_increment[1U] =
+      displacement_eta + half_thickness * zeta *
+                             DerivativeSum(shape.eta_derivatives, director_increment);
+  geometry.covariant_increment[2U] =
+      half_thickness * WeightedSum(shape.values, director_increment);

  geometry.reference_covariant[0U] =
      reference_xi + half_thickness * zeta * initial_director_xi;
@@ -215,10 +253,13 @@ 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[row][column] =
+          0.5 * (geometry.reference_covariant[row].Dot(
+                     geometry.covariant_increment[column]) +
+                 geometry.covariant_increment[row].Dot(
+                     geometry.reference_covariant[column]) +
+                 geometry.covariant_increment[row].Dot(
+                     geometry.covariant_increment[column]));
    }
  }
  strain[2U][2U] = 0.0;
@@ -278,8 +319,11 @@ EngineeringStress Stress(const EngineeringStrain& strain, const double youngs,
          shear * strain[4U]};
 }

+/// @brief Distinguishes invalid geometry from arithmetic recovery failure.
+enum class PointStatus { kValid, kInvalidGeometry, kNonfiniteRecovery };
+
 /// @brief Evaluates direct strain at one mandatory edge-midpoint tying point.
-bool EvaluateTyingSample(
+PointStatus EvaluateTyingSample(
    const std::array<Vector3, kNodeCount>& coordinates,
    const std::array<Vector3, kNodeCount>& initial_directors,
    const std::array<Vector3, kNodeCount>& directors,
@@ -288,14 +332,14 @@ bool EvaluateTyingSample(
  Geometry geometry{};
  if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
                        thickness, xi, eta, 0.0, geometry)) {
-    return false;
+    return PointStatus::kInvalidGeometry;
  }
  strain = DirectStrain(geometry);
-  return IsFinite(strain);
+  return IsFinite(strain) ? PointStatus::kValid : PointStatus::kNonfiniteRecovery;
 }

 /// @brief Computes one MITC-projected local strain and material stress point.
-bool EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
+PointStatus EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
                   const std::array<Vector3, kNodeCount>& initial_directors,
                   const std::array<Vector3, kNodeCount>& directors,
                   const NonlinearMitc4TrialDofs& trial,
@@ -305,7 +349,7 @@ bool EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
                   const double zeta, PointEvaluation& point) {
  if (!EvaluateGeometry(coordinates, initial_directors, directors, trial,
                        thickness, xi, eta, zeta, point.geometry)) {
-    return false;
+    return PointStatus::kInvalidGeometry;
  }
  Tensor covariant = DirectStrain(point.geometry);
  const Mitc4TyingWeights weights = Mitc4Shell::TyingWeights(xi, eta);
@@ -320,7 +364,9 @@ bool EvaluatePoint(const std::array<Vector3, kNodeCount>& coordinates,
  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);
+                 IsFinite(point.stress)
+             ? PointStatus::kValid
+             : PointStatus::kNonfiniteRecovery;
 }

 /// @brief Tests the fixed-order physical recovery inventory for finite values.
@@ -787,7 +833,8 @@ JetMatrix3 Inverse(const JetMatrix3& matrix) {
 }

 /// @brief Evaluates Rodrigues' map with an analytic near-zero series branch.
-JetMatrix3 ExpSkew(const JetVector3& rotation) {
+JetMatrix3 ExpSkew(const JetVector3& rotation, const bool increment_only = false) {
+  const JetMatrix3 identity = increment_only ? JetMatrix3{} : Identity();
  const Jet squared_norm = Dot(rotation, rotation);
  Jet sine_over_norm{};
  Jet one_minus_cosine_over_norm_squared{};
@@ -807,12 +854,12 @@ JetMatrix3 ExpSkew(const JetVector3& rotation) {
  }
  const JetMatrix3 skew = Skew(rotation);
  JetMatrix3 result =
-      Add(Add(Identity(), Scale(sine_over_norm, skew)),
+      Add(Add(identity, Scale(sine_over_norm, skew)),
          Scale(one_minus_cosine_over_norm_squared, Multiply(skew, skew)));
  if (squared_norm.value < 0.25) {
-    // Third derivatives amplify cancellation in 1-cos(theta). Preserve the
-    // established static coefficients and use a stable series for the new
-    // directional derivatives of the complete rotation tensor.
+    // Near the original small-angle boundary, cancellation in 1-cos(theta)
+    // also corrupts the physical-force gradient and Hessian. Use the existing
+    // stable series consistently for the complete Jet and rotation tensor.
    Jet sine_series = Constant(1.0);
    Jet cosine_series = Constant(0.5);
    Jet power = Constant(1.0);
@@ -825,18 +872,9 @@ JetMatrix3 ExpSkew(const JetVector3& rotation) {
      sine_series = sine_series + sine_coefficient * power;
      cosine_series = cosine_series + cosine_coefficient * power;
    }
-    const JetMatrix3 stable = Add(Add(Identity(), Scale(sine_series, skew)),
+    const JetMatrix3 stable = Add(Add(identity, Scale(sine_series, skew)),
                                  Scale(cosine_series, Multiply(skew, skew)));
-    for (std::size_t row = 0U; row < 3U; ++row) {
-      for (std::size_t column = 0U; column < 3U; ++column) {
-        result[row][column].directional_first =
-            stable[row][column].directional_first;
-        result[row][column].directional_second =
-            stable[row][column].directional_second;
-        result[row][column].mixed_first = stable[row][column].mixed_first;
-        result[row][column].mixed_second = stable[row][column].mixed_second;
-      }
-    }
+    result = stable;
  }
  return result;
 }
@@ -845,7 +883,7 @@ JetMatrix3 ExpSkew(const JetVector3& rotation) {
 struct JetTrial {
  std::array<JetVector3, kNodeCount> displacement;
  std::array<JetVector3, kNodeCount> rotation;
-  std::array<JetVector3, kNodeCount> director;
+  std::array<JetVector3, kNodeCount> director_increment;
 };

 /// @brief Lifts one validated nonlinear trial into independent Jet variables.
@@ -871,8 +909,8 @@ JetTrial MakeJetTrial(const NonlinearMitc4TrialDofs& trial,
    }
    const JetVector3 initial_director{frame[0U][2U], frame[1U][2U],
                                      frame[2U][2U]};
-    result.director[node] =
-        Multiply(ExpSkew(result.rotation[node]), initial_director);
+    result.director_increment[node] =
+        Multiply(ExpSkew(result.rotation[node], true), initial_director);
  }
  return result;
 }
@@ -882,13 +920,13 @@ std::array<JetVector3, 3> CovariantIncrement(
    const std::array<Vector3, kNodeCount>& initial_directors,
    const JetTrial& trial, const double thickness, const double xi,
    const double eta, const double zeta) {
+  static_cast<void>(initial_directors);
  const auto shape = Mitc4Shell::ShapeFunctions(xi, eta);
  std::array<JetVector3, 3> increment{};
  for (std::size_t node = 0U; node < kNodeCount; ++node) {
    const auto displacement =
        Subtract(trial.displacement[node], trial.displacement[0U]);
-    const auto director =
-        Subtract(trial.director[node], ConstantVector(initial_directors[node]));
+    const auto director = trial.director_increment[node];
    const auto layer =
        Add(displacement, Scale(0.5 * thickness * zeta, director));
    increment[0U] =
@@ -1043,8 +1081,8 @@ JetMatrix3 PolarRotation(const JetMatrix3& deformation) {
  return result;
 }

-/// @brief Returns the speed of physical spin for a unit chart null vector.
-Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) {
+/// @brief Returns the unnormalized physical drill-null chart vector.
+JetVector3 DrillNullPreimage(const JetVector3& phi, const Vector3& director) {
  const Jet square = Dot(phi, phi);
  Jet coefficient{};
  if (square.value < 0.25) {
@@ -1069,7 +1107,19 @@ Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) {
  const auto cross = Cross(phi, d0);
  const auto inverse_action =
      Add(Add(d0, Scale(0.5, cross)), Scale(coefficient, Cross(phi, cross)));
-  return Reciprocal(Sqrt(Dot(inverse_action, inverse_action)));
+  return inverse_action;
+}
+
+/// @brief Returns the normalized physical drill-null chart vector.
+JetVector3 DrillNull(const JetVector3& phi, const Vector3& director) {
+  const JetVector3 preimage = DrillNullPreimage(phi, director);
+  return Scale(Reciprocal(Sqrt(Dot(preimage, preimage))), preimage);
+}
+
+/// @brief Returns the speed of physical spin for a unit chart null vector.
+Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) {
+  const JetVector3 preimage = DrillNullPreimage(phi, director);
+  return Reciprocal(Sqrt(Dot(preimage, preimage)));
 }

 /// @brief Evaluates the objective numerical drilling potential and variations.
@@ -1332,14 +1382,17 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(
      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])) {
+    const auto point_status = EvaluateTyingSample(
+        coordinates_, initial_directors_, directors, trial, thickness_,
+        tying_locations[point][0U], tying_locations[point][1U],
+        tying_samples[point]);
+    if (point_status != PointStatus::kValid) {
      return Failure<NonlinearMitc4Recovery>(
-          "invalid-nonlinear-shell-geometry", source_location_, identity_,
-          "Nonlinear MITC4 tying-point reference or trial geometry is "
-          "invalid.");
+          point_status == PointStatus::kInvalidGeometry
+              ? "invalid-nonlinear-shell-geometry"
+              : "invalid-nonlinear-shell-recovery",
+          source_location_, identity_,
+          "Nonlinear MITC4 tying-point geometry or strain is invalid.");
    }
  }

@@ -1374,14 +1427,18 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(

    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)) {
+      const auto point_status = EvaluatePoint(
+          coordinates_, initial_directors_, directors, trial, tying_samples,
+          thickness_, youngs_modulus_, poisson_ratio_,
+          row.natural_coordinates[0U], row.natural_coordinates[1U],
+          thickness_sign * gauss, point);
+      if (point_status != PointStatus::kValid) {
        return Failure<NonlinearMitc4Recovery>(
-            "invalid-nonlinear-shell-geometry", source_location_, identity_,
-            "Nonlinear MITC4 volume recovery geometry is invalid.");
+            point_status == PointStatus::kInvalidGeometry
+                ? "invalid-nonlinear-shell-geometry"
+                : "invalid-nonlinear-shell-recovery",
+            source_location_, identity_,
+            "Nonlinear MITC4 volume recovery geometry or values are invalid.");
      }
      for (std::size_t component = 0U; component < 3U; ++component) {
        row.generalized_strain[component] += 0.5 * point.strain[component];
@@ -1410,14 +1467,18 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::EvaluateKinematics(

    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)) {
+      const auto point_status = EvaluatePoint(
+          coordinates_, initial_directors_, directors, trial, tying_samples,
+          thickness_, youngs_modulus_, poisson_ratio_,
+          row.natural_coordinates[0U], row.natural_coordinates[1U],
+          stress_zeta[position], point);
+      if (point_status != PointStatus::kValid) {
        return Failure<NonlinearMitc4Recovery>(
-            "invalid-nonlinear-shell-geometry", source_location_, identity_,
-            "Nonlinear MITC4 section-position recovery geometry is invalid.");
+            point_status == PointStatus::kInvalidGeometry
+                ? "invalid-nonlinear-shell-geometry"
+                : "invalid-nonlinear-shell-recovery",
+            source_location_, identity_,
+            "Nonlinear MITC4 section recovery geometry or values are invalid.");
      }
      row.stress[position] = {
          stress_position[position],
@@ -1439,6 +1500,97 @@ Result<NonlinearMitc4Recovery> NonlinearMitc4Shell::RecoverState(
  return EvaluateKinematics(trial);
 }

+Result<NonlinearMitc4ExternalDrillingConstraint>
+NonlinearMitc4Shell::EvaluateExternalDrillingConstraint(
+    const std::size_t node, const Vector3& phi, const Vector3& phi_dot,
+    const Vector3& force, const Vector3& force_rate,
+    const Vector3& force_acceleration) const {
+  if (node >= kNodeCount || !phi_dot.IsFinite() || !force.IsFinite() ||
+      !force_rate.IsFinite() || !force_acceleration.IsFinite()) {
+    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
+        "invalid-nonlinear-shell-external-constraint", source_location_,
+        identity_, "External drilling constraint inputs must be finite and "
+                   "identify one source-local node.");
+  }
+  const auto rotation = MakeFiniteRotationState(phi, initial_frames_[node]);
+  if (!rotation.HasValue()) {
+    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
+        "invalid-nonlinear-shell-rotation", source_location_, identity_,
+        "External drilling constraint rotations must remain in the principal "
+        "chart.");
+  }
+
+  JetVector3 jet_phi{};
+  for (std::size_t component = 0U; component < 3U; ++component) {
+    jet_phi[component] = Variable(phi[component], component,
+                                  phi_dot[component]);
+  }
+  const JetVector3 null_direction = DrillNull(jet_phi, initial_directors_[node]);
+  const Jet scalar = -1.0 * Dot(null_direction, ConstantVector(force));
+  const Jet scalar_time =
+      -1.0 * Dot(null_direction, ConstantVector(force_rate));
+  const Jet scalar_acceleration =
+      -1.0 * Dot(null_direction, ConstantVector(force_acceleration));
+  if (!IsFinite(scalar) || !IsFinite(scalar_time) ||
+      !IsFinite(scalar_acceleration)) {
+    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
+        "invalid-nonlinear-shell-external-constraint", source_location_,
+        identity_, "External drilling constraint derivatives must remain "
+                   "finite.");
+  }
+
+  NonlinearMitc4ExternalDrillingConstraint result{};
+  result.value = scalar.value;
+  result.gradient = Vector{3U};
+  result.hessian = Matrix{3U, 3U};
+  result.gradient_directional_first = Vector{3U};
+  result.gradient_directional_second = Vector{3U};
+  result.directional_first =
+      scalar.directional_first + scalar_time.value;
+  result.directional_second = scalar.directional_second +
+                              2.0 * scalar_time.directional_first +
+                              scalar_acceleration.value;
+  for (std::size_t row = 0U; row < 3U; ++row) {
+    result.gradient[row] = scalar.first[row];
+    result.gradient_directional_first[row] =
+        scalar.mixed_first[row] + scalar_time.first[row];
+    result.gradient_directional_second[row] =
+        scalar.mixed_second[row] +
+        2.0 * scalar_time.mixed_first[row] + scalar_acceleration.first[row];
+    for (std::size_t column = 0U; column < 3U; ++column) {
+      result.hessian(row, column) =
+          scalar.second[row * kGlobalDofCount + column];
+    }
+  }
+  if (!std::isfinite(result.value) || !std::isfinite(result.directional_first) ||
+      !std::isfinite(result.directional_second)) {
+    return Failure<NonlinearMitc4ExternalDrillingConstraint>(
+        "invalid-nonlinear-shell-external-constraint", source_location_,
+        identity_, "External drilling constraint derivatives must remain "
+                   "finite.");
+  }
+  for (std::size_t row = 0U; row < 3U; ++row) {
+    if (!std::isfinite(result.gradient[row]) ||
+        !std::isfinite(result.gradient_directional_first[row]) ||
+        !std::isfinite(result.gradient_directional_second[row])) {
+      return Failure<NonlinearMitc4ExternalDrillingConstraint>(
+          "invalid-nonlinear-shell-external-constraint", source_location_,
+          identity_, "External drilling constraint derivatives must remain "
+                     "finite.");
+    }
+    for (std::size_t column = 0U; column < 3U; ++column) {
+      if (!std::isfinite(result.hessian(row, column))) {
+        return Failure<NonlinearMitc4ExternalDrillingConstraint>(
+            "invalid-nonlinear-shell-external-constraint", source_location_,
+            identity_, "External drilling constraint derivatives must remain "
+                       "finite.");
+      }
+    }
+  }
+  return Result<NonlinearMitc4ExternalDrillingConstraint>::Success(
+      std::move(result));
+}
+
 Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
    const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial,
    const Vector* coordinate_direction) {
@@ -1526,9 +1678,13 @@ Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
      equilibrium.drilling_constraint[node] = constraints[node].value;
      equilibrium.drilling_constraint_hessian[node] =
          Matrix{kGlobalDofCount, kGlobalDofCount};
+      equilibrium.drilling_constraint_third_directional_gradient[node] =
+          Vector{kGlobalDofCount};
      for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
        equilibrium.drilling_constraint_gradient(node, i) =
            constraints[node].first[i];
+        equilibrium.drilling_constraint_third_directional_gradient[node][i] =
+            constraints[node].mixed_second[i];
        for (std::size_t j = 0U; j < kGlobalDofCount; ++j)
          equilibrium.drilling_constraint_hessian[node](i, j) =
              constraints[node].second[i * kGlobalDofCount + j];