fix(elements): compensate Hex8 gradient and nonlinear force contractions

src/fesa/elements/hex8_element.cpp

@@ -26,6 +26,17 @@ std::array<double, 6> Pack(const Tensor& value) {
 double Dot(const Triple& a, const Triple& b) {
   return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
 }
+/// @brief Retains product and cancellation errors in fixed-order contractions.
+void AccumulateProduct(double left, double right, double& sum,
+                       double& correction) {
+  const double product = left * right;
+  const double next = sum + product;
+  correction += std::fma(left, right, -product) +
+                (std::abs(sum) >= std::abs(product)
+                     ? (sum - next) + product
+                     : (product - next) + sum);
+  sum = next;
+}
 /// @brief Gives the finite-control failure without accepting a geometry guess.
 const char* CertificateMessage(hex8_internal::Certificate certificate) {
   switch (certificate) {
@@ -248,26 +259,25 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
     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)
-        const double product = stiffness_(i, j) * relative[j];
+        AccumulateProduct(stiffness_(i, j), relative[j], sum, correction);
-        const double next = sum + product;
-        correction += std::fma(stiffness_(i, j), relative[j], -product) +
-                      (std::abs(sum) >= std::abs(product)
-                           ? (sum - next) + product
-                           : (product - next) + sum);
-        sum = next;
-      }
       response.internal_force[i] = sum + correction;
     }
     response.strain_energy = 0.5 * relative.Dot(response.internal_force);
   }
+  Vector force_correction{24U};
   for (std::size_t g = 0; g < 8; ++g) {
     const auto& point = points_[g];
     Tensor gradient{};
+    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)
-          gradient[i][j] += (u[3U * a + i] - u[i]) * point.gradients[a][j];
+          AccumulateProduct(u[3U * a + i] - u[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];
     Tensor f = gradient;
     for (std::size_t i = 0; i < 3; ++i) f[i][i] += 1.0;
     Tensor strain{};
@@ -309,8 +319,19 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
             stress[i][j] += p[i][k] * f[j][k] / det_f;
       for (std::size_t a = 0; a < 8; ++a) {
         const auto& ga = point.gradients[a];
-        for (std::size_t i = 0; i < 3; ++i)
+        for (std::size_t i = 0; i < 3; ++i) {
-          response.internal_force[3U * a + i] += Dot(p[i], ga) * point.weight;
+          double contraction = 0.0;
+          double contraction_correction = 0.0;
+          for (std::size_t j = 0; j < 3; ++j)
+            AccumulateProduct(p[i][j], ga[j], contraction,
+                              contraction_correction);
+          AccumulateProduct(contraction, point.weight,
+                            response.internal_force[3U * a + i],
+                            force_correction[3U * a + i]);
+          AccumulateProduct(contraction_correction, point.weight,
+                            response.internal_force[3U * a + i],
+                            force_correction[3U * a + i]);
+        }
         for (std::size_t b = 0; b < 8; ++b) {
           const auto& gb = point.gradients[b];
           double geometric = 0.0;
@@ -339,6 +360,9 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
                           point.natural,  point.position,
                           Pack(stress),   Pack(strain)};
   }
+  if (nonlinear)
+    for (std::size_t i = 0; i < 24; ++i)
+      response.internal_force[i] += force_correction[i];
   if (!std::isfinite(response.strain_energy))
     return Result<Hex8Response>::Failure(
         Failure("Hex8 trial energy is nonfinite."));