perf: reuse owned shell Jet Hessians in accumulation

src/fesa/elements/nonlinear_mitc4_shell.cpp

@@ -386,6 +386,9 @@ class Hessian {
   /// @brief Identifies only the immutable zero Hessian, independently of seeds.
   bool IsImplicitZero() const { return !values_; }
 
+  /// @brief Allows a consumed temporary to reuse only its own materialized data.
+  bool HasUniqueStorage() const { return values_ && values_.use_count() == 1; }
+
   /// @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.
@@ -488,6 +491,30 @@ Jet operator+(const Jet& left, const Jet& right) {
   return result;
 }
 
+/// @brief Adds to a consumed temporary without allocating another Hessian.
+/// @note Shared/live aliases and unmasked exceptions retain the original path.
+Jet operator+(Jet&& left, const Jet& right) {
+  if ((_mm_getcsr() & _MM_MASK_MASK) != _MM_MASK_MASK || &left == &right ||
+      !left.second.HasUniqueStorage())
+    return static_cast<const Jet&>(left) + right;
+  left.value = left.value + right.value;
+  left.directional_first = left.directional_first + right.directional_first;
+  left.directional_second = left.directional_second + right.directional_second;
+  for (std::size_t index = 0U; index < kGlobalDofCount; ++index) {
+    left.first[index] = left.first[index] + right.first[index];
+    left.mixed_first[index] = left.mixed_first[index] + right.mixed_first[index];
+    left.mixed_second[index] = left.mixed_second[index] + right.mixed_second[index];
+  }
+  left.constant_derivatives = left.constant_derivatives &&
+                              right.constant_derivatives &&
+                              std::isfinite(left.value);
+  const double* right_second = right.second.ReadData();
+  double* left_second = left.second.Data();
+  for (std::size_t index = 0U; index < kHessianSize; ++index)
+    left_second[index] = left_second[index] + right_second[index];
+  return std::move(left);
+}
+
 /// @brief Subtracts two second-order trial-coordinate quantities.
 Jet operator-(const Jet& left, const Jet& right) {
   Jet result{};
@@ -1105,7 +1132,7 @@ Jet FrameComponent(const JetVector3& left, const JetTensor& tensor,
   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];
+      result = std::move(result) + left[row] * tensor[row][column] * right[column];
     }
   }
   return result;
@@ -1183,9 +1210,9 @@ Jet PhysicalEnergy(const std::array<Vector3, kNodeCount>& coordinates,
           shear * strain[4U]};
       Jet density{};
       for (std::size_t component = 0U; component < strain.size(); ++component) {
-        density = density + strain[component] * stress[component];
+        density = std::move(density) + strain[component] * stress[component];
       }
-      energy = energy + 0.5 * reference.reference_jacobian * density;
+      energy = std::move(energy) + 0.5 * reference.reference_jacobian * density;
     }
   }
   if (!IsFinite(energy)) {

src/fesa/math/jet_arithmetic.h

@@ -3,6 +3,7 @@
 
 #include <array>
 #include <cstddef>
+#include <utility>
 
 namespace fesa::math_internal {
 /// @brief Multiplies fixed Jet matrices in row/column/index order.
@@ -17,7 +18,7 @@ std::array<std::array<Jet, 3>, 3> MultiplyJetMatrix(
     for (std::size_t column = 0; column < 3; ++column) {
       result[row][column] = constant(0.0);
       for (std::size_t index = 0; index < 3; ++index)
-        result[row][column] = result[row][column] +
+        result[row][column] = std::move(result[row][column]) +
                               left[row][index] * right[index][column];
     }
   return result;
@@ -31,7 +32,7 @@ std::array<Jet, 3> MultiplyJetMatrixVector(
   std::array<Jet, 3> result{constant(0.0),constant(0.0),constant(0.0)};
   for (std::size_t row = 0; row < 3; ++row)
     for (std::size_t column = 0; column < 3; ++column)
-      result[row] = result[row] + matrix[row][column] * vector[column];
+      result[row] = std::move(result[row]) + matrix[row][column] * vector[column];
   return result;
 }
 }  // namespace fesa::math_internal