From 66b52ee3cb0c55a9b5b2ff6322c17383299e6f5b Mon Sep 17 00:00:00 2001
From: Alex Lindsay <alexlindsay239@gmail.com>
Date: Mon, 12 Sep 2022 10:47:50 -0700
Subject: [PATCH] Add partition-of-unity testing

This exercise was primarily undertaken to increase my understanding of our
implemented shape functions. I'm curious why in the literature (the literature
I've found at least) the partition-of-unity property (or lack thereof) is
relatively little discussed. It seems like this property (or lack thereof) has
very important consequences for being able to measure conservation ... or am I off here?
---
 tests/fe/fe_test.h | 76 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)

diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h
index 130e27bcf67..b9657c03a10 100644
--- a/tests/fe/fe_test.h
+++ b/tests/fe/fe_test.h
@@ -22,6 +22,7 @@
 
 #define FETEST                                  \
   CPPUNIT_TEST( testU );                        \
+  CPPUNIT_TEST( testPartitionOfUnity );         \
   CPPUNIT_TEST( testGradU );                    \
   CPPUNIT_TEST( testGradUComp );                \
   CPPUNIT_TEST( testHessU );                    \
@@ -564,6 +565,81 @@ class FETest : public FETestBase<order, family, elem_type, 1> {
 #endif // LIBMESH_ENABLE_EXCEPTIONS
   }
 
+  void testPartitionOfUnity()
+    {
+      if (!this->_elem)
+        return;
+
+      this->_fe->reinit(this->_elem);
+
+      bool satisfies_partition_of_unity = true;
+      for (const auto qp : make_range(this->_qrule->n_points()))
+      {
+        Real phi_sum = 0;
+        for (std::size_t d = 0; d != this->_dof_indices.size(); ++d)
+          phi_sum += this->_fe->get_phi()[d][qp];
+        if (phi_sum < (1 - TOLERANCE) || phi_sum > (1 + TOLERANCE))
+        {
+          satisfies_partition_of_unity = false;
+          break;
+        }
+      }
+
+      switch (this->_fe->get_family())
+      {
+        case MONOMIAL:
+        {
+          switch (this->_fe->get_order())
+          {
+            case CONSTANT:
+              CPPUNIT_ASSERT(satisfies_partition_of_unity);
+              break;
+
+            default:
+              CPPUNIT_ASSERT(!satisfies_partition_of_unity);
+              break;
+          }
+          break;
+        }
+        case SZABAB:
+        case HIERARCHIC:
+        case L2_HIERARCHIC:
+        {
+          switch (this->_fe->get_order())
+          {
+            case FIRST:
+              CPPUNIT_ASSERT(satisfies_partition_of_unity);
+              break;
+
+            default:
+              CPPUNIT_ASSERT(!satisfies_partition_of_unity);
+              break;
+          }
+          break;
+        }
+
+        case XYZ:
+        case CLOUGH:
+        case HERMITE:
+        {
+          CPPUNIT_ASSERT(!satisfies_partition_of_unity);
+          break;
+        }
+
+        case LAGRANGE:
+        case L2_LAGRANGE:
+        case BERNSTEIN:
+        case RATIONAL_BERNSTEIN:
+        {
+          CPPUNIT_ASSERT(satisfies_partition_of_unity);
+          break;
+        }
+
+        default:
+          CPPUNIT_FAIL("Uncovered FEFamily");
+      }
+    }
+
 
   void testU()
   {