fix(elements): retain low bits in Hex8 relative translations

src/fesa/elements/hex8_element.cpp

@@ -37,6 +37,12 @@ void AccumulateProduct(double left, double right, double& sum,
                      : (product - next) + sum);
   sum = next;
 }
+/// @brief Recovers low bits lost when subtracting unequal nodal translations.
+double DifferenceCorrection(double left, double right, double difference) {
+  const double right_virtual = left - difference;
+  const double left_virtual = difference + right_virtual;
+  return (left - left_virtual) + (right_virtual - right);
+}
 /// @brief Gives the finite-control failure without accepting a geometry guess.
 const char* CertificateMessage(hex8_internal::Certificate certificate) {
   switch (certificate) {
@@ -245,25 +251,41 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
     return Result<Hex8Response>::Failure(
         Failure(CertificateMessage(certificate)));
   Hex8Response response;
+  Vector relative{24U};
+  Vector relative_correction{24U};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i) {
+      const auto row = 3U * a + i;
+      relative[row] = u[row] - u[i];
+      relative_correction[row] =
+          DifferenceCorrection(u[row], u[i], relative[row]);
+    }
   if (!nonlinear) {
     response.material_tangent = stiffness_;
     response.tangent = stiffness_;
     // K0 annihilates uniform translation. Use the same relative field as
     // strain recovery to avoid cancellation of a large rigid displacement.
-    Vector relative{24U};
-    for (std::size_t a = 0; a < 8; ++a)
-      for (std::size_t i = 0; i < 3; ++i)
-        relative[3U * a + i] = u[3U * a + i] - u[i];
     // Retain low product bits and cancelled sum terms in bending-dominated
     // rows; the stored stiffness and its tangent are unchanged.
     for (std::size_t i = 0; i < 24; ++i) {
       double sum = 0.0;
       double correction = 0.0;
-      for (std::size_t j = 0; j < 24; ++j)
+      for (std::size_t j = 0; j < 24; ++j) {
         AccumulateProduct(stiffness_(i, j), relative[j], sum, correction);
+        AccumulateProduct(stiffness_(i, j), relative_correction[j], sum,
+                          correction);
+      }
       response.internal_force[i] = sum + correction;
     }
-    response.strain_energy = 0.5 * relative.Dot(response.internal_force);
+    double energy = 0.0;
+    double energy_correction = 0.0;
+    for (std::size_t i = 0; i < 24; ++i) {
+      AccumulateProduct(relative[i], response.internal_force[i], energy,
+                        energy_correction);
+      AccumulateProduct(relative_correction[i], response.internal_force[i],
+                        energy, energy_correction);
+    }
+    response.strain_energy = 0.5 * (energy + energy_correction);
   }
   Vector force_correction{24U};
   for (std::size_t g = 0; g < 8; ++g) {
@@ -272,9 +294,13 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
     Tensor gradient_correction{};
     for (std::size_t a = 0; a < 8; ++a)
       for (std::size_t i = 0; i < 3; ++i)
-        for (std::size_t j = 0; j < 3; ++j)
-          AccumulateProduct(u[3U * a + i] - u[i], point.gradients[a][j],
+        for (std::size_t j = 0; j < 3; ++j) {
+          AccumulateProduct(relative[3U * a + i], point.gradients[a][j],
                             gradient[i][j], gradient_correction[i][j]);
+          AccumulateProduct(relative_correction[3U * a + i],
+                            point.gradients[a][j],
+                            gradient[i][j], gradient_correction[i][j]);
+        }
     for (std::size_t i = 0; i < 3; ++i)
       for (std::size_t j = 0; j < 3; ++j)
         gradient[i][j] += gradient_correction[i][j];