perf(hex8): use MKL nodal mass and preserve rejected diagnostic priority

src/fesa/elements/hex8_element.cpp

@@ -180,6 +180,7 @@ Result<Hex8Element> Hex8Element::Create(EntityIndex index,
   // Zero Hessian of the same reference-relative mean-volume energy.
   std::array<Triple, 8> integrated_gradients{};
   double volume = 0.0;
+  bool backend_inputs_failed = false;
   for (const auto& point : candidate.points_) {
     auto gradients = point.gradients;
     gradients[0] = {};
@@ -192,22 +193,34 @@ Result<Hex8Element> Hex8Element::Create(EntityIndex index,
     Matrix deviatoric{0U, 0U};
     Matrix shear{0U, 0U};
     Matrix gram{0U, 0U};
+    bool point_products_failed = false;
     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(
+      backend_inputs_failed = true;
-          candidate.Failure("Hex8 reference stiffness is nonfinite."));
+      point_products_failed = true;
     }
     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) +=
+            if (point_products_failed) {
-                point.weight * (deviatoric(3U * a + i, 3U * b + k) +
+              // Diagnostic replay only: this candidate can never succeed.
-                                shear(3U * a + k, 3U * b + i) +
+              const auto& ga = gradients[a];
-                                (i == k ? candidate.mu_ * gram(a, b) : 0.0));
+              const auto& gb = gradients[b];
+              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(ga, gb) : 0.0));
+            } else {
+              candidate.stiffness_(3U * a + i, 3U * b + k) +=
+                  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));
+            }
           }
       }
   }
@@ -218,8 +231,15 @@ Result<Hex8Element> Hex8Element::Create(EntityIndex index,
   try {
     bulk = Outer(integrated_gradients, candidate.kappa_ / volume);
   } catch (const std::invalid_argument&) {
-    return Result<Hex8Element>::Failure(
+    backend_inputs_failed = true;
-        candidate.Failure("Hex8 reference stiffness is nonfinite."));
+    bulk = Matrix{24U, 24U};
+    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)
+            bulk(3U * a + i, 3U * b + k) =
+                (candidate.kappa_ / volume) * integrated_gradients[a][i] *
+                integrated_gradients[b][k];
   }
   for (std::size_t i = 0; i < 24; ++i)
     for (std::size_t j = 0; j < 24; ++j)
@@ -234,6 +254,9 @@ Result<Hex8Element> Hex8Element::Create(EntityIndex index,
       return Result<Hex8Element>::Failure(
           candidate.Failure("Hex8 reference stiffness has an unrepresentable "
                             "positive diagonal."));
+  if (backend_inputs_failed)
+    return Result<Hex8Element>::Failure(candidate.Failure(
+        "Hex8 dense MKL operands failed finite validation."));
   const Vector zero{24U};
   const auto paired = candidate.EvaluateHighLow(zero, zero, false);
   if (paired.HasValue()) {
@@ -262,17 +285,30 @@ Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
     return Result<ElementStiffnessContribution>::Failure(
         Failure("Hex8 mass requires reference density."));
   Matrix consistent{24U, 24U};
-  for (const auto& point : points_)
+  for (const auto& point : points_) {
+    Matrix left{8U, 1U};
+    Matrix right{1U, 8U};
+    for (std::size_t a = 0; a < 8; ++a) {
+      left(a, 0) = *density_ * point.shape[a];
+      right(0, a) = point.shape[a];
+    }
+    Matrix products{0U, 0U};
+    try {
+      products = left.Multiply(right);
+    } catch (const std::invalid_argument&) {
+      return Result<ElementStiffnessContribution>::Failure(
+          Failure("Hex8 reference mass is nonfinite."));
+    }
     for (std::size_t a = 0; a < 8; ++a)
       for (std::size_t b = 0; b < 8; ++b) {
-        const double value =
+        const double value = products(a, b) * point.weight;
-            *density_ * point.shape[a] * point.shape[b] * point.weight;
         if (!std::isfinite(value))
           return Result<ElementStiffnessContribution>::Failure(
               Failure("Hex8 reference mass is nonfinite."));
         for (std::size_t i = 0; i < 3; ++i)
           consistent(3U * a + i, 3U * b + i) += value;
       }
+  }
   Matrix mass{24U, 24U};
   for (std::size_t i = 0; i < 24; ++i)
     for (std::size_t j = 0; j < 24; ++j) {
@@ -387,6 +423,7 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
         Failure("Hex8 mean volume is nonfinite."));
   response.strain_energy = (0.5 * kappa_ * mean_trace) * (volume * mean_trace);
   std::array<Triple, 8> integrated_fg{};
+  bool backend_inputs_failed = false;
   for (std::size_t g = 0; g < 8; ++g) {
     const auto& point = points_[g];
     const auto& f = deformations[g];
@@ -428,6 +465,7 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
     Matrix shear{0U, 0U};
     Matrix gram{0U, 0U};
     Matrix geometric{0U, 0U};
+    bool point_products_failed = false;
     if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(p) ||
         !Finite(fft) || !Finite(stress))
       return Result<Hex8Response>::Failure(
@@ -443,8 +481,8 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
         gram = GradientProduct(gradients[g]);
         geometric = GradientProduct(gradients[g], &s);
       } catch (const std::invalid_argument&) {
-        return Result<Hex8Response>::Failure(
+        backend_inputs_failed = true;
-            Failure("Hex8 trial tangent is nonfinite."));
+        point_products_failed = true;
       }
     }
     for (std::size_t a = 0; a < 8; ++a) {
@@ -464,18 +502,34 @@ 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 diagnostic_geometric = 0.0;
+          if (point_products_failed)
+            for (std::size_t j = 0; j < 3; ++j)
+              for (std::size_t l = 0; l < 3; ++l)
+                diagnostic_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) +=
+              if (point_products_failed) {
-                  point.weight *
+                // Reproduce original diagnostic priority, never acceptance.
-                  (deviatoric(row, column) +
+                response.material_tangent(row, column) +=
-                   shear(3U * b + i, 3U * a + k) +
+                    point.weight *
-                   mu_ * fft[i][k] * gram(a, b));
+                    (-(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));
+              } else {
+                response.material_tangent(row, column) +=
+                    point.weight *
+                    (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(a, b) * point.weight;
+                    (point_products_failed ? diagnostic_geometric
+                                           : geometric(a, b)) * point.weight;
             }
         }
     }
@@ -490,8 +544,14 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
     try {
       bulk = Outer(integrated_fg, kappa_ / volume);
     } catch (const std::invalid_argument&) {
-      return Result<Hex8Response>::Failure(
+      backend_inputs_failed = true;
-          Failure("Hex8 trial tangent is nonfinite."));
+      bulk = Matrix{24U, 24U};
+      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)
+              bulk(3U * a + i, 3U * b + k) =
+                  (kappa_ / volume) * integrated_fg[a][i] * integrated_fg[b][k];
     }
     for (std::size_t i = 0; i < 24; ++i)
       for (std::size_t j = 0; j < 24; ++j)
@@ -515,6 +575,9 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
             Failure("Hex8 trial tangent is nonfinite."));
     }
   }
+  if (backend_inputs_failed)
+    return Result<Hex8Response>::Failure(
+        Failure("Hex8 dense MKL operands failed finite validation."));
   return Result<Hex8Response>::Success(std::move(response));
 }
 }  // namespace fesa