perf(hex8): use MKL Hessian products and exact zero pair addition

src/fesa/elements/hex8_element.cpp

@@ -2,6 +2,7 @@
 
 #include <cmath>
 #include <set>
+#include <stdexcept>
 #include <utility>
 
 #include "hex8_geometry.h"
@@ -26,6 +27,36 @@ 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 Forms a scaled 24-component outer product through dense MKL.
+Matrix Outer(const std::array<Triple, 8>& values, double coefficient) {
+  Matrix left{24U, 1U};
+  Matrix right{1U, 24U};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i) {
+      left(3U * a + i, 0) = coefficient * values[a][i];
+      right(0, 3U * a + i) = values[a][i];
+    }
+  return left.Multiply(right);
+}
+/// @brief Forms node-gradient Gram and stress congruence through dense MKL.
+Matrix GradientProduct(const std::array<Triple, 8>& gradients,
+                       const Tensor* stress = nullptr) {
+  Matrix left{8U, 3U};
+  Matrix transpose{3U, 8U};
+  for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t i = 0; i < 3; ++i) {
+      left(a, i) = gradients[a][i];
+      transpose(i, a) = gradients[a][i];
+    }
+  if (stress != nullptr) {
+    Matrix stress_transpose{3U, 3U};
+    for (std::size_t i = 0; i < 3; ++i)
+      for (std::size_t j = 0; j < 3; ++j)
+        stress_transpose(i, j) = (*stress)[j][i];
+    left = left.Multiply(stress_transpose);
+  }
+  return left.Multiply(transpose);
+}
 /// @brief Retains product and cancellation errors in fixed-order contractions.
 void AccumulateProduct(double left, double right, double& sum,
                        double& correction) {
@@ -158,30 +189,41 @@ Result<Hex8Element> Hex8Element::Create(EntityIndex index,
     for (std::size_t a = 0; a < 8; ++a)
       for (std::size_t i = 0; i < 3; ++i)
         integrated_gradients[a][i] += point.weight * gradients[a][i];
+    Matrix deviatoric{0U, 0U};
+    Matrix shear{0U, 0U};
+    Matrix gram{0U, 0U};
+    try {
+      deviatoric = Outer(gradients, -(2.0 / 3.0) * candidate.mu_);
+      shear = Outer(gradients, candidate.mu_);
+      gram = GradientProduct(gradients);
+    } catch (const std::invalid_argument&) {
+      return Result<Hex8Element>::Failure(
+          candidate.Failure("Hex8 reference stiffness is nonfinite."));
+    }
     for (std::size_t a = 0; a < 8; ++a)
       for (std::size_t b = 0; b < 8; ++b) {
-        const auto& ga = gradients[a];
-        const auto& gb = gradients[b];
-        const double dot = Dot(ga, gb);
         for (std::size_t i = 0; i < 3; ++i)
           for (std::size_t k = 0; k < 3; ++k) {
             candidate.stiffness_(3U * a + i, 3U * b + k) +=
-                point.weight * (-(2.0 / 3.0) * candidate.mu_ * ga[i] * gb[k] +
-                                candidate.mu_ * ga[k] * gb[i] +
-                                (i == k ? candidate.mu_ * dot : 0.0));
+                point.weight * (deviatoric(3U * a + i, 3U * b + k) +
+                                shear(3U * a + k, 3U * b + i) +
+                                (i == k ? candidate.mu_ * gram(a, b) : 0.0));
           }
       }
   }
   if (!std::isfinite(volume) || volume <= 0.0)
     return Result<Hex8Element>::Failure(
         candidate.Failure("Hex8 reference volume is nonfinite."));
-  for (std::size_t a = 0; a < 8; ++a)
-    for (std::size_t b = 0; b < 8; ++b)
-      for (std::size_t i = 0; i < 3; ++i)
-        for (std::size_t k = 0; k < 3; ++k)
-          candidate.stiffness_(3U * a + i, 3U * b + k) +=
-              (candidate.kappa_ / volume) * integrated_gradients[a][i] *
-              integrated_gradients[b][k];
+  Matrix bulk{0U, 0U};
+  try {
+    bulk = Outer(integrated_gradients, candidate.kappa_ / volume);
+  } catch (const std::invalid_argument&) {
+    return Result<Hex8Element>::Failure(
+        candidate.Failure("Hex8 reference stiffness is nonfinite."));
+  }
+  for (std::size_t i = 0; i < 24; ++i)
+    for (std::size_t j = 0; j < 24; ++j)
+      candidate.stiffness_(i, j) += bulk(i, j);
   for (std::size_t i = 0; i < 24; ++i)
     for (std::size_t j = 0; j < 24; ++j)
       if (!std::isfinite(candidate.stiffness_(i, j)))
@@ -382,6 +424,25 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
       for (std::size_t i = 0; i < 3; ++i)
         integrated_fg[a][i] += point.weight * Dot(f[i], gradients[g][a]);
     }
+    Matrix deviatoric{0U, 0U};
+    Matrix shear{0U, 0U};
+    Matrix gram{0U, 0U};
+    Matrix geometric{0U, 0U};
+    if (nonlinear) {
+      std::array<Triple, 8> fg{};
+      for (std::size_t a = 0; a < 8; ++a)
+        for (std::size_t i = 0; i < 3; ++i)
+          fg[a][i] = Dot(f[i], gradients[g][a]);
+      try {
+        deviatoric = Outer(fg, -(2.0 / 3.0) * mu_);
+        shear = Outer(fg, mu_);
+        gram = GradientProduct(gradients[g]);
+        geometric = GradientProduct(gradients[g], &s);
+      } catch (const std::invalid_argument&) {
+        return Result<Hex8Response>::Failure(
+            Failure("Hex8 trial tangent is nonfinite."));
+      }
+    }
     for (std::size_t a = 0; a < 8; ++a) {
       const auto& ga = gradients[g][a];
       for (std::size_t i = 0; i < 3; ++i) {
@@ -399,23 +460,18 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
       }
       if (nonlinear)
         for (std::size_t b = 0; b < 8; ++b) {
-          const auto& gb = gradients[g][b];
-          double geometric = 0.0;
-          for (std::size_t j = 0; j < 3; ++j)
-            for (std::size_t l = 0; l < 3; ++l)
-              geometric += ga[j] * s[l][j] * gb[l];
           for (std::size_t i = 0; i < 3; ++i)
             for (std::size_t k = 0; k < 3; ++k) {
               const auto row = 3U * a + i;
               const auto column = 3U * b + k;
               response.material_tangent(row, column) +=
                   point.weight *
-                  (-(2.0 / 3.0) * mu_ * Dot(f[i], ga) * Dot(f[k], gb) +
-                   mu_ * Dot(f[i], gb) * Dot(f[k], ga) +
-                   mu_ * fft[i][k] * Dot(ga, gb));
+                  (deviatoric(row, column) +
+                   shear(3U * b + i, 3U * a + k) +
+                   mu_ * fft[i][k] * gram(a, b));
               if (i == k)
                 response.geometric_tangent(row, column) +=
-                    geometric * point.weight;
+                    geometric(a, b) * point.weight;
             }
         }
     }
@@ -429,13 +485,18 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
   }
   for (std::size_t i = 0; i < 24; ++i)
     response.internal_force[i] += force_correction[i];
-  if (nonlinear)
-    for (std::size_t a = 0; a < 8; ++a)
-      for (std::size_t b = 0; b < 8; ++b)
-        for (std::size_t i = 0; i < 3; ++i)
-          for (std::size_t k = 0; k < 3; ++k)
-            response.material_tangent(3 * a + i, 3 * b + k) +=
-                (kappa_ / volume) * integrated_fg[a][i] * integrated_fg[b][k];
+  if (nonlinear) {
+    Matrix bulk{0U, 0U};
+    try {
+      bulk = Outer(integrated_fg, kappa_ / volume);
+    } catch (const std::invalid_argument&) {
+      return Result<Hex8Response>::Failure(
+          Failure("Hex8 trial tangent is nonfinite."));
+    }
+    for (std::size_t i = 0; i < 24; ++i)
+      for (std::size_t j = 0; j < 24; ++j)
+        response.material_tangent(i, j) += bulk(i, j);
+  }
   if (!std::isfinite(response.strain_energy))
     return Result<Hex8Response>::Failure(
         Failure("Hex8 trial energy is nonfinite."));