libMesh::isfinite() subroutines

include/base/libmesh_common.h

@@ -212,6 +212,15 @@ template <typename T>
 inline bool libmesh_isinf(std::complex<T> a)
 { return (std::isinf(std::real(a)) || std::isinf(std::imag(a))); }
 
+// std::isfinite() doesn't support complex, but we really want an
+// isfinite to support Number, at least in our own namespace
+template <typename T>
+inline bool isfinite(std::complex<T> a)
+{
+  using std::isfinite;
+  return (isfinite(std::real(a)) && isfinite(std::imag(a)));
+}
+
 // Define the value type for unknowns in simulations.
 // This is either Real or Complex, depending on how
 // the library was configures

include/numerics/dense_matrix.h

@@ -1218,6 +1218,18 @@ T DenseMatrix<T>::transpose (const unsigned int i,
 // }
 
 
+template <typename T>
+bool isfinite (const DenseMatrix<T> & var)
+{
+  using std::isfinite;
+  using libMesh::isfinite; // for T==complex
+  const auto m = var.m(), n = var.n();
+  for (auto i : make_range(m))
+    for (auto j : make_range(n))
+      if (!isfinite(var(i,j)))
+        return false;
+  return true;
+}
 
 
 } // namespace libMesh

include/numerics/dense_submatrix.h

@@ -229,6 +229,20 @@ T & DenseSubMatrix<T>::operator () (const unsigned int i,
 }
 
 
+template <typename T>
+bool isfinite (const DenseSubMatrix<T> & var)
+{
+  using std::isfinite;
+  using libMesh::isfinite; // for T==complex
+  const auto m = var.m(), n = var.n();
+  for (auto i : make_range(m))
+    for (auto j : make_range(n))
+      if (!isfinite(var(i,j)))
+        return false;
+  return true;
+}
+
+
 } // namespace libMesh
 
 

include/numerics/dense_vector.h

@@ -721,6 +721,21 @@ void DenseVector<T>::get_principal_subvector (unsigned int sub_n,
     dest(i) = _val[i];
 }
 
+
+template <typename T>
+bool isfinite (const DenseVector<T> & var)
+{
+  using std::isfinite;
+  using libMesh::isfinite; // for T==complex
+  for (auto i : make_range(var.size()))
+    if (!isfinite(var(i)))
+      return false;
+  return true;
+}
+
+
+
+
 } // namespace libMesh
 
 #ifdef LIBMESH_HAVE_METAPHYSICL

include/numerics/type_tensor.h

@@ -1435,6 +1435,20 @@ void TypeTensor<T>::print(std::ostream & os) const
 #endif
 }
 
+
+template <typename T>
+bool isfinite (const TypeTensor<T> & var)
+{
+  using std::isfinite;
+  using libMesh::isfinite; // for T==complex
+  for (unsigned int i=0; i<LIBMESH_DIM; i++)
+    for (unsigned int j=0; j<LIBMESH_DIM; j++)
+      if (!isfinite(var(i,j)))
+        return false;
+  return true;
+}
+
+
 template <typename T, typename T2>
 inline
 TypeTensor<typename CompareTypes<T, T2>::supertype>

include/numerics/type_vector.h

@@ -1186,6 +1186,19 @@ void TypeVector<T>::print(std::ostream & os) const
 #endif
 }
 
+
+template <typename T>
+bool isfinite (const TypeVector<T> & var)
+{
+  using std::isfinite;
+  using libMesh::isfinite; // for T==complex
+  for (unsigned int i=0; i<LIBMESH_DIM; i++)
+    if (!isfinite(var(i)))
+      return false;
+  return true;
+}
+
+
 template <typename T>
 struct CompareTypes<TypeVector<T>, TypeVector<T>>
 {

tests/libmesh_cppunit.h

@@ -19,12 +19,15 @@
          CPPUNIT_ASSERT_DOUBLES_EQUAL(expected,actual,tolerance)
 #endif
 
-#ifdef LIBMESH_USE_COMPLEX_NUMBERS
-# define LIBMESH_ASSERT_NUMBERS_EQUAL(expected,actual,tolerance) \
+# define LIBMESH_ASSERT_COMPLEX_EQUAL(expected,actual,tolerance) \
   do { \
     LIBMESH_ASSERT_FP_EQUAL(libMesh::libmesh_real(expected),libMesh::libmesh_real(actual),tolerance); \
     LIBMESH_ASSERT_FP_EQUAL(libMesh::libmesh_imag(expected),libMesh::libmesh_imag(actual),tolerance); \
   } while (0)
+
+#ifdef LIBMESH_USE_COMPLEX_NUMBERS
+# define LIBMESH_ASSERT_NUMBERS_EQUAL(expected,actual,tolerance) \
+         LIBMESH_ASSERT_COMPLEX_EQUAL(expected,actual,tolerance)
 #else
 # define LIBMESH_ASSERT_NUMBERS_EQUAL(expected,actual,tolerance) \
          LIBMESH_ASSERT_FP_EQUAL(expected,actual,tolerance)

tests/numerics/dense_matrix_test.C

@@ -1,5 +1,6 @@
 // libmesh includes
 #include <libmesh/dense_matrix.h>
+#include <libmesh/dense_submatrix.h>
 #include <libmesh/dense_vector.h>
 
 #ifdef LIBMESH_HAVE_PETSC
@@ -8,6 +9,7 @@
 
 #include "libmesh_cppunit.h"
 
+#include <limits>
 
 using namespace libMesh;
 
@@ -20,6 +22,7 @@ public:
 
   LIBMESH_CPPUNIT_TEST_SUITE(DenseMatrixTest);
 
+  CPPUNIT_TEST(testIsFinite);
   CPPUNIT_TEST(testOuterProduct);
   CPPUNIT_TEST(testSVD);
   CPPUNIT_TEST(testEVDreal);
@@ -32,6 +35,33 @@ public:
 
 private:
 
+  void testIsFinite()
+  {
+    LOG_UNIT_TEST;
+
+    DenseMatrix<Real> a(2, 2, {1, 2,
+                               3, 4});
+    CPPUNIT_ASSERT(isfinite(a));
+
+    DenseVector<Real> diag = a.diagonal();
+    CPPUNIT_ASSERT(isfinite(diag));
+
+    DenseSubMatrix suba1(a,0,0,1,1);
+    CPPUNIT_ASSERT(isfinite(suba1));
+
+    DenseSubMatrix suba2(a,1,1,1,1);
+    CPPUNIT_ASSERT(isfinite(suba2));
+
+    a(0, 0) = std::numeric_limits<Real>::infinity();
+    CPPUNIT_ASSERT(!isfinite(a));
+    CPPUNIT_ASSERT(!isfinite(suba1));
+    CPPUNIT_ASSERT(isfinite(suba2));
+
+    CPPUNIT_ASSERT(isfinite(diag));
+    diag = a.diagonal();
+    CPPUNIT_ASSERT(!isfinite(diag));
+  }
+
   void testOuterProduct()
   {
     LOG_UNIT_TEST;

tests/numerics/type_tensor_test.C

@@ -23,6 +23,7 @@ public:
   CPPUNIT_TEST(testRotation);
 #endif
   CPPUNIT_TEST(testRowCol);
+  CPPUNIT_TEST(testIsFinite);
   CPPUNIT_TEST(testIsZero);
   CPPUNIT_TEST(testIsHPD);
 #ifdef LIBMESH_HAVE_METAPHYSICL
@@ -91,6 +92,21 @@ private:
     LIBMESH_ASSERT_FP_EQUAL(28, product(2, 2), tol);
   }
 
+  void testIsFinite()
+  {
+    LOG_UNIT_TEST;
+
+    {
+      TensorValue<double> tensor;
+      CPPUNIT_ASSERT(isfinite(tensor));
+    }
+    {
+      TensorValue<double> tensor(0,1,0,2,3);
+      tensor(0,0) = std::numeric_limits<double>::infinity();
+      CPPUNIT_ASSERT(!isfinite(tensor));
+    }
+  }
+
   void testIsZero()
   {
     LOG_UNIT_TEST;