fix(eigen): validate basis operands before dense dispatch

tests/unit/analysis/eigen_reduction_test.cpp

@@ -6,6 +6,7 @@
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
+#include <limits>
 #include <utility>
 #include <vector>
 
@@ -92,6 +93,38 @@ TEST(EigenReduction, RemovesOnlyActualFreeDirectorNull) {
                    .HasValue());
 }
 
+TEST(EigenReduction, RejectsNonfiniteUsedBasisAndPreservesPriorCandidate) {
+  const fesa::ShellNodeInitialFrame frame{0U, {0, 0, 1}, {1, 0, 0}, {0, 1, 0}};
+  const auto mass = NodalMass(frame.director);
+  fesa::Matrix stiffness(6U, 6U);
+  for (std::size_t i = 0; i < 5U; ++i) stiffness(i, i) = mass(i, i) * 7.0;
+  stiffness(5U, 5U) = .25;
+  const auto matrix = Sparse(stiffness);
+  const std::vector<std::uint8_t> mask(6U, 0U);
+  const auto prior = fesa::eigen_internal::BuildEigenReduction(
+      matrix, {mass}, {frame}, mask, {.25}, 2.0);
+  ASSERT_TRUE(prior.HasValue());
+  const auto prior_values = prior.Value().whitened_operator.Values();
+  for (const double invalid : {std::numeric_limits<double>::quiet_NaN(),
+                                std::numeric_limits<double>::infinity()}) {
+    for (const bool first_tangent : {false, true}) {
+      auto bad_frame = frame;
+      (first_tangent ? bad_frame.tangent_a : bad_frame.tangent_b)[0] = invalid;
+      const auto rejected = fesa::eigen_internal::BuildEigenReduction(
+          matrix, {mass}, {bad_frame}, mask, {.25}, 2.0);
+      ASSERT_FALSE(rejected.HasValue());
+      EXPECT_EQ(rejected.GetStatus().Category(), fesa::FailureCategory::kSolver);
+      EXPECT_EQ(rejected.GetStatus().Diagnostics()[0].code,
+                "invalid-eigen-mass-space");
+      EXPECT_EQ(prior.Value().whitened_operator.Values(), prior_values);
+    }
+  }
+  const auto rerun = fesa::eigen_internal::BuildEigenReduction(
+      matrix, {mass}, {frame}, mask, {.25}, 2.0);
+  ASSERT_TRUE(rerun.HasValue());
+  EXPECT_EQ(rerun.Value().whitened_operator.Values(), prior_values);
+}
+
 TEST(EigenReduction, RejectsPositiveButUnresolvableMassWithoutDeletingDofs) {
   const double p = 1e-8;
   const std::array<double, 3> d{2 * p / (1 + p * p), 0,