From 8f1c5702cd3091c36209039e8fdd443c29bede1f Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 5 Dec 2020 12:13:23 +0100
Subject: [PATCH] ScalarWedgePolynomial: Make polynomials for triangle and line
 class field

---
 include/deal.II/simplex/polynomials.h | 11 ++++++++
 source/simplex/polynomials.cc         | 40 ++++++++++++---------------
 2 files changed, 29 insertions(+), 22 deletions(-)

diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h
index a1c484aaa5..0fc13fdcdd 100644
--- a/include/deal.II/simplex/polynomials.h
+++ b/include/deal.II/simplex/polynomials.h
@@ -262,6 +262,17 @@ namespace Simplex
      */
     virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
     clone() const override;
+
+  private:
+    /**
+     * Scalar polynomial defined on a triangle.
+     */
+    const ScalarPolynomial<2> poly_tri;
+
+    /**
+     * Scalar polynomial defined on a line.
+     */
+    const ScalarPolynomial<1> poly_line;
   };
 
 
diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc
index bae98e5f95..1682c3d4f5 100644
--- a/source/simplex/polynomials.cc
+++ b/source/simplex/polynomials.cc
@@ -508,20 +508,26 @@ namespace Simplex
   ScalarWedgePolynomial<dim>::ScalarWedgePolynomial(const unsigned int degree)
     : ScalarPolynomialsBase<dim>(degree,
                                  compute_n_polynomials_wedge(dim, degree))
+    , poly_tri(degree)
+    , poly_line(degree)
   {}
 
 
   namespace
   {
     /**
-     * TODO
+     * Decompose the shape-function index of a linear wedge into an index
+     * to access the right shape function within the triangle and and within
+     * the line.
      */
     static const constexpr std::array<std::array<unsigned int, 2>, 6>
       wedge_table_1{
         {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
 
     /**
-     * TODO
+     * Decompose the shape-function index of a quadratic wedge into an index
+     * to access the right shape function within the triangle and and within
+     * the line.
      */
     static const constexpr std::array<std::array<unsigned int, 2>, 18>
       wedge_table_2{{{{0, 0}},
@@ -550,18 +556,13 @@ namespace Simplex
   ScalarWedgePolynomial<dim>::compute_value(const unsigned int i,
                                             const Point<dim> & p) const
   {
-    AssertDimension(dim, 3);
-    AssertIndexRange(this->degree(), 3);
-
     const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
 
-    const ScalarPolynomial<2> poly_tri(this->degree());
-    const Point<2>            p_tri(p[0], p[1]);
-    const auto                v_tri = poly_tri.compute_value(pair[0], p_tri);
+    const Point<2> p_tri(p[0], p[1]);
+    const auto     v_tri = poly_tri.compute_value(pair[0], p_tri);
 
-    const ScalarPolynomial<1> poly_line(this->degree());
-    const Point<1>            p_line(p[2]);
-    const auto                v_line = poly_line.compute_value(pair[1], p_line);
+    const Point<1> p_line(p[2]);
+    const auto     v_line = poly_line.compute_value(pair[1], p_line);
 
     return v_tri * v_line;
   }
@@ -573,20 +574,15 @@ namespace Simplex
   ScalarWedgePolynomial<dim>::compute_grad(const unsigned int i,
                                            const Point<dim> & p) const
   {
-    AssertDimension(dim, 3);
-    AssertIndexRange(this->degree(), 3);
-
     const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
 
-    const ScalarPolynomial<2> poly_tri(this->degree());
-    const Point<2>            p_tri(p[0], p[1]);
-    const auto                v_tri = poly_tri.compute_value(pair[0], p_tri);
-    const auto                g_tri = poly_tri.compute_grad(pair[0], p_tri);
+    const Point<2> p_tri(p[0], p[1]);
+    const auto     v_tri = poly_tri.compute_value(pair[0], p_tri);
+    const auto     g_tri = poly_tri.compute_grad(pair[0], p_tri);
 
-    const ScalarPolynomial<1> poly_line(this->degree());
-    const Point<1>            p_line(p[2]);
-    const auto                v_line = poly_line.compute_value(pair[1], p_line);
-    const auto                g_line = poly_line.compute_grad(pair[1], p_line);
+    const Point<1> p_line(p[2]);
+    const auto     v_line = poly_line.compute_value(pair[1], p_line);
+    const auto     g_line = poly_line.compute_grad(pair[1], p_line);
 
     Tensor<1, dim> grad;
     grad[0] = g_tri[0] * v_line;
-- 
2.39.5