From 6f38c52acb17c10d7b563deed04b38e014c12909 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Fri, 30 Jun 2023 18:10:42 -0400
Subject: [PATCH] Implement a cubic version of BarycentricPolynomials.

---
 source/base/polynomials_barycentric.cc | 51 ++++++++++++++++++++++
 tests/simplex/barycentric_01.cc        | 59 +++++++++++++++++++++++++-
 tests/simplex/barycentric_01.output    |  6 ++-
 3 files changed, 114 insertions(+), 2 deletions(-)

diff --git a/source/base/polynomials_barycentric.cc b/source/base/polynomials_barycentric.cc
index bbda850c11..7c6c086507 100644
--- a/source/base/polynomials_barycentric.cc
+++ b/source/base/polynomials_barycentric.cc
@@ -92,6 +92,57 @@ BarycentricPolynomials<dim>::get_fe_p_basis(const unsigned int degree)
             }
           break;
         }
+      case 3:
+        {
+          // vertices, then lines, then quads:
+          for (const unsigned int v : reference_cell.vertex_indices())
+            polys.push_back(
+              0.5 * BarycentricPolynomial<dim, double>::monomial(v) *
+              (3 * BarycentricPolynomial<dim, double>::monomial(v) - 1) *
+              (3 * BarycentricPolynomial<dim, double>::monomial(v) - 2));
+          for (unsigned int l : reference_cell.line_indices())
+            {
+              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.5 * BarycentricPolynomial<dim, double>::monomial(v0) *
+                (3 * BarycentricPolynomial<dim, double>::monomial(v0) - 1) *
+                BarycentricPolynomial<dim, double>::monomial(v1));
+              polys.push_back(
+                4.5 * BarycentricPolynomial<dim, double>::monomial(v0) *
+                (3 * BarycentricPolynomial<dim, double>::monomial(v1) - 1) *
+                BarycentricPolynomial<dim, double>::monomial(v1));
+            }
+
+          if (dim == 2)
+            {
+              polys.push_back(27 *
+                              BarycentricPolynomial<dim, double>::monomial(0) *
+                              BarycentricPolynomial<dim, double>::monomial(1) *
+                              BarycentricPolynomial<dim, double>::monomial(2));
+            }
+          else if (dim == 3)
+            {
+              polys.push_back(27 *
+                              BarycentricPolynomial<dim, double>::monomial(0) *
+                              BarycentricPolynomial<dim, double>::monomial(1) *
+                              BarycentricPolynomial<dim, double>::monomial(2));
+              polys.push_back(27 *
+                              BarycentricPolynomial<dim, double>::monomial(0) *
+                              BarycentricPolynomial<dim, double>::monomial(1) *
+                              BarycentricPolynomial<dim, double>::monomial(3));
+              polys.push_back(27 *
+                              BarycentricPolynomial<dim, double>::monomial(0) *
+                              BarycentricPolynomial<dim, double>::monomial(2) *
+                              BarycentricPolynomial<dim, double>::monomial(3));
+              polys.push_back(27 *
+                              BarycentricPolynomial<dim, double>::monomial(1) *
+                              BarycentricPolynomial<dim, double>::monomial(2) *
+                              BarycentricPolynomial<dim, double>::monomial(3));
+            }
+
+          break;
+        }
       default:
         Assert(false, ExcNotImplemented());
     }
diff --git a/tests/simplex/barycentric_01.cc b/tests/simplex/barycentric_01.cc
index 2d9e199117..b766fbf865 100644
--- a/tests/simplex/barycentric_01.cc
+++ b/tests/simplex/barycentric_01.cc
@@ -107,10 +107,39 @@ main()
           }
         deallog << std::endl;
       }
+    deallog << "Test with TRI6 - Success" << std::endl;
   }
 
   {
-    deallog << std::endl << "Test with TET4" << std::endl;
+    deallog << "Test with TRI10" << std::endl;
+    const auto tri10 = BarycentricPolynomials<2>::get_fe_p_basis(3);
+
+    FE_SimplexP<2> fe(3);
+    const auto    &points = fe.get_unit_support_points();
+    for (unsigned int i = 0; i < 10; ++i)
+      {
+        Assert(points.size() == 10, ExcInternalError());
+        for (unsigned int j = 0; j < 10; ++j)
+          {
+            Assert(std::abs(tri10.compute_value(i, points[j]) -
+                            double(i == j)) < 1e-12,
+                   ExcInternalError());
+
+            // third derivatives should be constant
+            Assert((tri10.compute_3rd_derivative(i, points[0]) -
+                    tri10.compute_3rd_derivative(i, points[j]))
+                       .norm() == 0.0,
+                   ExcInternalError());
+
+            Assert(tri10.compute_4th_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+          }
+      }
+    deallog << "Test with TRI10 - Success" << std::endl;
+  }
+
+  {
+    deallog << "Test with TET4" << std::endl;
     const auto tet4 = BarycentricPolynomials<3>::get_fe_p_basis(1);
 
     FE_SimplexP<3> fe(1);
@@ -169,4 +198,32 @@ main()
       }
     deallog << "Test with TET10 - Success" << std::endl;
   }
+
+  {
+    deallog << "Test with TET20" << std::endl;
+    const auto tet20 = BarycentricPolynomials<3>::get_fe_p_basis(3);
+
+    FE_SimplexP<3> fe(3);
+    const auto    &points = fe.get_unit_support_points();
+    for (unsigned int i = 0; i < 20; ++i)
+      {
+        Assert(points.size() == 20, ExcInternalError());
+        for (unsigned int j = 0; j < 20; ++j)
+          {
+            Assert(std::abs(tet20.compute_value(i, points[j]) -
+                            double(i == j)) < 1e-12,
+                   ExcInternalError());
+
+            // third derivatives should be constant
+            Assert((tet20.compute_3rd_derivative(i, points[0]) -
+                    tet20.compute_3rd_derivative(i, points[j]))
+                       .norm() == 0.0,
+                   ExcInternalError());
+
+            Assert(tet20.compute_4th_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+          }
+      }
+    deallog << "Test with TET20 - Success" << std::endl;
+  }
 }
diff --git a/tests/simplex/barycentric_01.output b/tests/simplex/barycentric_01.output
index f083f6c228..0cc6aaa4b3 100644
--- a/tests/simplex/barycentric_01.output
+++ b/tests/simplex/barycentric_01.output
@@ -44,8 +44,12 @@ DEAL::p = 4.00000 * t0^1 * t2^1
 DEAL::p_x = -4.00000 * t2^1
 DEAL::p_y = -4.00000 * t2^1 + 4.00000 * t0^1
 DEAL::
-DEAL::
+DEAL::Test with TRI6 - Success
+DEAL::Test with TRI10
+DEAL::Test with TRI10 - Success
 DEAL::Test with TET4
 DEAL::Test with TET4 - Success
 DEAL::Test with TET10
 DEAL::Test with TET10 - Success
+DEAL::Test with TET20
+DEAL::Test with TET20 - Success
-- 
2.39.5