Merge reviewed MKL shell products and complete-trial reuse

src/fesa/elements/mitc4_shell.cpp

@@ -8,6 +8,7 @@
 #include <utility>
 
 #include "fesa/model/model_types.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 
@@ -223,28 +224,18 @@ Matrix ScaledMatrix(const Matrix& source, double factor) {
   return result;
 }
 
-/// @brief Lifts a local matrix by a deterministic serial congruence product.
-/// @note The loop and accumulation order are part of reproducible stiffness.
+/// @brief Lifts a local matrix through MKL and mirrors the upper triangle.
 Matrix Congruence(const Matrix& local, const Matrix& transformation) {
   if (local.Rows() != local.Columns() ||
       local.Rows() != transformation.Rows()) {
     throw std::invalid_argument{
         "MITC4 congruence dimensions are incompatible."};
   }
-  Matrix result{transformation.Columns(), transformation.Columns()};
+  Matrix result = dense_blas_internal::MultiplyTransposeLeft(
+      transformation, local.Multiply(transformation));
   for (std::size_t row = 0U; row < result.Rows(); ++row) {
     for (std::size_t column = row; column < result.Columns(); ++column) {
-      double value = 0.0;
-      for (std::size_t local_row = 0U; local_row < local.Rows(); ++local_row) {
-        for (std::size_t local_column = 0U; local_column < local.Columns();
-             ++local_column) {
-          value += transformation(local_row, row) *
-                   local(local_row, local_column) *
-                   transformation(local_column, column);
-        }
-      }
-      result(row, column) = value;
-      result(column, row) = value;
+      result(column, row) = result(row, column);
     }
   }
   return result;
@@ -611,7 +602,8 @@ Matrix Mitc4Shell::TransverseShearSectionMatrix() const {
   return result;
 }
 
-Result<Mitc4Stiffness> Mitc4Shell::Stiffness() const {
+Result<Mitc4Stiffness> Mitc4Shell::Stiffness() const try {
+  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
   const Matrix constitutive = MaterialConstitutive5();
   const Matrix tying_samples = CovariantTyingShearSamples20();
   Matrix physical_local{kPhysicalDofCount, kPhysicalDofCount};
@@ -627,17 +619,12 @@ Result<Mitc4Stiffness> Mitc4Shell::Stiffness() const {
     const Matrix strain = StrainDisplacement(
         point.natural_coordinates[0], point.natural_coordinates[1],
         point.natural_coordinates[2], &tying_samples);
+    const Matrix integrand = dense_blas_internal::MultiplyTransposeLeft(
+        strain, constitutive.Multiply(strain));
     for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) {
       for (std::size_t column = row; column < kPhysicalDofCount; ++column) {
-        double integrand = 0.0;
-        for (std::size_t first = 0U; first < 5U; ++first) {
-          for (std::size_t second = 0U; second < 5U; ++second) {
-            integrand += strain(first, row) * constitutive(first, second) *
-                         strain(second, column);
-          }
-        }
         const double contribution =
-            integrand * geometry.jacobian * point.weight;
+            integrand(row, column) * geometry.jacobian * point.weight;
         physical_local(row, column) += contribution;
         if (row != column) {
           physical_local(column, row) += contribution;
@@ -714,6 +701,8 @@ Result<Mitc4Stiffness> Mitc4Shell::Stiffness() const {
       Mitc4Stiffness{std::move(physical_local), std::move(physical_global),
                      std::move(drilling_global), std::move(stabilized_global),
                      drilling_stiffness});
+} catch (const std::exception& error) {
+  return StiffnessFailure(source_location_, identity_, error.what());
 }
 
 Result<Mitc4Mass> Mitc4Shell::Mass(const double density) const {

src/fesa/elements/nonlinear_mitc4_configuration_internal.h

@@ -0,0 +1,24 @@
+#ifndef FESA_SRC_FESA_ELEMENTS_NONLINEAR_MITC4_CONFIGURATION_INTERNAL_H_
+#define FESA_SRC_FESA_ELEMENTS_NONLINEAR_MITC4_CONFIGURATION_INTERNAL_H_
+
+#include <array>
+#include <string>
+
+#include "fesa/elements/nonlinear_mitc4_shell.h"
+
+namespace fesa::nonlinear_mitc4_internal {
+
+/// @brief Owns reference values without depending on an analysis or identity
+/// key.
+struct ConfigurationSnapshot {
+  std::array<double, 64> values{};
+  SourceLocation source{};
+  std::string identity;
+
+  /// @brief Copies all immutable inputs of the internal equilibrium evaluator.
+  static ConfigurationSnapshot Read(const NonlinearMitc4Shell& element);
+};
+
+}  // namespace fesa::nonlinear_mitc4_internal
+
+#endif  // FESA_SRC_FESA_ELEMENTS_NONLINEAR_MITC4_CONFIGURATION_INTERNAL_H_

src/fesa/elements/nonlinear_mitc4_inertia.cpp

@@ -4,11 +4,13 @@
 #include <array>
 #include <cmath>
 #include <cstddef>
+#include <exception>
 #include <limits>
 #include <string>
 #include <utility>
 
 #include "fesa/elements/finite_rotation.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 namespace {
@@ -35,6 +37,16 @@ class InertiaSum {
   double correction_ = 0.0;
 };
 
+/// @brief Contracts two nodal 6-by-3 coefficient blocks through dense MKL.
+/// @note Nonfinite coefficients remain failures of the final inertia gate.
+void DotBlocks(const std::array<double, 18U>& left,
+               const std::array<double, 18U>& right,
+               std::array<double, 36U>& result) {
+  if (!dense_blas_internal::MultiplyTransposeRight(
+          left.data(), 6U, right.data(), 6U, 3U, result.data()))
+    result.fill(std::numeric_limits<double>::quiet_NaN());
+}
+
 Vector3 WeightedSum(const std::array<double, kNodes>& weights,
                     const std::array<Vector3, kNodes>& values) {
   Vector3 result{};
@@ -442,12 +454,11 @@ Status CheckPhysicalMass(const Matrix& mass,
   }
   double gram_error = 0.0;
   double orthogonality_error = 0.0;
+  const Matrix gram =
+      dense_blas_internal::MultiplyTransposeLeft(complement, complement);
   for (std::size_t row = 0U; row < kPhysical; ++row) {
     for (std::size_t column = 0U; column < kPhysical; ++column) {
-      double entry = row == column ? -1.0 : 0.0;
-      for (std::size_t inner = 0U; inner < kDofs; ++inner) {
-        entry += complement(inner, row) * complement(inner, column);
-      }
+      const double entry = gram(row, column) - (row == column ? 1.0 : 0.0);
       gram_error = std::hypot(gram_error, entry);
     }
     for (std::size_t node = 0U; node < kNodes; ++node) {
@@ -472,55 +483,42 @@ Status CheckPhysicalMass(const Matrix& mass,
     return fail("invalid-nonlinear-inertia-null",
                 "Scaled mass symmetry and analytic null identities failed.");
   }
-  Matrix block{kPhysical, kPhysical};
-  Matrix mass_complement{kDofs, kPhysical};
-  for (std::size_t row = 0U; row < kDofs; ++row) {
-    for (std::size_t column = 0U; column < kPhysical; ++column) {
-      for (std::size_t inner = 0U; inner < kDofs; ++inner) {
-        mass_complement(row, column) +=
-            scaled(row, inner) * complement(inner, column);
-      }
-    }
-  }
+  const Matrix mass_complement = scaled.Multiply(complement);
+  const Matrix block =
+      dense_blas_internal::MultiplyTransposeLeft(complement, mass_complement);
   double block_norm = 0.0;
   for (std::size_t row = 0U; row < kPhysical; ++row) {
     double row_sum = 0.0;
     for (std::size_t column = 0U; column < kPhysical; ++column) {
-      for (std::size_t inner = 0U; inner < kDofs; ++inner) {
-        block(row, column) +=
-            complement(inner, row) * mass_complement(inner, column);
-      }
       row_sum += std::abs(block(row, column));
     }
     block_norm = std::max(block_norm, row_sum);
   }
   Matrix lower{kPhysical, kPhysical};
-  for (std::size_t row = 0U; row < kPhysical; ++row) {
-    for (std::size_t column = 0U; column <= row; ++column) {
-      double entry = block(row, column);
-      for (std::size_t inner = 0U; inner < column; ++inner) {
-        entry -= lower(row, inner) * lower(column, inner);
-      }
-      if (!std::isfinite(entry) || (row == column && !(entry > 0.0))) {
-        return fail(
-            "invalid-nonlinear-inertia-positive-mass",
-            "Physical mass requires strictly positive Cholesky pivots.");
-      }
-      lower(row, column) =
-          row == column ? std::sqrt(entry) : entry / lower(column, column);
-    }
+  Matrix inverse_lower{kPhysical, kPhysical};
+  auto factor_failure = dense_blas_internal::CholeskyFailure::kNone;
+  if (!dense_blas_internal::CholeskyLowerAndInverse(block, lower, inverse_lower,
+                                                    &factor_failure)) {
+    if (factor_failure == dense_blas_internal::CholeskyFailure::kInverse)
+      return fail(
+          "ill-conditioned-nonlinear-inertia-positive-mass",
+          "Physical mass reciprocal infinity condition is below 2^-32.");
+    return fail("invalid-nonlinear-inertia-positive-mass",
+                "Physical mass requires strictly positive Cholesky pivots.");
   }
   std::array<double, kPhysical> inverse_row_sums{};
   double reconstruction_error = 0.0;
   double block_frobenius = 0.0;
+  Matrix transpose{kPhysical, kPhysical};
+  for (std::size_t row = 0U; row < kPhysical; ++row)
+    for (std::size_t column = 0U; column < kPhysical; ++column)
+      transpose(row, column) = lower(column, row);
+  const Matrix reconstructed = lower.Multiply(transpose);
   for (std::size_t row = 0U; row < kPhysical; ++row) {
     for (std::size_t column = 0U; column < kPhysical; ++column) {
-      double reconstructed = 0.0;
-      for (std::size_t inner = 0U; inner <= std::min(row, column); ++inner) {
-        reconstructed += lower(row, inner) * lower(column, inner);
-      }
       reconstruction_error =
-          std::hypot(reconstruction_error, reconstructed - block(row, column));
+          std::hypot(reconstruction_error,
+                     reconstructed(row, column) - block(row, column));
       block_frobenius = std::hypot(block_frobenius, block(row, column));
     }
   }
@@ -528,25 +526,11 @@ Status CheckPhysicalMass(const Matrix& mass,
     return fail("invalid-nonlinear-inertia-positive-mass",
                 "Physical mass Cholesky reconstruction failed.");
   }
-  for (std::size_t column = 0U; column < kPhysical; ++column) {
-    std::array<double, kPhysical> solution{};
-    for (std::size_t row = 0U; row < kPhysical; ++row) {
-      double rhs = row == column ? 1.0 : 0.0;
-      for (std::size_t inner = 0U; inner < row; ++inner) {
-        rhs -= lower(row, inner) * solution[inner];
-      }
-      solution[row] = rhs / lower(row, row);
-    }
-    for (std::size_t reverse = kPhysical; reverse > 0U; --reverse) {
-      const std::size_t row = reverse - 1U;
-      double rhs = solution[row];
-      for (std::size_t inner = row + 1U; inner < kPhysical; ++inner) {
-        rhs -= lower(inner, row) * solution[inner];
-      }
-      solution[row] = rhs / lower(row, row);
-      inverse_row_sums[row] += std::abs(solution[row]);
-    }
-  }
+  const Matrix inverse =
+      dense_blas_internal::MultiplyTransposeLeft(inverse_lower, inverse_lower);
+  for (std::size_t column = 0U; column < kPhysical; ++column)
+    for (std::size_t row = 0U; row < kPhysical; ++row)
+      inverse_row_sums[row] += std::abs(inverse(row, column));
   const double inverse_norm =
       *std::max_element(inverse_row_sums.begin(), inverse_row_sums.end());
   const double reciprocal_condition = 1.0 / block_norm / inverse_norm;
@@ -712,6 +696,7 @@ Result<NonlinearMitc4Inertia> NonlinearMitc4Inertia::Create(
 Result<NonlinearMitc4InertiaEvaluation> NonlinearMitc4Inertia::Evaluate(
     const NonlinearMitc4TrialDofs& trial, const Vector& velocity,
     const Vector& acceleration) const {
+  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
   if (velocity.Size() != kDofs || acceleration.Size() != kDofs ||
       !IsFinite(velocity) || !IsFinite(acceleration)) {
     return Failure<NonlinearMitc4InertiaEvaluation>(
@@ -899,6 +884,42 @@ Result<NonlinearMitc4InertiaEvaluation> NonlinearMitc4Inertia::Evaluate(
           result.momentum_rate_scale += scale * acceleration_point.Norm();
           result.angular_momentum_rate_scale +=
               scale * position.Cross(acceleration_point).Norm();
+          std::array<double, 18U> h_block;
+          std::array<double, 18U> velocity_block;
+          std::array<double, 18U> position_block;
+          for (std::size_t local = 0U; local < 6U; ++local) {
+            const auto field = 6U * mass_node + local;
+            for (std::size_t axis = 0U; axis < 3U; ++axis) {
+              h_block[3U * local + axis] = h[field][axis];
+              velocity_block[3U * local + axis] =
+                  velocity_derivative[field][axis];
+              position_block[3U * local + axis] =
+                  position_derivative[field][axis];
+            }
+          }
+          std::array<double, 36U> mass_block;
+          std::array<double, 36U> velocity_product;
+          std::array<double, 36U> position_product;
+          DotBlocks(h_block, h_block, mass_block);
+          DotBlocks(h_block, velocity_block, velocity_product);
+          DotBlocks(h_block, position_block, position_product);
+          std::array<std::array<double, 36U>, 3U> mass_variation;
+          for (std::size_t direction = 0U; direction < 3U; ++direction) {
+            std::array<double, 18U> derivative_block;
+            for (std::size_t local = 0U; local < 6U; ++local) {
+              const Vector3 derivative =
+                  local < 3U ? Vector3{}
+                             : z * rotations[mass_node]
+                                       .director_second[local - 3U][direction];
+              for (std::size_t axis = 0U; axis < 3U; ++axis)
+                derivative_block[3U * local + axis] = derivative[axis];
+            }
+            std::array<double, 36U> reverse;
+            DotBlocks(derivative_block, h_block, mass_variation[direction]);
+            DotBlocks(h_block, derivative_block, reverse);
+            for (std::size_t index = 0U; index < 36U; ++index)
+              mass_variation[direction][index] += reverse[index];
+          }
           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.
@@ -906,10 +927,11 @@ Result<NonlinearMitc4InertiaEvaluation> NonlinearMitc4Inertia::Evaluate(
                 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]);
+              const auto local_index = 6U * (row % 6U) + column % 6U;
+              result.mass(row, column) += scale * mass_block[local_index];
               result.velocity_tangent(row, column) +=
-                  scale * h[row].Dot(velocity_derivative[column]);
-              double position_entry = h[row].Dot(position_derivative[column]);
+                  scale * velocity_product[local_index];
+              double position_entry = position_product[local_index];
               if (row / 6U == column / 6U && row % 6U >= 3U &&
                   column % 6U >= 3U) {
                 const std::size_t node = row / 6U;
@@ -922,20 +944,10 @@ Result<NonlinearMitc4InertiaEvaluation> NonlinearMitc4Inertia::Evaluate(
               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));
+                      scale *
+                      mass_variation[direction][6U * (row % 6U) + local];
                 }
               }
             }
@@ -1004,8 +1016,15 @@ Result<NonlinearMitc4InertiaEvaluation> NonlinearMitc4Inertia::Evaluate(
         "invalid-nonlinear-inertia-evaluation", source_location_, identity_,
         "Nonlinear MITC4 inertia evaluation must remain finite.");
   }
-  const Status physical_mass = CheckPhysicalMass(result.mass, nulls, length,
-                                                 source_location_, identity_);
+  Status physical_mass = Status::Ok();
+  try {
+    physical_mass = CheckPhysicalMass(result.mass, nulls, length,
+                                      source_location_, identity_);
+  } catch (const std::exception& error) {
+    return Failure<NonlinearMitc4InertiaEvaluation>(
+        "invalid-nonlinear-inertia-positive-mass", source_location_, identity_,
+        error.what());
+  }
   if (!physical_mass.IsOk()) {
     return Result<NonlinearMitc4InertiaEvaluation>::Failure(physical_mass);
   }

src/fesa/elements/nonlinear_mitc4_shell.cpp

@@ -10,8 +10,34 @@
 #include <utility>
 
 #include "fesa/elements/finite_rotation.h"
+#include "elements/nonlinear_mitc4_configuration_internal.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
+
+nonlinear_mitc4_internal::ConfigurationSnapshot
+nonlinear_mitc4_internal::ConfigurationSnapshot::Read(
+    const NonlinearMitc4Shell& element) {
+  ConfigurationSnapshot result{};
+  std::size_t index = 0U;
+  for (const auto& coordinates : element.coordinates_)
+    for (std::size_t c = 0U; c < 3U; ++c)
+      result.values[index++] = coordinates[c];
+  for (const auto& frame : element.initial_frames_)
+    for (std::size_t i = 0U; i < 3U; ++i)
+      for (std::size_t j = 0U; j < 3U; ++j)
+        result.values[index++] = frame(i, j);
+  for (const auto& director : element.initial_directors_)
+    for (std::size_t c = 0U; c < 3U; ++c) result.values[index++] = director[c];
+  result.values[index++] = element.thickness_;
+  result.values[index++] = element.youngs_modulus_;
+  result.values[index++] = element.poisson_ratio_;
+  result.values[index] = element.drilling_stiffness_;
+  result.source = element.source_location_;
+  result.identity = element.identity_;
+  return result;
+}
+
 namespace {
 
 constexpr std::size_t kNodeCount = 4U;
@@ -394,22 +420,66 @@ 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>>()} {}
+  Hessian() = default;
 
-  double& operator[](const std::size_t index) { return (*values_)[index]; }
+  double& operator[](const std::size_t index) { return Data()[index]; }
 
   const double& operator[](const std::size_t index) const {
-    return (*values_)[index];
+    return Storage()[index];
   }
 
-  auto begin() const { return values_->begin(); }
+  auto begin() const { return Storage().begin(); }
 
-  auto end() const { return values_->end(); }
+  auto end() const { return Storage().end(); }
+
+  double* Data() {
+    if (!values_) {
+      values_ = std::make_shared<Values>();
+      values_->entries.fill(0.0);
+    }
+    return values_->entries.data();
+  }
+
+  const double* ReadData() const { return Storage().data(); }
+
+  /// @brief Identifies only the immutable zero Hessian, independently of seeds.
+  bool IsImplicitZero() const { return !values_; }
+
+  /// @brief Gets fresh output storage without redundant zero initialization.
+  /// @note Only fresh operator results use this pointer; all 576 entries must
+  /// be written before any read or copy of the result.
+  double* OverwriteData() {
+    if (!values_) values_ = std::make_shared<Values>();
+    return values_->entries.data();
+  }
 
  private:
-  std::shared_ptr<std::array<double, kHessianSize>> values_;
+  struct Values {
+    // User-provided construction leaves entries uninitialized for a complete
+    // output overwrite while make_shared retains the single allocation.
+    Values() noexcept {}
+    std::array<double, kHessianSize> entries;
+  };
+
+  const std::array<double, kHessianSize>& Storage() const {
+    // Constants borrow only this immutable zero. Mutable coefficients retain
+    // the existing shared ownership; no trial-dependent cache is introduced.
+    static const std::array<double, kHessianSize> kZero{};
+    return values_ ? values_->entries : kZero;
+  }
+
+  std::shared_ptr<Values> values_;
 };
 
+/// @brief Forms an ordinary derivative outer product without a Matrix copy.
+/// @note Rejected nonfinite operands propagate to the existing Jet finite gate.
+void Outer(const std::array<double, kGlobalDofCount>& left,
+           const std::array<double, kGlobalDofCount>& right, double* result) {
+  if (!dense_blas_internal::OuterProduct(left.data(), left.size(), right.data(),
+                                         right.size(), result))
+    std::fill_n(result, kHessianSize, std::numeric_limits<double>::quiet_NaN());
+}
+
 /// @brief Carries exact first and second trial-coordinate variations.
 struct Jet {
   double value{};
@@ -419,6 +489,9 @@ struct Jet {
   double directional_second{};
   std::array<double, kGlobalDofCount> mixed_first{};
   std::array<double, kGlobalDofCount> mixed_second{};
+  // Provenance for numerically zero derivatives; a null Hessian alone is not
+  // constant provenance because Variable retains a nonzero gradient seed.
+  bool constant_derivatives{false};
 };
 
 using JetVector3 = std::array<Jet, 3>;
@@ -429,6 +502,7 @@ using JetTensor = std::array<std::array<Jet, 3>, 3>;
 Jet Constant(const double value) {
   Jet result{};
   result.value = value;
+  result.constant_derivatives = true;
   return result;
 }
 
@@ -436,6 +510,7 @@ Jet Constant(const double value) {
 Jet Variable(const double value, const std::size_t index,
              const double direction = 0.0) {
   Jet result = Constant(value);
+  result.constant_derivatives = false;
   result.first[index] = 1.0;
   result.directional_first = direction;
   return result;
@@ -455,8 +530,16 @@ Jet operator+(const Jet& left, const Jet& right) {
     result.mixed_second[index] =
         left.mixed_second[index] + right.mixed_second[index];
   }
+  result.constant_derivatives = left.constant_derivatives &&
+                                right.constant_derivatives &&
+                                std::isfinite(result.value);
+  if (left.second.IsImplicitZero() && right.second.IsImplicitZero())
+    return result;
+  const double* left_second = left.second.ReadData();
+  const double* right_second = right.second.ReadData();
+  double* result_second = result.second.OverwriteData();
   for (std::size_t index = 0U; index < kHessianSize; ++index) {
-    result.second[index] = left.second[index] + right.second[index];
+    result_second[index] = left_second[index] + right_second[index];
   }
   return result;
 }
@@ -475,8 +558,16 @@ Jet operator-(const Jet& left, const Jet& right) {
     result.mixed_second[index] =
         left.mixed_second[index] - right.mixed_second[index];
   }
+  result.constant_derivatives = left.constant_derivatives &&
+                                right.constant_derivatives &&
+                                std::isfinite(result.value);
+  if (left.second.IsImplicitZero() && right.second.IsImplicitZero())
+    return result;
+  const double* left_second = left.second.ReadData();
+  const double* right_second = right.second.ReadData();
+  double* result_second = result.second.OverwriteData();
   for (std::size_t index = 0U; index < kHessianSize; ++index) {
-    result.second[index] = left.second[index] - right.second[index];
+    result_second[index] = left_second[index] - right_second[index];
   }
   return result;
 }
@@ -506,13 +597,38 @@ Jet operator*(const Jet& left, const Jet& right) {
         2.0 * left.directional_first * right.mixed_first[index] +
         left.value * right.mixed_second[index];
   }
+  const bool finite_values = std::isfinite(left.value) &&
+                             std::isfinite(right.value) &&
+                             std::isfinite(result.value);
+  const auto finite_gradient = [](const auto& gradient) {
+    return std::all_of(gradient.begin(), gradient.end(),
+                       [](const double value) { return std::isfinite(value); });
+  };
+  const bool zero_outer = finite_values &&
+      ((left.constant_derivatives && finite_gradient(right.first)) ||
+       (right.constant_derivatives && finite_gradient(left.first)));
+  result.constant_derivatives = finite_values && left.constant_derivatives &&
+                                right.constant_derivatives;
+  if (zero_outer && left.second.IsImplicitZero() &&
+      right.second.IsImplicitZero()) return result;
+  const double* left_second = left.second.ReadData();
+  const double* right_second = right.second.ReadData();
+  double* result_second = result.second.OverwriteData();
+  if (zero_outer) {
+    for (std::size_t index = 0U; index < kHessianSize; ++index)
+      result_second[index] = left_second[index] * right.value + 0.0 + 0.0 +
+                             left.value * right_second[index];
+    return result;
+  }
+  Outer(left.first, right.first, result_second);
+  std::array<double, kHessianSize> reverse;
+  Outer(right.first, left.first, reverse.data());
   for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
     for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
       const std::size_t index = row * kGlobalDofCount + column;
-      result.second[index] = left.second[index] * right.value +
-                             left.first[row] * right.first[column] +
-                             left.first[column] * right.first[row] +
-                             left.value * right.second[index];
+      result_second[index] = left_second[index] * right.value +
+                             result_second[index] + reverse[index] +
+                             left.value * right_second[index];
     }
   }
   return result;
@@ -529,8 +645,14 @@ Jet operator*(const double scalar, const Jet& value) {
     result.mixed_first[index] = scalar * value.mixed_first[index];
     result.mixed_second[index] = scalar * value.mixed_second[index];
   }
+  result.constant_derivatives = value.constant_derivatives &&
+                                std::isfinite(scalar) &&
+                                std::isfinite(result.value);
+  if (std::isfinite(scalar) && value.second.IsImplicitZero()) return result;
+  const double* value_second = value.second.ReadData();
+  double* result_second = result.second.OverwriteData();
   for (std::size_t index = 0U; index < kHessianSize; ++index) {
-    result.second[index] = scalar * value.second[index];
+    result_second[index] = scalar * value_second[index];
   }
   return result;
 }
@@ -559,11 +681,20 @@ Jet Compose(const Jet& input, const double value, const double first,
                   input.first[index] * input.directional_second) +
         first * input.mixed_second[index];
   }
+  result.constant_derivatives = input.constant_derivatives &&
+      std::isfinite(input.value) && std::isfinite(value) &&
+      std::isfinite(first) && std::isfinite(second) && std::isfinite(third);
+  if (result.constant_derivatives) return result;
+  std::array<double, kGlobalDofCount> scaled_first{};
+  for (std::size_t row = 0U; row < kGlobalDofCount; ++row)
+    scaled_first[row] = second * input.first[row];
+  const double* input_second = input.second.ReadData();
+  double* result_second = result.second.OverwriteData();
+  Outer(scaled_first, input.first, result_second);
   for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
     for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
       const std::size_t index = row * kGlobalDofCount + column;
-      result.second[index] = first * input.second[index] +
-                             second * input.first[row] * input.first[column];
+      result_second[index] = first * input_second[index] + result_second[index];
     }
   }
   return result;
@@ -613,19 +744,38 @@ Jet Atan2(const Jet& y, const Jet& x) {
   const double derivative_yx =
       (y.value * y.value - x.value * x.value) / radius_fourth;
   const double derivative_xx = 2.0 * x.value * y.value / radius_fourth;
+  const bool zero_hessian = x.constant_derivatives && y.constant_derivatives &&
+      std::isfinite(result.value) && std::isfinite(derivative_y) &&
+      std::isfinite(derivative_x) && std::isfinite(derivative_yy) &&
+      std::isfinite(derivative_yx) && std::isfinite(derivative_xx);
   for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
     result.first[index] =
         derivative_y * y.first[index] + derivative_x * x.first[index];
   }
-  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];
+  if (!zero_hessian) {
+    std::array<double, kGlobalDofCount> scaled_y{};
+    std::array<double, kGlobalDofCount> scaled_x{};
+    for (std::size_t i = 0U; i < kGlobalDofCount; ++i) {
+      scaled_y[i] = derivative_yy * y.first[i];
+      scaled_x[i] = derivative_xx * x.first[i];
+    }
+    std::array<std::array<double, kHessianSize>, 4U> products;
+    Outer(scaled_y, y.first, products[0U].data());
+    Outer(y.first, x.first, products[1U].data());
+    Outer(x.first, y.first, products[2U].data());
+    Outer(scaled_x, x.first, products[3U].data());
+    const double* y_second = y.second.ReadData();
+    const double* x_second = x.second.ReadData();
+    double* result_second = result.second.OverwriteData();
+    for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
+      for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
+        const std::size_t index = row * kGlobalDofCount + column;
+        result_second[index] =
+            derivative_y * y_second[index] + derivative_x * x_second[index] +
+            products[0U][index] +
+            derivative_yx * (products[1U][index] + products[2U][index]) +
+            products[3U][index];
+      }
     }
   }
   // Use a bounded ratio for the directional third derivative only. Retain
@@ -644,6 +794,7 @@ Jet Atan2(const Jet& y, const Jet& x) {
     result.mixed_first[i] = sign * angle.mixed_first[i];
     result.mixed_second[i] = sign * angle.mixed_second[i];
   }
+  result.constant_derivatives = zero_hessian && angle.constant_derivatives;
   return result;
 }
 
@@ -1181,10 +1332,14 @@ Jet DrillingEnergy(const std::array<Vector3, kNodeCount>& coordinates,
           DrillSpinRate(trial.rotation[node], initial_directors[node]) * gamma;
     }
     energy = energy + 0.5 * drilling_stiffness * gamma * gamma;
+    std::array<double, kGlobalDofCount> scaled_gradient{};
+    for (std::size_t row = 0U; row < kGlobalDofCount; ++row)
+      scaled_gradient[row] = drilling_stiffness * gamma.first[row];
+    std::array<double, kHessianSize> contribution;
+    Outer(scaled_gradient, gamma.first, contribution.data());
     for (std::size_t row = 0U; row < kGlobalDofCount; ++row) {
       for (std::size_t column = 0U; column < kGlobalDofCount; ++column) {
-        (*metric)(row, column) +=
-            drilling_stiffness * gamma.first[row] * gamma.first[column];
+        (*metric)(row, column) += contribution[row * kGlobalDofCount + column];
       }
     }
   }
@@ -1505,6 +1660,7 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint(
     const std::size_t node, const Vector3& phi, const Vector3& phi_dot,
     const Vector3& force, const Vector3& force_rate,
     const Vector3& force_acceleration) const {
+  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
   if (node >= kNodeCount || !phi_dot.IsFinite() || !force.IsFinite() ||
       !force_rate.IsFinite() || !force_acceleration.IsFinite()) {
     return Failure<NonlinearMitc4ExternalDrillingConstraint>(
@@ -1594,6 +1750,9 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint(
 Result<NonlinearMitc4InternalEquilibrium> EvaluateInternalResidualAndTangent(
     const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial,
     const Vector* coordinate_direction) {
+  // Element-local MKL work uses one worker only when the caller has no TLS
+  // override. The guard restores that request before returning to global solve.
+  const dense_blas_internal::ScopedLocalDenseThreads local_threads{};
   if (coordinate_direction != nullptr) {
     if (coordinate_direction->Size() != kGlobalDofCount) {
       return Failure<NonlinearMitc4InternalEquilibrium>(