From: grahambenharper <grahambenharper@gmail.com>
Date: Mon, 9 Sep 2019 04:55:33 +0000 (-0600)
Subject: Change TensorProductPolynomials and AnisotropicPolynomials to derive from ScalarPolyn... 
X-Git-Tag: v9.2.0-rc1~1103^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8723%2Fhead;p=dealii.git

Change TensorProductPolynomials and AnisotropicPolynomials to derive from ScalarPolynomialsBase
---

diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h
index 3fd35d075c..1022310303 100644
--- a/include/deal.II/base/tensor_product_polynomials.h
+++ b/include/deal.II/base/tensor_product_polynomials.h
@@ -22,6 +22,7 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/polynomial.h>
+#include <deal.II/base/scalar_polynomials_base.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/utilities.h>
 
@@ -65,7 +66,7 @@ class TensorProductPolynomialsBubbles;
  * 2003
  */
 template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
-class TensorProductPolynomials
+class TensorProductPolynomials : public ScalarPolynomialsBase<dim>
 {
 public:
   /**
@@ -126,7 +127,7 @@ public:
            std::vector<Tensor<1, dim>> &grads,
            std::vector<Tensor<2, dim>> &grad_grads,
            std::vector<Tensor<3, dim>> &third_derivatives,
-           std::vector<Tensor<4, dim>> &fourth_derivatives) const;
+           std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
 
   /**
    * Compute the value of the <tt>i</tt>th tensor product polynomial at
@@ -192,12 +193,16 @@ public:
   compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
 
   /**
-   * Return the number of tensor product polynomials. For <i>n</i> 1d
-   * polynomials this is <i>n<sup>dim</sup></i>.
+   * Return the name of the space, which is <tt>TensorProductPolynomials</tt>.
    */
-  unsigned int
-  n() const;
+  std::string
+  name() const override;
 
+  /**
+   * @copydoc ScalarPolynomialsBase<dim>::clone()
+   */
+  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+  clone() const override;
 
 protected:
   /**
@@ -205,11 +210,6 @@ protected:
    */
   std::vector<PolynomialType> polynomials;
 
-  /**
-   * Number of tensor product polynomials. See n().
-   */
-  unsigned int n_tensor_pols;
-
   /**
    * Index map for reordering the polynomials.
    */
@@ -268,7 +268,7 @@ protected:
  * @author Wolfgang Bangerth 2003
  */
 template <int dim>
-class AnisotropicPolynomials
+class AnisotropicPolynomials : public ScalarPolynomialsBase<dim>
 {
 public:
   /**
@@ -282,6 +282,9 @@ public:
    * Since we want to build <i>anisotropic</i> polynomials, the @p dim
    * sets of polynomials passed in as arguments may of course be
    * different, and may also vary in number.
+   *
+   * The number of tensor product polynomials is <tt>Nx*Ny*Nz</tt>, or with
+   * terms dropped if the number of space dimensions is less than 3.
    */
   AnisotropicPolynomials(
     const std::vector<std::vector<Polynomials::Polynomial<double>>>
@@ -292,7 +295,7 @@ public:
    * product polynomial at <tt>unit_point</tt>.
    *
    * The size of the vectors must either be equal <tt>0</tt> or equal
-   * <tt>n_tensor_pols</tt>.  In the first case, the function will not compute
+   * <tt>this->n()</tt>.  In the first case, the function will not compute
    * these values.
    *
    * If you need values or derivatives of all tensor product polynomials then
@@ -301,12 +304,12 @@ public:
    * in a loop over all tensor product polynomials.
    */
   void
-  compute(const Point<dim> &           unit_point,
-          std::vector<double> &        values,
-          std::vector<Tensor<1, dim>> &grads,
-          std::vector<Tensor<2, dim>> &grad_grads,
-          std::vector<Tensor<3, dim>> &third_derivatives,
-          std::vector<Tensor<4, dim>> &fourth_derivatives) const;
+  evaluate(const Point<dim> &           unit_point,
+           std::vector<double> &        values,
+           std::vector<Tensor<1, dim>> &grads,
+           std::vector<Tensor<2, dim>> &grad_grads,
+           std::vector<Tensor<3, dim>> &third_derivatives,
+           std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
 
   /**
    * Compute the value of the <tt>i</tt>th tensor product polynomial at
@@ -317,7 +320,7 @@ public:
    * polynomials is not efficient, because then each point value of the
    * underlying (one-dimensional) polynomials is (unnecessarily) computed
    * several times.  Instead use the <tt>compute</tt> function, see above,
-   * with <tt>values.size()==n_tensor_pols</tt> to get the point values of all
+   * with <tt>values.size()==this->n()</tt> to get the point values of all
    * tensor polynomials all at once and in a much more efficient way.
    */
   double
@@ -350,7 +353,7 @@ public:
    * polynomials is not efficient, because then each derivative value of the
    * underlying (one-dimensional) polynomials is (unnecessarily) computed
    * several times.  Instead use the <tt>compute</tt> function, see above,
-   * with <tt>grads.size()==n_tensor_pols</tt> to get the point value of all
+   * with <tt>grads.size()==this->n()</tt> to get the point value of all
    * tensor polynomials all at once and in a much more efficient way.
    */
   Tensor<1, dim>
@@ -365,30 +368,29 @@ public:
    * polynomials is not efficient, because then each derivative value of the
    * underlying (one-dimensional) polynomials is (unnecessarily) computed
    * several times.  Instead use the <tt>compute</tt> function, see above,
-   * with <tt>grad_grads.size()==n_tensor_pols</tt> to get the point value of
+   * with <tt>grad_grads.size()==this->n()</tt> to get the point value of
    * all tensor polynomials all at once and in a much more efficient way.
    */
   Tensor<2, dim>
   compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
 
   /**
-   * Return the number of tensor product polynomials. It is the product of
-   * the number of polynomials in each coordinate direction.
+   * Return the name of the space, which is <tt>AnisotropicPolynomials</tt>.
    */
-  unsigned int
-  n() const;
+  std::string
+  name() const override;
 
-private:
   /**
-   * Copy of the vector <tt>pols</tt> of polynomials given to the constructor.
+   * @copydoc ScalarPolynomialsBase<dim>::clone()
    */
-  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
+  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+  clone() const override;
 
+private:
   /**
-   * Number of tensor product polynomials. This is <tt>Nx*Ny*Nz</tt>, or with
-   * terms dropped if the number of space dimensions is less than 3.
+   * Copy of the vector <tt>pols</tt> of polynomials given to the constructor.
    */
-  const unsigned int n_tensor_pols;
+  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
 
   /**
    * Each tensor product polynomial @þ{i} is a product of one-dimensional
@@ -400,7 +402,7 @@ private:
   compute_index(const unsigned int i, unsigned int (&indices)[dim]) const;
 
   /**
-   * Given the input to the constructor, compute <tt>n_tensor_pols</tt>.
+   * Given the input to the constructor, compute <tt>n_pols</tt>.
    */
   static unsigned int
   get_n_tensor_pols(
@@ -419,14 +421,14 @@ template <int dim, typename PolynomialType>
 template <class Pol>
 inline TensorProductPolynomials<dim, PolynomialType>::TensorProductPolynomials(
   const std::vector<Pol> &pols)
-  : polynomials(pols.begin(), pols.end())
-  , n_tensor_pols(Utilities::fixed_power<dim>(pols.size()))
-  , index_map(n_tensor_pols)
-  , index_map_inverse(n_tensor_pols)
+  : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
+  , polynomials(pols.begin(), pols.end())
+  , index_map(this->n())
+  , index_map_inverse(this->n())
 {
   // per default set this index map to identity. This map can be changed by
   // the user through the set_numbering() function
-  for (unsigned int i = 0; i < n_tensor_pols; ++i)
+  for (unsigned int i = 0; i < this->n(); ++i)
     {
       index_map[i]         = i;
       index_map_inverse[i] = i;
@@ -434,19 +436,6 @@ inline TensorProductPolynomials<dim, PolynomialType>::TensorProductPolynomials(
 }
 
 
-
-template <int dim, typename PolynomialType>
-inline unsigned int
-TensorProductPolynomials<dim, PolynomialType>::n() const
-{
-  if (dim == 0)
-    return numbers::invalid_unsigned_int;
-  else
-    return n_tensor_pols;
-}
-
-
-
 template <int dim, typename PolynomialType>
 inline const std::vector<unsigned int> &
 TensorProductPolynomials<dim, PolynomialType>::get_numbering() const
@@ -462,6 +451,15 @@ TensorProductPolynomials<dim, PolynomialType>::get_numbering_inverse() const
   return index_map_inverse;
 }
 
+
+template <int dim, typename PolynomialType>
+inline std::string
+TensorProductPolynomials<dim, PolynomialType>::name() const
+{
+  return "TensorProductPolynomials";
+}
+
+
 template <int dim, typename PolynomialType>
 template <int order>
 Tensor<order, dim>
@@ -710,6 +708,14 @@ AnisotropicPolynomials<dim>::compute_derivative(const unsigned int i,
 }
 
 
+template <int dim>
+inline std::string
+AnisotropicPolynomials<dim>::name() const
+{
+  return "AnisotropicPolynomials";
+}
+
+
 
 #endif // DOXYGEN
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h
index af66095a6d..77fddfaa41 100644
--- a/include/deal.II/base/tensor_product_polynomials_bubbles.h
+++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h
@@ -44,7 +44,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Daniel Arndt, 2015
  */
 template <int dim>
-class TensorProductPolynomialsBubbles
+class TensorProductPolynomialsBubbles : public ScalarPolynomialsBase<dim>
 {
 public:
   /**
@@ -105,7 +105,7 @@ public:
            std::vector<Tensor<1, dim>> &grads,
            std::vector<Tensor<2, dim>> &grad_grads,
            std::vector<Tensor<3, dim>> &third_derivatives,
-           std::vector<Tensor<4, dim>> &fourth_derivatives) const;
+           std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
 
   /**
    * Compute the value of the <tt>i</tt>th tensor product polynomial at
@@ -177,6 +177,18 @@ public:
   unsigned int
   n() const;
 
+  /**
+   * Return the name of the space, which is <tt>AnisotropicPolynomials</tt>.
+   */
+  std::string
+  name() const override;
+
+  /**
+   * @copydoc ScalarPolynomialsBase<dim>::clone()
+   */
+  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+  clone() const override;
+
 private:
   /**
    * The TensorProductPolynomials object
@@ -205,7 +217,9 @@ template <int dim>
 template <class Pol>
 inline TensorProductPolynomialsBubbles<dim>::TensorProductPolynomialsBubbles(
   const std::vector<Pol> &pols)
-  : tensor_polys(pols)
+  : ScalarPolynomialsBase<dim>(1,
+                               Utilities::fixed_power<dim>(pols.size()) + dim)
+  , tensor_polys(pols)
   , index_map(tensor_polys.n() +
               ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
   , index_map_inverse(tensor_polys.n() +
@@ -254,6 +268,14 @@ TensorProductPolynomialsBubbles<dim>::get_numbering_inverse() const
 }
 
 
+template <int dim>
+inline std::string
+TensorProductPolynomialsBubbles<dim>::name() const
+{
+  return "TensorProductPolynomialsBubbles";
+}
+
+
 template <int dim>
 template <int order>
 Tensor<order, dim>
@@ -264,7 +286,6 @@ TensorProductPolynomialsBubbles<dim>::compute_derivative(
   const unsigned int q_degree      = tensor_polys.polynomials.size() - 1;
   const unsigned int max_q_indices = tensor_polys.n();
   const unsigned int n_bubbles     = ((q_degree <= 1) ? 1 : dim);
-  (void)n_bubbles;
   Assert(i < max_q_indices + n_bubbles, ExcInternalError());
 
   // treat the regular basis functions
diff --git a/include/deal.II/base/tensor_product_polynomials_const.h b/include/deal.II/base/tensor_product_polynomials_const.h
index 1671532506..a9b1e131f8 100644
--- a/include/deal.II/base/tensor_product_polynomials_const.h
+++ b/include/deal.II/base/tensor_product_polynomials_const.h
@@ -80,7 +80,7 @@ public:
            std::vector<Tensor<1, dim>> &grads,
            std::vector<Tensor<2, dim>> &grad_grads,
            std::vector<Tensor<3, dim>> &third_derivatives,
-           std::vector<Tensor<4, dim>> &fourth_derivatives) const;
+           std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
 
   /**
    * Compute the value of the <tt>i</tt>th tensor product polynomial at
@@ -167,8 +167,8 @@ inline TensorProductPolynomialsConst<dim>::TensorProductPolynomialsConst(
   : TensorProductPolynomials<dim>(pols)
 {
   // append index for renumbering
-  this->index_map.push_back(this->n_tensor_pols);
-  this->index_map_inverse.push_back(this->n_tensor_pols);
+  this->index_map.push_back(TensorProductPolynomials<dim>::n());
+  this->index_map_inverse.push_back(TensorProductPolynomials<dim>::n());
 }
 
 
@@ -177,7 +177,7 @@ template <int dim>
 inline unsigned int
 TensorProductPolynomialsConst<dim>::n() const
 {
-  return this->n_tensor_pols + 1;
+  return TensorProductPolynomials<dim>::n() + 1;
 }
 
 
@@ -196,7 +196,7 @@ TensorProductPolynomialsConst<dim>::compute_derivative(
   const unsigned int i,
   const Point<dim> & p) const
 {
-  const unsigned int max_indices = this->n_tensor_pols;
+  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
diff --git a/source/base/polynomials_abf.cc b/source/base/polynomials_abf.cc
index 7f737e5f56..05ad5480ec 100644
--- a/source/base/polynomials_abf.cc
+++ b/source/base/polynomials_abf.cc
@@ -112,12 +112,12 @@ PolynomialsABF<dim>::evaluate(
       for (unsigned int c = 0; c < dim; ++c)
         p(c) = unit_point((c + d) % dim);
 
-      polynomial_space.compute(p,
-                               p_values,
-                               p_grads,
-                               p_grad_grads,
-                               p_third_derivatives,
-                               p_fourth_derivatives);
+      polynomial_space.evaluate(p,
+                                p_values,
+                                p_grads,
+                                p_grad_grads,
+                                p_third_derivatives,
+                                p_fourth_derivatives);
 
       for (unsigned int i = 0; i < p_values.size(); ++i)
         values[i + d * n_sub][d] = p_values[i];
diff --git a/source/base/polynomials_bernardi_raugel.cc b/source/base/polynomials_bernardi_raugel.cc
index ddeab2a319..a42ac870e4 100644
--- a/source/base/polynomials_bernardi_raugel.cc
+++ b/source/base/polynomials_bernardi_raugel.cc
@@ -137,7 +137,7 @@ PolynomialsBernardiRaugel<dim>::evaluate(
     }
 
   // set indices for the anisotropic polynomials to find
-  // them after polynomial_space_bubble.compute is called
+  // them after polynomial_space_bubble.evaluate is called
   std::vector<int> aniso_indices;
   if (dim == 2)
     {
@@ -156,18 +156,18 @@ PolynomialsBernardiRaugel<dim>::evaluate(
       aniso_indices.push_back(17);
     }
 
-  polynomial_space_bubble.compute(unit_point,
-                                  bubble_values,
-                                  bubble_grads,
-                                  bubble_grad_grads,
-                                  bubble_third_derivatives,
-                                  bubble_fourth_derivatives);
-  polynomial_space_Q.compute(unit_point,
-                             Q_values,
-                             Q_grads,
-                             Q_grad_grads,
-                             Q_third_derivatives,
-                             Q_fourth_derivatives);
+  polynomial_space_bubble.evaluate(unit_point,
+                                   bubble_values,
+                                   bubble_grads,
+                                   bubble_grad_grads,
+                                   bubble_third_derivatives,
+                                   bubble_fourth_derivatives);
+  polynomial_space_Q.evaluate(unit_point,
+                              Q_values,
+                              Q_grads,
+                              Q_grad_grads,
+                              Q_third_derivatives,
+                              Q_fourth_derivatives);
 
   // first dim*vertices_per_cell functions are Q_1^2 functions
   for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
diff --git a/source/base/polynomials_nedelec.cc b/source/base/polynomials_nedelec.cc
index 143806310d..4151ee8557 100644
--- a/source/base/polynomials_nedelec.cc
+++ b/source/base/polynomials_nedelec.cc
@@ -95,12 +95,12 @@ PolynomialsNedelec<dim>::evaluate(
     {
       case 1:
         {
-          polynomial_space.compute(unit_point,
-                                   unit_point_values,
-                                   unit_point_grads,
-                                   unit_point_grad_grads,
-                                   empty_vector_of_3rd_order_tensors,
-                                   empty_vector_of_4th_order_tensors);
+          polynomial_space.evaluate(unit_point,
+                                    unit_point_values,
+                                    unit_point_grads,
+                                    unit_point_grad_grads,
+                                    empty_vector_of_3rd_order_tensors,
+                                    empty_vector_of_4th_order_tensors);
 
           // Assign the correct values to the
           // corresponding shape functions.
@@ -121,12 +121,12 @@ PolynomialsNedelec<dim>::evaluate(
 
       case 2:
         {
-          polynomial_space.compute(unit_point,
-                                   unit_point_values,
-                                   unit_point_grads,
-                                   unit_point_grad_grads,
-                                   empty_vector_of_3rd_order_tensors,
-                                   empty_vector_of_4th_order_tensors);
+          polynomial_space.evaluate(unit_point,
+                                    unit_point_values,
+                                    unit_point_grads,
+                                    unit_point_grad_grads,
+                                    empty_vector_of_3rd_order_tensors,
+                                    empty_vector_of_4th_order_tensors);
 
           // Declare the values, derivatives and
           // second derivatives vectors of
@@ -144,12 +144,12 @@ PolynomialsNedelec<dim>::evaluate(
           std::vector<Tensor<2, dim>> p_grad_grads(
             (grad_grads.size() == 0) ? 0 : n_basis);
 
-          polynomial_space.compute(p,
-                                   p_values,
-                                   p_grads,
-                                   p_grad_grads,
-                                   empty_vector_of_3rd_order_tensors,
-                                   empty_vector_of_4th_order_tensors);
+          polynomial_space.evaluate(p,
+                                    p_values,
+                                    p_grads,
+                                    p_grad_grads,
+                                    empty_vector_of_3rd_order_tensors,
+                                    empty_vector_of_4th_order_tensors);
 
           // Assign the correct values to the
           // corresponding shape functions.
@@ -307,12 +307,12 @@ PolynomialsNedelec<dim>::evaluate(
 
       case 3:
         {
-          polynomial_space.compute(unit_point,
-                                   unit_point_values,
-                                   unit_point_grads,
-                                   unit_point_grad_grads,
-                                   empty_vector_of_3rd_order_tensors,
-                                   empty_vector_of_4th_order_tensors);
+          polynomial_space.evaluate(unit_point,
+                                    unit_point_values,
+                                    unit_point_grads,
+                                    unit_point_grad_grads,
+                                    empty_vector_of_3rd_order_tensors,
+                                    empty_vector_of_4th_order_tensors);
 
           // Declare the values, derivatives
           // and second derivatives vectors of
@@ -335,21 +335,21 @@ PolynomialsNedelec<dim>::evaluate(
           p1(0) = unit_point(1);
           p1(1) = unit_point(2);
           p1(2) = unit_point(0);
-          polynomial_space.compute(p1,
-                                   p1_values,
-                                   p1_grads,
-                                   p1_grad_grads,
-                                   empty_vector_of_3rd_order_tensors,
-                                   empty_vector_of_4th_order_tensors);
+          polynomial_space.evaluate(p1,
+                                    p1_values,
+                                    p1_grads,
+                                    p1_grad_grads,
+                                    empty_vector_of_3rd_order_tensors,
+                                    empty_vector_of_4th_order_tensors);
           p2(0) = unit_point(2);
           p2(1) = unit_point(0);
           p2(2) = unit_point(1);
-          polynomial_space.compute(p2,
-                                   p2_values,
-                                   p2_grads,
-                                   p2_grad_grads,
-                                   empty_vector_of_3rd_order_tensors,
-                                   empty_vector_of_4th_order_tensors);
+          polynomial_space.evaluate(p2,
+                                    p2_values,
+                                    p2_grads,
+                                    p2_grad_grads,
+                                    empty_vector_of_3rd_order_tensors,
+                                    empty_vector_of_4th_order_tensors);
 
           // Assign the correct values to the
           // corresponding shape functions.
diff --git a/source/base/polynomials_raviart_thomas.cc b/source/base/polynomials_raviart_thomas.cc
index 25a018c601..0affb59ff5 100644
--- a/source/base/polynomials_raviart_thomas.cc
+++ b/source/base/polynomials_raviart_thomas.cc
@@ -127,12 +127,12 @@ PolynomialsRaviartThomas<dim>::evaluate(
       for (unsigned int c = 0; c < dim; ++c)
         p(c) = unit_point((c + d) % dim);
 
-      polynomial_space.compute(p,
-                               p_values,
-                               p_grads,
-                               p_grad_grads,
-                               p_third_derivatives,
-                               p_fourth_derivatives);
+      polynomial_space.evaluate(p,
+                                p_values,
+                                p_grads,
+                                p_grad_grads,
+                                p_third_derivatives,
+                                p_fourth_derivatives);
 
       for (unsigned int i = 0; i < p_values.size(); ++i)
         values[i + d * n_sub][d] = p_values[i];
diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc
index 71d5f38186..49ce046c0d 100644
--- a/source/base/tensor_product_polynomials.cc
+++ b/source/base/tensor_product_polynomials.cc
@@ -15,6 +15,7 @@
 
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/polynomials_piecewise.h>
+#include <deal.II/base/std_cxx14/memory.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 
 #include <boost/container/small_vector.hpp>
@@ -98,7 +99,7 @@ TensorProductPolynomials<dim, PolynomialType>::output_indices(
   std::ostream &out) const
 {
   unsigned int ix[dim];
-  for (unsigned int i = 0; i < n_tensor_pols; ++i)
+  for (unsigned int i = 0; i < this->n(); ++i)
     {
       compute_index(i, ix);
       out << i << "\t";
@@ -250,26 +251,23 @@ TensorProductPolynomials<dim, PolynomialType>::evaluate(
   std::vector<Tensor<4, dim>> &fourth_derivatives) const
 {
   Assert(dim <= 3, ExcNotImplemented());
-  Assert(values.size() == n_tensor_pols || values.size() == 0,
-         ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
-  Assert(grads.size() == n_tensor_pols || grads.size() == 0,
-         ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
-  Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
-         ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
-  Assert(third_derivatives.size() == n_tensor_pols ||
-           third_derivatives.size() == 0,
-         ExcDimensionMismatch2(third_derivatives.size(), n_tensor_pols, 0));
-  Assert(fourth_derivatives.size() == n_tensor_pols ||
+  Assert(values.size() == this->n() || values.size() == 0,
+         ExcDimensionMismatch2(values.size(), this->n(), 0));
+  Assert(grads.size() == this->n() || grads.size() == 0,
+         ExcDimensionMismatch2(grads.size(), this->n(), 0));
+  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
+         ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
+  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
+         ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
+  Assert(fourth_derivatives.size() == this->n() ||
            fourth_derivatives.size() == 0,
-         ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));
+         ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
 
-  const bool update_values     = (values.size() == n_tensor_pols),
-             update_grads      = (grads.size() == n_tensor_pols),
-             update_grad_grads = (grad_grads.size() == n_tensor_pols),
-             update_3rd_derivatives =
-               (third_derivatives.size() == n_tensor_pols),
-             update_4th_derivatives =
-               (fourth_derivatives.size() == n_tensor_pols);
+  const bool update_values          = (values.size() == this->n()),
+             update_grads           = (grads.size() == this->n()),
+             update_grad_grads      = (grad_grads.size() == this->n()),
+             update_3rd_derivatives = (third_derivatives.size() == this->n()),
+             update_4th_derivatives = (fourth_derivatives.size() == this->n());
 
   // check how many values/derivatives we have to compute
   unsigned int n_values_and_derivatives = 0;
@@ -405,14 +403,24 @@ TensorProductPolynomials<dim, PolynomialType>::evaluate(
 
 
 
+template <int dim, typename PolynomialType>
+std::unique_ptr<ScalarPolynomialsBase<dim>>
+TensorProductPolynomials<dim, PolynomialType>::clone() const
+{
+  return std_cxx14::make_unique<TensorProductPolynomials<dim, PolynomialType>>(
+    *this);
+}
+
+
+
 /* ------------------- AnisotropicPolynomials -------------- */
 
 
 template <int dim>
 AnisotropicPolynomials<dim>::AnisotropicPolynomials(
   const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
-  : polynomials(pols)
-  , n_tensor_pols(get_n_tensor_pols(pols))
+  : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
+  , polynomials(pols)
 {
   Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
   for (unsigned int d = 0; d < dim; ++d)
@@ -531,7 +539,7 @@ AnisotropicPolynomials<dim>::compute_grad_grad(const unsigned int i,
 
 template <int dim>
 void
-AnisotropicPolynomials<dim>::compute(
+AnisotropicPolynomials<dim>::evaluate(
   const Point<dim> &           p,
   std::vector<double> &        values,
   std::vector<Tensor<1, dim>> &grads,
@@ -539,26 +547,23 @@ AnisotropicPolynomials<dim>::compute(
   std::vector<Tensor<3, dim>> &third_derivatives,
   std::vector<Tensor<4, dim>> &fourth_derivatives) const
 {
-  Assert(values.size() == n_tensor_pols || values.size() == 0,
-         ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
-  Assert(grads.size() == n_tensor_pols || grads.size() == 0,
-         ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
-  Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
-         ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
-  Assert(third_derivatives.size() == n_tensor_pols ||
-           third_derivatives.size() == 0,
-         ExcDimensionMismatch2(third_derivatives.size(), n_tensor_pols, 0));
-  Assert(fourth_derivatives.size() == n_tensor_pols ||
+  Assert(values.size() == this->n() || values.size() == 0,
+         ExcDimensionMismatch2(values.size(), this->n(), 0));
+  Assert(grads.size() == this->n() || grads.size() == 0,
+         ExcDimensionMismatch2(grads.size(), this->n(), 0));
+  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
+         ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
+  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
+         ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
+  Assert(fourth_derivatives.size() == this->n() ||
            fourth_derivatives.size() == 0,
-         ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));
+         ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
 
-  const bool update_values     = (values.size() == n_tensor_pols),
-             update_grads      = (grads.size() == n_tensor_pols),
-             update_grad_grads = (grad_grads.size() == n_tensor_pols),
-             update_3rd_derivatives =
-               (third_derivatives.size() == n_tensor_pols),
-             update_4th_derivatives =
-               (fourth_derivatives.size() == n_tensor_pols);
+  const bool update_values          = (values.size() == this->n()),
+             update_grads           = (grads.size() == this->n()),
+             update_grad_grads      = (grad_grads.size() == this->n()),
+             update_3rd_derivatives = (third_derivatives.size() == this->n()),
+             update_4th_derivatives = (fourth_derivatives.size() == this->n());
 
   // check how many
   // values/derivatives we have to
@@ -590,7 +595,7 @@ AnisotropicPolynomials<dim>::compute(
         }
     }
 
-  for (unsigned int i = 0; i < n_tensor_pols; ++i)
+  for (unsigned int i = 0; i < this->n(); ++i)
     {
       // first get the
       // one-dimensional indices of
@@ -679,15 +684,6 @@ AnisotropicPolynomials<dim>::compute(
 }
 
 
-
-template <int dim>
-unsigned int
-AnisotropicPolynomials<dim>::n() const
-{
-  return n_tensor_pols;
-}
-
-
 template <int dim>
 unsigned int
 AnisotropicPolynomials<dim>::get_n_tensor_pols(
@@ -700,6 +696,14 @@ AnisotropicPolynomials<dim>::get_n_tensor_pols(
 }
 
 
+template <int dim>
+std::unique_ptr<ScalarPolynomialsBase<dim>>
+AnisotropicPolynomials<dim>::clone() const
+{
+  return std_cxx14::make_unique<AnisotropicPolynomials<dim>>(*this);
+}
+
+
 
 /* ------------------- explicit instantiations -------------- */
 template class TensorProductPolynomials<1, Polynomials::Polynomial<double>>;
diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc
index dd231f8304..bd70e3e2d2 100644
--- a/source/base/tensor_product_polynomials_bubbles.cc
+++ b/source/base/tensor_product_polynomials_bubbles.cc
@@ -15,6 +15,7 @@
 
 
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/std_cxx14/memory.h>
 #include <deal.II/base/tensor_product_polynomials_bubbles.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -70,7 +71,6 @@ TensorProductPolynomialsBubbles<dim>::compute_value(const unsigned int i,
   const unsigned int q_degree      = tensor_polys.polynomials.size() - 1;
   const unsigned int max_q_indices = tensor_polys.n();
   const unsigned int n_bubbles     = ((q_degree <= 1) ? 1 : dim);
-  (void)n_bubbles;
   Assert(i < max_q_indices + n_bubbles, ExcInternalError());
 
   // treat the regular basis functions
@@ -110,7 +110,6 @@ TensorProductPolynomialsBubbles<dim>::compute_grad(const unsigned int i,
   const unsigned int q_degree      = tensor_polys.polynomials.size() - 1;
   const unsigned int max_q_indices = tensor_polys.n();
   const unsigned int n_bubbles     = ((q_degree <= 1) ? 1 : dim);
-  (void)n_bubbles;
   Assert(i < max_q_indices + n_bubbles, ExcInternalError());
 
   // treat the regular basis functions
@@ -158,7 +157,6 @@ TensorProductPolynomialsBubbles<dim>::compute_grad_grad(
   const unsigned int q_degree      = tensor_polys.polynomials.size() - 1;
   const unsigned int max_q_indices = tensor_polys.n();
   const unsigned int n_bubbles     = ((q_degree <= 1) ? 1 : dim);
-  (void)n_bubbles;
   Assert(i < max_q_indices + n_bubbles, ExcInternalError());
 
   // treat the regular basis functions
@@ -266,8 +264,7 @@ TensorProductPolynomialsBubbles<dim>::evaluate(
 {
   const unsigned int q_degree      = tensor_polys.polynomials.size() - 1;
   const unsigned int max_q_indices = tensor_polys.n();
-  (void)max_q_indices;
-  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
+  const unsigned int n_bubbles     = ((q_degree <= 1) ? 1 : dim);
   Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
          ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
   Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
@@ -333,6 +330,15 @@ TensorProductPolynomialsBubbles<dim>::evaluate(
 }
 
 
+
+template <int dim>
+std::unique_ptr<ScalarPolynomialsBase<dim>>
+TensorProductPolynomialsBubbles<dim>::clone() const
+{
+  return std_cxx14::make_unique<TensorProductPolynomialsBubbles<dim>>(*this);
+}
+
+
 /* ------------------- explicit instantiations -------------- */
 template class TensorProductPolynomialsBubbles<1>;
 template class TensorProductPolynomialsBubbles<2>;
diff --git a/source/base/tensor_product_polynomials_const.cc b/source/base/tensor_product_polynomials_const.cc
index 63bce45b98..54286ba84c 100644
--- a/source/base/tensor_product_polynomials_const.cc
+++ b/source/base/tensor_product_polynomials_const.cc
@@ -29,7 +29,7 @@ double
 TensorProductPolynomialsConst<dim>::compute_value(const unsigned int i,
                                                   const Point<dim> & p) const
 {
-  const unsigned int max_indices = this->n_tensor_pols;
+  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
@@ -57,7 +57,7 @@ Tensor<1, dim>
 TensorProductPolynomialsConst<dim>::compute_grad(const unsigned int i,
                                                  const Point<dim> & p) const
 {
-  const unsigned int max_indices = this->n_tensor_pols;
+  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
@@ -73,7 +73,7 @@ Tensor<2, dim>
 TensorProductPolynomialsConst<dim>::compute_grad_grad(const unsigned int i,
                                                       const Point<dim> &p) const
 {
-  const unsigned int max_indices = this->n_tensor_pols;
+  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
@@ -94,21 +94,30 @@ TensorProductPolynomialsConst<dim>::evaluate(
   std::vector<Tensor<3, dim>> &third_derivatives,
   std::vector<Tensor<4, dim>> &fourth_derivatives) const
 {
-  Assert(values.size() == this->n_tensor_pols + 1 || values.size() == 0,
-         ExcDimensionMismatch2(values.size(), this->n_tensor_pols + 1, 0));
-  Assert(grads.size() == this->n_tensor_pols + 1 || grads.size() == 0,
-         ExcDimensionMismatch2(grads.size(), this->n_tensor_pols + 1, 0));
-  Assert(grad_grads.size() == this->n_tensor_pols + 1 || grad_grads.size() == 0,
-         ExcDimensionMismatch2(grad_grads.size(), this->n_tensor_pols + 1, 0));
-  Assert(third_derivatives.size() == this->n_tensor_pols + 1 ||
+  Assert(values.size() == TensorProductPolynomials<dim>::n() + 1 ||
+           values.size() == 0,
+         ExcDimensionMismatch2(values.size(),
+                               TensorProductPolynomials<dim>::n() + 1,
+                               0));
+  Assert(grads.size() == TensorProductPolynomials<dim>::n() + 1 ||
+           grads.size() == 0,
+         ExcDimensionMismatch2(grads.size(),
+                               TensorProductPolynomials<dim>::n() + 1,
+                               0));
+  Assert(grad_grads.size() == TensorProductPolynomials<dim>::n() + 1 ||
+           grad_grads.size() == 0,
+         ExcDimensionMismatch2(grad_grads.size(),
+                               TensorProductPolynomials<dim>::n() + 1,
+                               0));
+  Assert(third_derivatives.size() == TensorProductPolynomials<dim>::n() + 1 ||
            third_derivatives.size() == 0,
          ExcDimensionMismatch2(third_derivatives.size(),
-                               this->n_tensor_pols + 1,
+                               TensorProductPolynomials<dim>::n() + 1,
                                0));
-  Assert(fourth_derivatives.size() == this->n_tensor_pols + 1 ||
+  Assert(fourth_derivatives.size() == TensorProductPolynomials<dim>::n() + 1 ||
            fourth_derivatives.size() == 0,
          ExcDimensionMismatch2(fourth_derivatives.size(),
-                               this->n_tensor_pols + 1,
+                               TensorProductPolynomials<dim>::n() + 1,
                                0));
 
   // remove slot for const value, go into the base class compute method and
@@ -132,12 +141,12 @@ TensorProductPolynomialsConst<dim>::evaluate(
     }
   if (third_derivatives.empty() == false)
     {
-      third_derivatives.resize(this->n_tensor_pols);
+      third_derivatives.resize(TensorProductPolynomials<dim>::n());
       do_3rd_derivatives = true;
     }
   if (fourth_derivatives.empty() == false)
     {
-      fourth_derivatives.resize(this->n_tensor_pols);
+      fourth_derivatives.resize(TensorProductPolynomials<dim>::n());
       do_4th_derivatives = true;
     }
 
diff --git a/tests/base/anisotropic_1.cc b/tests/base/anisotropic_1.cc
index e6f3327a91..63cd366303 100644
--- a/tests/base/anisotropic_1.cc
+++ b/tests/base/anisotropic_1.cc
@@ -43,7 +43,7 @@ check_poly(const Point<dim> &     x,
   std::vector<Tensor<4, dim>> fourth1(n), fourth2(n);
 
   p.evaluate(x, values1, gradients1, second1, third1, fourth1);
-  q.compute(x, values2, gradients2, second2, third2, fourth2);
+  q.evaluate(x, values2, gradients2, second2, third2, fourth2);
 
   for (unsigned int k = 0; k < n; ++k)
     {