Add immutable MITC4 drilling selection and direct director springs

src/fesa/elements/mitc4_drilling_internal.h

@@ -0,0 +1,28 @@
+#ifndef FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_
+#define FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_
+
+#include "fesa/elements/nonlinear_mitc4_shell.h"
+
+namespace fesa::mitc4_drilling_internal {
+
+// The sole production selection. Change this constant and rebuild all targets.
+inline constexpr Mitc4DrillingMode kProductionDrillingMode =
+    Mitc4DrillingMode::kDirectSpring;
+
+/// @brief Provides explicit alternate formulations only for internal tests.
+/// @note Production factories use the single compile-time selection above.
+struct TestAccess {
+  static Mitc4Shell WithMode(Mitc4Shell shell, Mitc4DrillingMode mode) {
+    shell.drilling_mode_ = mode;
+    return shell;
+  }
+  static NonlinearMitc4Shell WithMode(NonlinearMitc4Shell shell,
+                                      Mitc4DrillingMode mode) {
+    shell.drilling_mode_ = mode;
+    return shell;
+  }
+};
+
+}  // namespace fesa::mitc4_drilling_internal
+
+#endif  // FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_

src/fesa/elements/mitc4_shell.cpp

@@ -12,6 +12,7 @@
 #include "math/finite_values.h"
 #include "math/vector3_operations.h"
 #include "elements/shell_kernel_values.h"
+#include "elements/mitc4_drilling_internal.h"
 
 namespace fesa {
 
@@ -1002,6 +1003,7 @@ Mitc4Shell::Mitc4Shell(std::array<Vector3, 4> coordinates,
       thickness_{thickness},
       youngs_modulus_{youngs_modulus},
       poisson_ratio_{poisson_ratio},
+      drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode},
       source_location_{std::move(source_location)},
       identity_{std::move(identity)} {}
 

src/fesa/elements/nonlinear_mitc4_configuration_internal.h

@@ -12,6 +12,7 @@ namespace fesa::nonlinear_mitc4_internal {
 /// key.
 struct ConfigurationSnapshot {
   std::array<double, 64> values{};
+  Mitc4DrillingMode drilling_mode = Mitc4DrillingMode::kDirectSpring;
   SourceLocation source{};
   std::string identity;
 

src/fesa/elements/nonlinear_mitc4_shell.cpp

@@ -12,6 +12,7 @@
 
 #include "fesa/elements/finite_rotation.h"
 #include "elements/nonlinear_mitc4_configuration_internal.h"
+#include "elements/mitc4_drilling_internal.h"
 #include "math/dense_blas_internal.h"
 #include "math/jet_arithmetic.h"
 #include "math/jet_linear_solve.h"
@@ -39,6 +40,7 @@ nonlinear_mitc4_internal::ConfigurationSnapshot::Read(
   result.values[index++] = element.youngs_modulus_;
   result.values[index++] = element.poisson_ratio_;
   result.values[index] = element.drilling_stiffness_;
+  result.drilling_mode = element.drilling_mode_;
   result.source = element.source_location_;
   result.identity = element.identity_;
   return result;
@@ -1227,8 +1229,91 @@ Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) {
   return Reciprocal(Sqrt(Dot(preimage, preimage)));
 }
 
+/// @brief Uses precisely the current inertia kernel's cubic rotation maps.
+/// @note This is deliberately distinct from the legacy analytic null. Keep the
+/// cubic Exp/Jacobian order synchronized with EvaluateDynamicRotation; it is
+/// not an approximation to be silently replaced by the exact exponential
+/// inverse.
+JetVector3 InertiaDrillNull(const JetVector3& phi, const Vector3& director) {
+  const JetMatrix3 skew = Skew(phi);
+  JetMatrix3 rotation = Identity();
+  JetMatrix3 jacobian = Identity();
+  JetMatrix3 term = Identity();
+  double rotation_factorial = 1.0;
+  double jacobian_factorial = 1.0;
+  for (std::size_t order = 1U; order <= 3U; ++order) {
+    term = Multiply(term, skew);
+    rotation_factorial *= static_cast<double>(order);
+    jacobian_factorial *= static_cast<double>(order + 1U);
+    rotation = Add(rotation, Scale(1.0 / rotation_factorial, term));
+    jacobian = Add(jacobian, Scale(1.0 / jacobian_factorial, term));
+  }
+  const JetVector3 rhs = Multiply(rotation, ConstantVector(director));
+  Matrix primal{3U, 3U};
+  for (std::size_t row = 0U; row < 3U; ++row)
+    for (std::size_t column = 0U; column < 3U; ++column)
+      primal(row, column) = jacobian[row][column].value;
+  const auto factor = dense_blas_internal::DenseFactor::Lu(primal);
+  const std::array<std::size_t, 4> channels{1U, kGlobalDofCount + 1U,
+                                            kHessianSize + 1U + kGlobalDofCount,
+                                            kGlobalDofCount};
+  JetVector3 preimage{};
+  if (!factor.HasValue() ||
+      !math_internal::SolveJetAction(jacobian, rhs, factor.Value(), channels,
+                                     Constant, ReadChannel, WriteChannel,
+                                     preimage)) {
+    preimage.fill(Constant(std::numeric_limits<double>::quiet_NaN()));
+    return preimage;
+  }
+  return Scale(Reciprocal(Sqrt(Dot(preimage, preimage))), preimage);
+}
+
+/// @brief Evaluates the fixed initial-director spring without corotation.
+Jet DirectDrillingEnergy(
+    const std::array<Vector3, kNodeCount>& initial_directors,
+    const JetTrial& trial, double drilling_stiffness, Matrix* metric,
+    std::array<Jet, 4>* constraints) {
+  Jet energy{};
+  double* tangent = energy.second.Data();
+  for (std::size_t node = 0U; node < kNodeCount; ++node) {
+    const auto d0 = ConstantVector(initial_directors[node]);
+    const Jet gamma = Dot(d0, trial.rotation[node]);
+    // Retain the validated spring scale through subnormal intermediates.
+    // Fill its linear derivatives explicitly, without differentiating a
+    // floating half-energy that can round kd before doubling it again.
+    energy.value += 0.5 * gamma.value * (drilling_stiffness * gamma.value);
+    energy.directional_first +=
+        drilling_stiffness * gamma.value * gamma.directional_first;
+    energy.directional_second +=
+        drilling_stiffness * gamma.directional_first * gamma.directional_first;
+    for (std::size_t i = 0U; i < 3U; ++i) {
+      const auto row = 6U * node + 3U + i;
+      energy.first[row] =
+          drilling_stiffness * gamma.value * initial_directors[node][i];
+      energy.mixed_first[row] = drilling_stiffness * gamma.directional_first *
+                                initial_directors[node][i];
+    }
+    // D2 f_drill[v,v] is identically zero; mixed_second retains exact zeros.
+    if (constraints != nullptr) {
+      const auto null =
+          InertiaDrillNull(trial.rotation[node], initial_directors[node]);
+      (*constraints)[node] = drilling_stiffness * gamma * Dot(null, d0);
+    }
+    for (std::size_t i = 0U; i < 3U; ++i)
+      for (std::size_t j = 0U; j < 3U; ++j) {
+        const auto row = 6U * node + 3U + i;
+        const auto column = 6U * node + 3U + j;
+        const double value = drilling_stiffness * initial_directors[node][i] *
+                             initial_directors[node][j];
+        tangent[row * kGlobalDofCount + column] = value;
+        (*metric)(row, column) = value;
+      }
+  }
+  return energy;
+}
+
 /// @brief Evaluates the objective numerical drilling potential and variations.
-Jet DrillingEnergy(const std::array<Vector3, kNodeCount>& coordinates,
+Jet LegacyDrillingEnergy(const std::array<Vector3, kNodeCount>& coordinates,
                    const std::array<Vector3, kNodeCount>& initial_directors,
                    const JetTrial& trial, const double thickness,
                    const double drilling_stiffness, bool& valid, Matrix* metric,
@@ -1323,6 +1408,7 @@ NonlinearMitc4Shell::NonlinearMitc4Shell(
       youngs_modulus_{youngs_modulus},
       poisson_ratio_{poisson_ratio},
       drilling_stiffness_{drilling_stiffness},
+      drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode},
       source_location_{std::move(source_location)},
       identity_{std::move(identity)} {}
 
@@ -1640,7 +1726,10 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint(
     jet_phi[component] = Variable(phi[component], component,
                                   phi_dot[component]);
   }
-  const JetVector3 null_direction = DrillNull(jet_phi, initial_directors_[node]);
+  const JetVector3 null_direction =
+      drilling_mode_ == Mitc4DrillingMode::kDirectSpring
+          ? InertiaDrillNull(jet_phi, initial_directors_[node])
+          : 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));
@@ -1741,10 +1830,17 @@ Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
                      element.poisson_ratio_, valid);
   Matrix drilling_metric{kGlobalDofCount, kGlobalDofCount};
   std::array<Jet, 4> constraints{};
-  const Jet drilling = DrillingEnergy(
-      element.coordinates_, element.initial_directors_, jet_trial,
-      element.thickness_, element.drilling_stiffness_, valid, &drilling_metric,
-      coordinate_direction != nullptr ? &constraints : nullptr);
+  const Jet drilling =
+      element.drilling_mode_ == Mitc4DrillingMode::kDirectSpring
+          ? DirectDrillingEnergy(
+                element.initial_directors_, jet_trial,
+                element.drilling_stiffness_, &drilling_metric,
+                coordinate_direction != nullptr ? &constraints : nullptr)
+          : LegacyDrillingEnergy(
+                element.coordinates_, element.initial_directors_, jet_trial,
+                element.thickness_, element.drilling_stiffness_, valid,
+                &drilling_metric,
+                coordinate_direction != nullptr ? &constraints : nullptr);
   const Jet internal_energy = physical + drilling;
   if (!valid || !IsFinite(internal_energy)) {
     return Failure<NonlinearMitc4InternalEquilibrium>(
@@ -1815,10 +1911,13 @@ Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
   equilibrium.drilling_metric = std::move(drilling_metric);
   equilibrium.drilling_energy = drilling.value;
   equilibrium.force_directional_second = Vector{kGlobalDofCount};
+  equilibrium.drilling_force_directional_second = Vector{kGlobalDofCount};
   for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
     equilibrium.internal_force[index] = equilibrium.internal_residual[index];
     equilibrium.physical_force[index] = physical.first[index];
     equilibrium.drilling_force[index] = drilling.first[index];
+    equilibrium.drilling_force_directional_second[index] =
+        drilling.mixed_second[index];
     equilibrium.force_directional_second[index] =
         internal_energy.mixed_second[index];
     for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {