From 8293d5652b90404b0e7e04fc2c8d236964f66d53 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Fri, 16 Dec 2022 19:21:07 -0500
Subject: [PATCH] BarycentricPolynomials: clean up implementation.

---
 source/base/polynomials_barycentric.cc | 38 +++++++++++++-------------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/source/base/polynomials_barycentric.cc b/source/base/polynomials_barycentric.cc
index 4c899e2ee2..eb5601ccb8 100644
--- a/source/base/polynomials_barycentric.cc
+++ b/source/base/polynomials_barycentric.cc
@@ -13,9 +13,10 @@
 //
 // ---------------------------------------------------------------------
 
-
 #include <deal.II/base/polynomials_barycentric.h>
 
+#include <deal.II/grid/reference_cell.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace internal
@@ -60,35 +61,34 @@ BarycentricPolynomials<dim>::get_fe_p_basis(const unsigned int degree)
 {
   std::vector<PolyType> polys;
 
-  auto M = [](const unsigned int d) {
-    return BarycentricPolynomial<dim, double>::monomial(d);
-  };
+  const auto reference_cell = ReferenceCells::get_simplex<dim>();
+
   switch (degree)
     {
       case 0:
-        polys.push_back(0 * M(0) + 1);
+        polys.push_back(0 * BarycentricPolynomial<dim, double>::monomial(0) +
+                        1);
         break;
       case 1:
         {
-          for (unsigned int d = 0; d < dim + 1; ++d)
-            polys.push_back(M(d));
+          for (unsigned int v : reference_cell.vertex_indices())
+            polys.push_back(BarycentricPolynomial<dim, double>::monomial(v));
           break;
         }
       case 2:
         {
-          for (unsigned int d = 0; d < dim + 1; ++d)
-            polys.push_back(M(d) * (2 * M(d) - 1));
-          polys.push_back(4 * M(1) * M(0));
-          if (dim >= 2)
-            {
-              polys.push_back(4 * M(1) * M(2));
-              polys.push_back(4 * M(2) * M(0));
-            }
-          if (dim == 3)
+          // vertices, then lines:
+          for (unsigned int v : reference_cell.vertex_indices())
+            polys.push_back(
+              BarycentricPolynomial<dim, double>::monomial(v) *
+              (2 * BarycentricPolynomial<dim, double>::monomial(v) - 1));
+          for (unsigned int l : reference_cell.line_indices())
             {
-              polys.push_back(4 * M(3) * M(0));
-              polys.push_back(4 * M(1) * M(3));
-              polys.push_back(4 * M(2) * M(3));
+              const auto v0 = reference_cell.line_to_cell_vertices(l, 0);
+              const auto v1 = reference_cell.line_to_cell_vertices(l, 1);
+              polys.push_back(4 *
+                              BarycentricPolynomial<dim, double>::monomial(v0) *
+                              BarycentricPolynomial<dim, double>::monomial(v1));
             }
           break;
         }
-- 
2.39.5