From 3b3e9b2b4d1f609718b99d449ce8793df087ed97 Mon Sep 17 00:00:00 2001
From: grahambenharper <grahambenharper@gmail.com>
Date: Sat, 7 Sep 2019 17:49:59 -0600
Subject: [PATCH] Remove PolynomialType template argument from FE_Poly

---
 include/deal.II/base/polynomial_space.h       |   6 +-
 include/deal.II/base/polynomials_adini.h      |   6 +-
 .../base/polynomials_rannacher_turek.h        |   6 +-
 .../deal.II/base/scalar_polynomials_base.h    |  58 +++++
 .../deal.II/base/tensor_product_polynomials.h |  12 +-
 .../base/tensor_product_polynomials_bubbles.h |   6 +-
 .../base/tensor_product_polynomials_const.h   |   6 +-
 include/deal.II/fe/fe_dgp.h                   |   2 +-
 include/deal.II/fe/fe_dgp_monomial.h          |   2 +-
 include/deal.II/fe/fe_dgq.h                   |   2 +-
 include/deal.II/fe/fe_poly.h                  |  29 +--
 include/deal.II/fe/fe_poly.templates.h        | 221 ++++++++----------
 include/deal.II/fe/fe_q_base.h                |   2 +-
 include/deal.II/fe/fe_q_hierarchical.h        |   2 +-
 include/deal.II/fe/fe_rannacher_turek.h       |   2 +-
 .../matrix_free/matrix_free.templates.h       |   8 +-
 .../matrix_free/shape_info.templates.h        |  19 +-
 source/fe/fe_dgp.cc                           |   2 +-
 source/fe/fe_dgp_monomial.cc                  |  22 +-
 source/fe/fe_dgq.cc                           |  23 +-
 source/fe/fe_poly.cc                          | 137 +----------
 source/fe/fe_poly.inst.in                     |  24 +-
 source/fe/fe_q.cc                             |   4 +-
 source/fe/fe_q_base.cc                        |  74 ++++--
 source/fe/fe_q_hierarchical.cc                |  18 +-
 source/fe/fe_rannacher_turek.cc               |  15 +-
 26 files changed, 316 insertions(+), 392 deletions(-)

diff --git a/include/deal.II/base/polynomial_space.h b/include/deal.II/base/polynomial_space.h
index ed1a0a5cfd..4e050b6225 100644
--- a/include/deal.II/base/polynomial_space.h
+++ b/include/deal.II/base/polynomial_space.h
@@ -157,7 +157,7 @@ public:
    * Consider using evaluate() instead.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the <tt>order</tt>th derivative of the <tt>i</tt>th polynomial
@@ -178,7 +178,7 @@ public:
    * Consider using evaluate() instead.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the second derivative (grad_grad) of the <tt>i</tt>th polynomial
@@ -187,7 +187,7 @@ public:
    * Consider using evaluate() instead.
    */
   Tensor<2, dim>
-  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Return the number of polynomials spanning the space represented by this
diff --git a/include/deal.II/base/polynomials_adini.h b/include/deal.II/base/polynomials_adini.h
index 4b9c0b6a06..64daffad89 100644
--- a/include/deal.II/base/polynomials_adini.h
+++ b/include/deal.II/base/polynomials_adini.h
@@ -78,7 +78,7 @@ public:
    * Consider using evaluate() instead.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the gradient of the <tt>i</tt>th polynomial at
@@ -87,7 +87,7 @@ public:
    * Consider using evaluate() instead.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the second derivative (grad_grad) of the <tt>i</tt>th polynomial
@@ -96,7 +96,7 @@ public:
    * Consider using evaluate() instead.
    */
   Tensor<2, dim>
-  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Return the name of the space, which is <tt>PolynomialsAdini</tt>.
diff --git a/include/deal.II/base/polynomials_rannacher_turek.h b/include/deal.II/base/polynomials_rannacher_turek.h
index 0b56598ba1..4fb378458b 100644
--- a/include/deal.II/base/polynomials_rannacher_turek.h
+++ b/include/deal.II/base/polynomials_rannacher_turek.h
@@ -58,7 +58,7 @@ public:
    * Value of basis function @p i at @p p.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * <tt>order</tt>-th of basis function @p i at @p p.
@@ -73,13 +73,13 @@ public:
    * Gradient of basis function @p i at @p p.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Gradient of gradient of basis function @p i at @p p.
    */
   Tensor<2, dim>
-  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute values and derivatives of all basis functions at @p unit_point.
diff --git a/include/deal.II/base/scalar_polynomials_base.h b/include/deal.II/base/scalar_polynomials_base.h
index 5b192a885c..dc702bbc06 100644
--- a/include/deal.II/base/scalar_polynomials_base.h
+++ b/include/deal.II/base/scalar_polynomials_base.h
@@ -106,6 +106,64 @@ public:
            std::vector<Tensor<3, dim>> &third_derivatives,
            std::vector<Tensor<4, dim>> &fourth_derivatives) const = 0;
 
+  /**
+   * Compute the value of the <tt>i</tt>th polynomial at unit point
+   * <tt>p</tt>.
+   *
+   * Consider using compute() instead.
+   */
+  virtual double
+  compute_value(const unsigned int /*i*/, const Point<dim> & /*p*/) const
+  {
+    Assert(false, ExcNotImplemented());
+    return 0;
+  }
+
+  /**
+   * Compute the <tt>order</tt>th derivative of the <tt>i</tt>th polynomial
+   * at unit point <tt>p</tt>.
+   *
+   * Consider using compute() instead.
+   *
+   * @tparam order The order of the derivative.
+   */
+  template <int order>
+  Tensor<order, dim>
+  compute_derivative(const unsigned int /*i*/, const Point<dim> & /*p*/) const
+  {
+    Assert(false, ExcNotImplemented());
+    Tensor<order, dim> empty;
+    return empty;
+  }
+
+  /**
+   * Compute the gradient of the <tt>i</tt>th polynomial at unit point
+   * <tt>p</tt>.
+   *
+   * Consider using compute() instead.
+   */
+  virtual Tensor<1, dim>
+  compute_grad(const unsigned int /*i*/, const Point<dim> & /*p*/) const
+  {
+    Assert(false, ExcNotImplemented());
+    Tensor<1, dim> empty;
+    return empty;
+  }
+
+  /**
+   * Compute the second derivative (grad_grad) of the <tt>i</tt>th polynomial
+   * at unit point <tt>p</tt>.
+   *
+   * Consider using compute() instead.
+   */
+  virtual Tensor<2, dim>
+  compute_grad_grad(const unsigned int /*i*/, const Point<dim> & /*p*/) const
+  {
+    Assert(false, ExcNotImplemented());
+    Tensor<2, dim> empty;
+    return empty;
+  }
+
   /**
    * Return the number of polynomials.
    */
diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h
index c6bec0062d..e024553d69 100644
--- a/include/deal.II/base/tensor_product_polynomials.h
+++ b/include/deal.II/base/tensor_product_polynomials.h
@@ -147,7 +147,7 @@ public:
    * polynomials all at once and in a much more efficient way.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the <tt>order</tt>th derivative of the <tt>i</tt>th tensor
@@ -180,7 +180,7 @@ public:
    * polynomials all at once and in a much more efficient way.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the second derivative (grad_grad) of the <tt>i</tt>th tensor
@@ -195,7 +195,7 @@ public:
    * 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;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Return the name of the space, which is <tt>TensorProductPolynomials</tt>.
@@ -340,7 +340,7 @@ public:
    * tensor polynomials all at once and in a much more efficient way.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the <tt>order</tt>th derivative of the <tt>i</tt>th tensor
@@ -373,7 +373,7 @@ public:
    * tensor polynomials all at once and in a much more efficient way.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the second derivative (grad_grad) of the <tt>i</tt>th tensor
@@ -388,7 +388,7 @@ public:
    * 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;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Return the name of the space, which is <tt>AnisotropicPolynomials</tt>.
diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h
index 1bfcaaf4d4..0e8d28823f 100644
--- a/include/deal.II/base/tensor_product_polynomials_bubbles.h
+++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h
@@ -120,7 +120,7 @@ public:
    * polynomials all at once and in a much more efficient way.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the order @p order derivative of the <tt>i</tt>th tensor product
@@ -151,7 +151,7 @@ public:
    * polynomials all at once and in a much more efficient way.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the second derivative (grad_grad) of the <tt>i</tt>th tensor
@@ -166,7 +166,7 @@ public:
    * 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;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Return the number of tensor product polynomials plus the bubble
diff --git a/include/deal.II/base/tensor_product_polynomials_const.h b/include/deal.II/base/tensor_product_polynomials_const.h
index dc069adef4..0c75e2f59c 100644
--- a/include/deal.II/base/tensor_product_polynomials_const.h
+++ b/include/deal.II/base/tensor_product_polynomials_const.h
@@ -121,7 +121,7 @@ public:
    * polynomials all at once and in a much more efficient way.
    */
   double
-  compute_value(const unsigned int i, const Point<dim> &p) const;
+  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the <tt>order</tt>th derivative of the <tt>i</tt>th tensor
@@ -154,7 +154,7 @@ public:
    * polynomials all at once and in a much more efficient way.
    */
   Tensor<1, dim>
-  compute_grad(const unsigned int i, const Point<dim> &p) const;
+  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Compute the second derivative (grad_grad) of the <tt>i</tt>th tensor
@@ -169,7 +169,7 @@ public:
    * 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;
+  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
   /**
    * Return the number of tensor product polynomials plus the constant
diff --git a/include/deal.II/fe/fe_dgp.h b/include/deal.II/fe/fe_dgp.h
index 355252111a..dcfbac266c 100644
--- a/include/deal.II/fe/fe_dgp.h
+++ b/include/deal.II/fe/fe_dgp.h
@@ -308,7 +308,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Guido Kanschat, 2001, 2002, Ralf Hartmann 2004
  */
 template <int dim, int spacedim = dim>
-class FE_DGP : public FE_Poly<PolynomialSpace<dim>, dim, spacedim>
+class FE_DGP : public FE_Poly<dim, spacedim>
 {
 public:
   /**
diff --git a/include/deal.II/fe/fe_dgp_monomial.h b/include/deal.II/fe/fe_dgp_monomial.h
index 16dffa8ed8..fca49e7da9 100644
--- a/include/deal.II/fe/fe_dgp_monomial.h
+++ b/include/deal.II/fe/fe_dgp_monomial.h
@@ -283,7 +283,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Ralf Hartmann, 2004
  */
 template <int dim>
-class FE_DGPMonomial : public FE_Poly<PolynomialsP<dim>, dim>
+class FE_DGPMonomial : public FE_Poly<dim>
 {
 public:
   /**
diff --git a/include/deal.II/fe/fe_dgq.h b/include/deal.II/fe/fe_dgq.h
index 26705d2c3c..2c9a697e44 100644
--- a/include/deal.II/fe/fe_dgq.h
+++ b/include/deal.II/fe/fe_dgq.h
@@ -109,7 +109,7 @@ class Quadrature;
  * @author Ralf Hartmann, Guido Kanschat 2001, 2004
  */
 template <int dim, int spacedim = dim>
-class FE_DGQ : public FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>
+class FE_DGQ : public FE_Poly<dim, spacedim>
 {
 public:
   /**
diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h
index 77efe6a4a5..cfa174d72d 100644
--- a/include/deal.II/fe/fe_poly.h
+++ b/include/deal.II/fe/fe_poly.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/quadrature.h>
+#include <deal.II/base/scalar_polynomials_base.h>
 #include <deal.II/base/std_cxx14/memory.h>
 
 #include <deal.II/fe/fe.h>
@@ -70,20 +71,23 @@ DEAL_II_NAMESPACE_OPEN
  * @author Ralf Hartmann 2004, Guido Kanschat, 2009
  */
 
-template <class PolynomialType,
-          int dim      = PolynomialType::dimension,
-          int spacedim = dim>
+template <int dim, int spacedim = dim>
 class FE_Poly : public FiniteElement<dim, spacedim>
 {
 public:
   /**
    * Constructor.
    */
-  FE_Poly(const PolynomialType &            poly_space,
+  FE_Poly(const ScalarPolynomialsBase<dim> &poly_space,
           const FiniteElementData<dim> &    fe_data,
           const std::vector<bool> &         restriction_is_additive_flags,
           const std::vector<ComponentMask> &nonzero_components);
 
+  /**
+   * Copy constructor.
+   */
+  FE_Poly(const FE_Poly &fe);
+
   /**
    * Return the polynomial degree of this finite element, i.e. the value
    * passed to the constructor.
@@ -311,12 +315,12 @@ protected:
                         update_3rd_derivatives))
       for (unsigned int i = 0; i < n_q_points; ++i)
         {
-          poly_space.evaluate(quadrature.point(i),
-                              values,
-                              grads,
-                              grad_grads,
-                              third_derivatives,
-                              fourth_derivatives);
+          poly_space->evaluate(quadrature.point(i),
+                               values,
+                               grads,
+                               grad_grads,
+                               third_derivatives,
+                               fourth_derivatives);
 
           // the values of shape functions at quadrature points don't change.
           // consequently, write these values right into the output array if
@@ -510,10 +514,9 @@ protected:
 
 
   /**
-   * The polynomial space. Its type is given by the template parameter
-   * PolynomialType.
+   * The polynomial space.
    */
-  PolynomialType poly_space;
+  const std::unique_ptr<ScalarPolynomialsBase<dim>> poly_space;
 };
 
 /*@}*/
diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h
index 3cea1551b1..4285aeae0a 100644
--- a/include/deal.II/fe/fe_poly.templates.h
+++ b/include/deal.II/fe/fe_poly.templates.h
@@ -32,42 +32,40 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-template <class PolynomialType, int dim, int spacedim>
-FE_Poly<PolynomialType, dim, spacedim>::FE_Poly(
-  const PolynomialType &            poly_space,
+template <int dim, int spacedim>
+FE_Poly<dim, spacedim>::FE_Poly(
+  const ScalarPolynomialsBase<dim> &poly_space,
   const FiniteElementData<dim> &    fe_data,
   const std::vector<bool> &         restriction_is_additive_flags,
   const std::vector<ComponentMask> &nonzero_components)
   : FiniteElement<dim, spacedim>(fe_data,
                                  restriction_is_additive_flags,
                                  nonzero_components)
-  , poly_space(poly_space)
-{
-  AssertDimension(dim, PolynomialType::dimension);
-}
+  , poly_space(poly_space.clone())
+{}
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 unsigned int
-FE_Poly<PolynomialType, dim, spacedim>::get_degree() const
+FE_Poly<dim, spacedim>::get_degree() const
 {
   return this->degree;
 }
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 double
-FE_Poly<PolynomialType, dim, spacedim>::shape_value(const unsigned int i,
-                                                    const Point<dim> & p) const
+FE_Poly<dim, spacedim>::shape_value(const unsigned int i,
+                                    const Point<dim> & p) const
 {
   AssertIndexRange(i, this->dofs_per_cell);
-  return poly_space.compute_value(i, p);
+  return poly_space->compute_value(i, p);
 }
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 double
-FE_Poly<PolynomialType, dim, spacedim>::shape_value_component(
+FE_Poly<dim, spacedim>::shape_value_component(
   const unsigned int i,
   const Point<dim> & p,
   const unsigned int component) const
@@ -75,52 +73,50 @@ FE_Poly<PolynomialType, dim, spacedim>::shape_value_component(
   (void)component;
   AssertIndexRange(i, this->dofs_per_cell);
   AssertIndexRange(component, 1);
-  return poly_space.compute_value(i, p);
+  return poly_space->compute_value(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<1, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_grad(const unsigned int i,
-                                                   const Point<dim> & p) const
+FE_Poly<dim, spacedim>::shape_grad(const unsigned int i,
+                                   const Point<dim> & p) const
 {
   AssertIndexRange(i, this->dofs_per_cell);
-  return poly_space.template compute_derivative<1>(i, p);
+  return poly_space->template compute_derivative<1>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<1, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_grad_component(
-  const unsigned int i,
-  const Point<dim> & p,
-  const unsigned int component) const
+FE_Poly<dim, spacedim>::shape_grad_component(const unsigned int i,
+                                             const Point<dim> & p,
+                                             const unsigned int component) const
 {
   (void)component;
   AssertIndexRange(i, this->dofs_per_cell);
   AssertIndexRange(component, 1);
-  return poly_space.template compute_derivative<1>(i, p);
+  return poly_space->template compute_derivative<1>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<2, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_grad_grad(
-  const unsigned int i,
-  const Point<dim> & p) const
+FE_Poly<dim, spacedim>::shape_grad_grad(const unsigned int i,
+                                        const Point<dim> & p) const
 {
   AssertIndexRange(i, this->dofs_per_cell);
-  return poly_space.template compute_derivative<2>(i, p);
+  return poly_space->template compute_derivative<2>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<2, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_grad_grad_component(
+FE_Poly<dim, spacedim>::shape_grad_grad_component(
   const unsigned int i,
   const Point<dim> & p,
   const unsigned int component) const
@@ -128,26 +124,25 @@ FE_Poly<PolynomialType, dim, spacedim>::shape_grad_grad_component(
   (void)component;
   AssertIndexRange(i, this->dofs_per_cell);
   AssertIndexRange(component, 1);
-  return poly_space.template compute_derivative<2>(i, p);
+  return poly_space->template compute_derivative<2>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<3, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_3rd_derivative(
-  const unsigned int i,
-  const Point<dim> & p) const
+FE_Poly<dim, spacedim>::shape_3rd_derivative(const unsigned int i,
+                                             const Point<dim> & p) const
 {
   AssertIndexRange(i, this->dofs_per_cell);
-  return poly_space.template compute_derivative<3>(i, p);
+  return poly_space->template compute_derivative<3>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<3, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_3rd_derivative_component(
+FE_Poly<dim, spacedim>::shape_3rd_derivative_component(
   const unsigned int i,
   const Point<dim> & p,
   const unsigned int component) const
@@ -155,26 +150,25 @@ FE_Poly<PolynomialType, dim, spacedim>::shape_3rd_derivative_component(
   (void)component;
   AssertIndexRange(i, this->dofs_per_cell);
   AssertIndexRange(component, 1);
-  return poly_space.template compute_derivative<3>(i, p);
+  return poly_space->template compute_derivative<3>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<4, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_4th_derivative(
-  const unsigned int i,
-  const Point<dim> & p) const
+FE_Poly<dim, spacedim>::shape_4th_derivative(const unsigned int i,
+                                             const Point<dim> & p) const
 {
   AssertIndexRange(i, this->dofs_per_cell);
-  return poly_space.template compute_derivative<4>(i, p);
+  return poly_space->template compute_derivative<4>(i, p);
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<4, dim>
-FE_Poly<PolynomialType, dim, spacedim>::shape_4th_derivative_component(
+FE_Poly<dim, spacedim>::shape_4th_derivative_component(
   const unsigned int i,
   const Point<dim> & p,
   const unsigned int component) const
@@ -182,7 +176,7 @@ FE_Poly<PolynomialType, dim, spacedim>::shape_4th_derivative_component(
   (void)component;
   AssertIndexRange(i, this->dofs_per_cell);
   AssertIndexRange(component, 1);
-  return poly_space.template compute_derivative<4>(i, p);
+  return poly_space->template compute_derivative<4>(i, p);
 }
 
 
@@ -192,10 +186,9 @@ FE_Poly<PolynomialType, dim, spacedim>::shape_4th_derivative_component(
 //---------------------------------------------------------------------------
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 UpdateFlags
-FE_Poly<PolynomialType, dim, spacedim>::requires_update_flags(
-  const UpdateFlags flags) const
+FE_Poly<dim, spacedim>::requires_update_flags(const UpdateFlags flags) const
 {
   UpdateFlags out = update_default;
 
@@ -264,9 +257,9 @@ higher_derivatives_need_correcting(
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 void
-FE_Poly<PolynomialType, dim, spacedim>::fill_fe_values(
+FE_Poly<dim, spacedim>::fill_fe_values(
   const typename Triangulation<dim, spacedim>::cell_iterator &,
   const CellSimilarity::Similarity                         cell_similarity,
   const Quadrature<dim> &                                  quadrature,
@@ -336,9 +329,9 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_values(
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 void
-FE_Poly<PolynomialType, dim, spacedim>::fill_fe_face_values(
+FE_Poly<dim, spacedim>::fill_fe_face_values(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
   const unsigned int                                          face_no,
   const Quadrature<dim - 1> &                                 quadrature,
@@ -427,9 +420,9 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_face_values(
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 void
-FE_Poly<PolynomialType, dim, spacedim>::fill_fe_subface_values(
+FE_Poly<dim, spacedim>::fill_fe_subface_values(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
   const unsigned int                                          face_no,
   const unsigned int                                          sub_no,
@@ -521,9 +514,9 @@ FE_Poly<PolynomialType, dim, spacedim>::fill_fe_subface_values(
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 inline void
-FE_Poly<PolynomialType, dim, spacedim>::correct_hessians(
+FE_Poly<dim, spacedim>::correct_hessians(
   internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
     &output_data,
   const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
@@ -540,9 +533,9 @@ FE_Poly<PolynomialType, dim, spacedim>::correct_hessians(
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 inline void
-FE_Poly<PolynomialType, dim, spacedim>::correct_third_derivatives(
+FE_Poly<dim, spacedim>::correct_third_derivatives(
   internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
     &output_data,
   const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
@@ -569,82 +562,62 @@ FE_Poly<PolynomialType, dim, spacedim>::correct_third_derivatives(
               }
 }
 
-namespace internal
-{
-  template <class PolynomialType>
-  inline std::vector<unsigned int>
-  get_poly_space_numbering(const PolynomialType &)
-  {
-    Assert(false, ExcNotImplemented());
-    return std::vector<unsigned int>();
-  }
-
-  template <class PolynomialType>
-  inline std::vector<unsigned int>
-  get_poly_space_numbering_inverse(const PolynomialType &)
-  {
-    Assert(false, ExcNotImplemented());
-    return std::vector<unsigned int>();
-  }
-
-  template <int dim, typename PolynomialType>
-  inline std::vector<unsigned int>
-  get_poly_space_numbering(
-    const TensorProductPolynomials<dim, PolynomialType> &poly)
-  {
-    return poly.get_numbering();
-  }
-
-  template <int dim, typename PolynomialType>
-  inline std::vector<unsigned int>
-  get_poly_space_numbering_inverse(
-    const TensorProductPolynomials<dim, PolynomialType> &poly)
-  {
-    return poly.get_numbering_inverse();
-  }
-
-  template <int dim>
-  inline std::vector<unsigned int>
-  get_poly_space_numbering(const TensorProductPolynomialsConst<dim> &poly)
-  {
-    return poly.get_numbering();
-  }
-
-  template <int dim>
-  inline std::vector<unsigned int>
-  get_poly_space_numbering_inverse(
-    const TensorProductPolynomialsConst<dim> &poly)
-  {
-    return poly.get_numbering_inverse();
-  }
-} // namespace internal
-
-
-
-template <class PolynomialType, int dim, int spacedim>
+
+
+template <int dim, int spacedim>
 std::vector<unsigned int>
-FE_Poly<PolynomialType, dim, spacedim>::get_poly_space_numbering() const
+FE_Poly<dim, spacedim>::get_poly_space_numbering() const
 {
-  return internal::get_poly_space_numbering(poly_space);
+  auto *const space_tensor_prod =
+    dynamic_cast<TensorProductPolynomials<dim, ScalarPolynomialsBase<dim>> *>(
+      this->poly_space.get());
+
+  if (space_tensor_prod != nullptr)
+    return space_tensor_prod->get_numbering();
+
+  auto *const space_tensor_prod_const =
+    dynamic_cast<TensorProductPolynomials<dim, ScalarPolynomialsBase<dim>> *>(
+      this->poly_space.get());
+
+  if (space_tensor_prod_const != nullptr)
+    return space_tensor_prod_const->get_numbering();
+
+  Assert(false, ExcNotImplemented());
+  return std::vector<unsigned int>();
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 std::vector<unsigned int>
-FE_Poly<PolynomialType, dim, spacedim>::get_poly_space_numbering_inverse() const
+FE_Poly<dim, spacedim>::get_poly_space_numbering_inverse() const
 {
-  return internal::get_poly_space_numbering_inverse(poly_space);
+  auto *const space_tensor_prod =
+    dynamic_cast<TensorProductPolynomials<dim, ScalarPolynomialsBase<dim>> *>(
+      this->poly_space.get());
+
+  if (space_tensor_prod != nullptr)
+    return space_tensor_prod->get_numbering_inverse();
+
+  auto *const space_tensor_prod_const =
+    dynamic_cast<TensorProductPolynomials<dim, ScalarPolynomialsBase<dim>> *>(
+      this->poly_space.get());
+
+  if (space_tensor_prod_const != nullptr)
+    return space_tensor_prod_const->get_numbering_inverse();
+
+  Assert(false, ExcNotImplemented());
+  return std::vector<unsigned int>();
 }
 
 
 
-template <class PolynomialType, int dim, int spacedim>
+template <int dim, int spacedim>
 std::size_t
-FE_Poly<PolynomialType, dim, spacedim>::memory_consumption() const
+FE_Poly<dim, spacedim>::memory_consumption() const
 {
   return FiniteElement<dim, spacedim>::memory_consumption() +
-         poly_space.memory_consumption();
+         poly_space->memory_consumption();
 }
 
 
diff --git a/include/deal.II/fe/fe_q_base.h b/include/deal.II/fe/fe_q_base.h
index 5828c6e970..998c928f56 100644
--- a/include/deal.II/fe/fe_q_base.h
+++ b/include/deal.II/fe/fe_q_base.h
@@ -41,7 +41,7 @@ DEAL_II_NAMESPACE_OPEN
 template <class PolynomialType,
           int dim      = PolynomialType::dimension,
           int spacedim = dim>
-class FE_Q_Base : public FE_Poly<PolynomialType, dim, spacedim>
+class FE_Q_Base : public FE_Poly<dim, spacedim>
 {
 public:
   /**
diff --git a/include/deal.II/fe/fe_q_hierarchical.h b/include/deal.II/fe/fe_q_hierarchical.h
index 250e98bb59..738ce3aa5d 100644
--- a/include/deal.II/fe/fe_q_hierarchical.h
+++ b/include/deal.II/fe/fe_q_hierarchical.h
@@ -540,7 +540,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Brian Carnes, 2002, Ralf Hartmann 2004, 2005, Denis Davydov, 2015
  */
 template <int dim>
-class FE_Q_Hierarchical : public FE_Poly<TensorProductPolynomials<dim>, dim>
+class FE_Q_Hierarchical : public FE_Poly<dim>
 {
 public:
   /**
diff --git a/include/deal.II/fe/fe_rannacher_turek.h b/include/deal.II/fe/fe_rannacher_turek.h
index fedaa1e8a0..0ca8a0204f 100644
--- a/include/deal.II/fe/fe_rannacher_turek.h
+++ b/include/deal.II/fe/fe_rannacher_turek.h
@@ -81,7 +81,7 @@ DEAL_II_NAMESPACE_OPEN
  * @date 2015
  */
 template <int dim>
-class FE_RannacherTurek : public FE_Poly<PolynomialsRannacherTurek<dim>, dim>
+class FE_RannacherTurek : public FE_Poly<dim>
 {
 public:
   /**
diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index ca57050fc0..5827480fba 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -555,13 +555,9 @@ MatrixFree<dim, Number, VectorizedArrayType>::is_supported(
     return false;
 
   // then check of the base element is supported
-  if (dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>
-                     *>(fe_ptr) != nullptr)
+  if (dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr) != nullptr)
     return true;
-  if (dynamic_cast<const FE_Poly<
-        TensorProductPolynomials<dim, Polynomials::PiecewisePolynomial<double>>,
-        dim,
-        spacedim> *>(fe_ptr) != nullptr)
+  if (dynamic_cast<const FE_Poly<dim, spacedim> *>(fe_ptr) != nullptr)
     return true;
   if (dynamic_cast<const FE_DGP<dim, spacedim> *>(fe_ptr) != nullptr)
     return true;
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 542bf4d6b2..4749787ebf 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -136,20 +136,11 @@ namespace internal
         Assert(fe->n_components() == 1,
                ExcMessage("Expected a scalar element"));
 
-        const FE_Poly<TensorProductPolynomials<dim>, dim, dim> *fe_poly =
-          dynamic_cast<
-            const FE_Poly<TensorProductPolynomials<dim>, dim, dim> *>(fe);
-
-        const FE_Poly<
-          TensorProductPolynomials<dim,
-                                   Polynomials::PiecewisePolynomial<double>>,
-          dim,
-          dim> *fe_poly_piece =
-          dynamic_cast<const FE_Poly<
-            TensorProductPolynomials<dim,
-                                     Polynomials::PiecewisePolynomial<double>>,
-            dim,
-            dim> *>(fe);
+        const FE_Poly<dim, dim> *fe_poly =
+          dynamic_cast<const FE_Poly<dim, dim> *>(fe);
+
+        const FE_Poly<dim, dim> *fe_poly_piece =
+          dynamic_cast<const FE_Poly<dim, dim> *>(fe);
 
         const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim> *>(fe);
 
diff --git a/source/fe/fe_dgp.cc b/source/fe/fe_dgp.cc
index 208f8705cb..04fe320a73 100644
--- a/source/fe/fe_dgp.cc
+++ b/source/fe/fe_dgp.cc
@@ -27,7 +27,7 @@ DEAL_II_NAMESPACE_OPEN
 
 template <int dim, int spacedim>
 FE_DGP<dim, spacedim>::FE_DGP(const unsigned int degree)
-  : FE_Poly<PolynomialSpace<dim>, dim, spacedim>(
+  : FE_Poly<dim, spacedim>(
       PolynomialSpace<dim>(
         Polynomials::Legendre::generate_complete_basis(degree)),
       FiniteElementData<dim>(get_dpo_vector(degree),
diff --git a/source/fe/fe_dgp_monomial.cc b/source/fe/fe_dgp_monomial.cc
index 314950e33d..f505bb986f 100644
--- a/source/fe/fe_dgp_monomial.cc
+++ b/source/fe/fe_dgp_monomial.cc
@@ -130,7 +130,7 @@ namespace internal
 
 template <int dim>
 FE_DGPMonomial<dim>::FE_DGPMonomial(const unsigned int degree)
-  : FE_Poly<PolynomialsP<dim>, dim>(
+  : FE_Poly<dim>(
       PolynomialsP<dim>(degree),
       FiniteElementData<dim>(get_dpo_vector(degree),
                              1,
@@ -143,8 +143,8 @@ FE_DGPMonomial<dim>::FE_DGPMonomial(const unsigned int degree)
         FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
         std::vector<bool>(1, true)))
 {
-  Assert(this->poly_space.n() == this->dofs_per_cell, ExcInternalError());
-  Assert(this->poly_space.degree() == this->degree, ExcInternalError());
+  Assert(this->poly_space->n() == this->dofs_per_cell, ExcInternalError());
+  Assert(this->poly_space->degree() == this->degree, ExcInternalError());
 
   // DG doesn't have constraints, so
   // leave them empty
@@ -233,9 +233,13 @@ FE_DGPMonomial<dim>::get_interpolation_matrix(
           source_fe_matrix(k, j) = source_fe.shape_value(j, unit_points[k]);
 
       FullMatrix<double> this_matrix(this->dofs_per_cell, this->dofs_per_cell);
+      auto *const        polynomial_space_p =
+        dynamic_cast<PolynomialsP<dim> *>(this->poly_space.get());
+      Assert(polynomial_space_p != nullptr, ExcInternalError());
       for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
         for (unsigned int k = 0; k < unit_points.size(); ++k)
-          this_matrix(k, j) = this->poly_space.compute_value(j, unit_points[k]);
+          this_matrix(k, j) =
+            polynomial_space_p->compute_value(j, unit_points[k]);
 
       this_matrix.gauss_jordan();
 
@@ -454,8 +458,11 @@ FE_DGPMonomial<2>::has_support_on_face(const unsigned int shape_index,
     support_on_face = true;
   else
     {
+      auto *const polynomial_space_p =
+        dynamic_cast<PolynomialsP<2> *>(this->poly_space.get());
+      Assert(polynomial_space_p != nullptr, ExcInternalError());
       const std::array<unsigned int, 2> degrees =
-        this->poly_space.directional_degrees(shape_index);
+        polynomial_space_p->directional_degrees(shape_index);
 
       if ((face_index == 0 && degrees[1] == 0) ||
           (face_index == 3 && degrees[0] == 0))
@@ -476,8 +483,11 @@ FE_DGPMonomial<3>::has_support_on_face(const unsigned int shape_index,
     support_on_face = true;
   else
     {
+      auto *const polynomial_space_p =
+        dynamic_cast<PolynomialsP<3> *>(this->poly_space.get());
+      Assert(polynomial_space_p != nullptr, ExcInternalError());
       const std::array<unsigned int, 3> degrees =
-        this->poly_space.directional_degrees(shape_index);
+        polynomial_space_p->directional_degrees(shape_index);
 
       if ((face_index == 0 && degrees[1] == 0) ||
           (face_index == 2 && degrees[2] == 0) ||
diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc
index 818ddb0b6f..f54154595d 100644
--- a/source/fe/fe_dgq.cc
+++ b/source/fe/fe_dgq.cc
@@ -60,7 +60,7 @@ namespace internal
 
 template <int dim, int spacedim>
 FE_DGQ<dim, spacedim>::FE_DGQ(const unsigned int degree)
-  : FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>(
+  : FE_Poly<dim, spacedim>(
       TensorProductPolynomials<dim>(
         Polynomials::generate_complete_Lagrange_basis(
           internal::FE_DGQ::get_QGaussLobatto_points(degree))),
@@ -93,7 +93,7 @@ FE_DGQ<dim, spacedim>::FE_DGQ(const unsigned int degree)
 template <int dim, int spacedim>
 FE_DGQ<dim, spacedim>::FE_DGQ(
   const std::vector<Polynomials::Polynomial<double>> &polynomials)
-  : FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>(
+  : FE_Poly<dim, spacedim>(
       TensorProductPolynomials<dim>(polynomials),
       FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
                              1,
@@ -303,11 +303,14 @@ FE_DGQ<dim, spacedim>::get_interpolation_matrix(
       // cell and evaluate the
       // shape functions there
       const Point<dim> p = this->unit_support_points[j];
+      auto *const      polynomial_space_p =
+        dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+      Assert(polynomial_space_p != nullptr, ExcInternalError());
       for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
-        cell_interpolation(j, i) = this->poly_space.compute_value(i, p);
+        cell_interpolation(j, i) = polynomial_space_p->compute_value(i, p);
 
       for (unsigned int i = 0; i < source_fe.dofs_per_cell; ++i)
-        source_interpolation(j, i) = source_fe.poly_space.compute_value(i, p);
+        source_interpolation(j, i) = polynomial_space_p->compute_value(i, p);
     }
 
   // then compute the
@@ -832,8 +835,11 @@ FE_DGQArbitraryNodes<dim, spacedim>::get_name() const
   bool                equidistant = true;
   std::vector<double> points(this->degree + 1);
 
+  auto *const polynomial_space_p =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+  Assert(polynomial_space_p != nullptr, ExcInternalError());
   std::vector<unsigned int> lexicographic =
-    this->poly_space.get_numbering_inverse();
+    polynomial_space_p->get_numbering_inverse();
   for (unsigned int j = 0; j <= this->degree; j++)
     points[j] = this->unit_support_points[lexicographic[j]][0];
 
@@ -945,9 +951,12 @@ std::unique_ptr<FiniteElement<dim, spacedim>>
 FE_DGQArbitraryNodes<dim, spacedim>::clone() const
 {
   // Construct a dummy quadrature formula containing the FE's nodes:
-  std::vector<Point<1>>     qpoints(this->degree + 1);
+  std::vector<Point<1>> qpoints(this->degree + 1);
+  auto *const           polynomial_space_p =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+  Assert(polynomial_space_p != nullptr, ExcInternalError());
   std::vector<unsigned int> lexicographic =
-    this->poly_space.get_numbering_inverse();
+    polynomial_space_p->get_numbering_inverse();
   for (unsigned int i = 0; i <= this->degree; ++i)
     qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]);
   Quadrature<1> pquadrature(qpoints);
diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc
index c17fd45f7d..bee64acba8 100644
--- a/source/fe/fe_poly.cc
+++ b/source/fe/fe_poly.cc
@@ -29,139 +29,12 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-template <>
-void
-FE_Poly<TensorProductPolynomials<1>, 1, 2>::fill_fe_values(
-  const Triangulation<1, 2>::cell_iterator &,
-  const CellSimilarity::Similarity       cell_similarity,
-  const Quadrature<1> &                  quadrature,
-  const Mapping<1, 2> &                  mapping,
-  const Mapping<1, 2>::InternalDataBase &mapping_internal,
-  const dealii::internal::FEValuesImplementation::MappingRelatedData<1, 2>
-    &                                          mapping_data,
-  const FiniteElement<1, 2>::InternalDataBase &fe_internal,
-  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<1, 2>
-    &output_data) const
-{
-  // convert data object to internal data for this class. fails with an
-  // exception if that is not possible
-  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
-         ExcInternalError());
-  const InternalData &fe_data =
-    static_cast<const InternalData &>(fe_internal); // NOLINT
+template <int dim, int spacedim>
+FE_Poly<dim, spacedim>::FE_Poly(const FE_Poly &fe)
+  : FiniteElement<dim, spacedim>(fe)
+  , poly_space(fe.poly_space->clone())
+{}
 
-  const bool need_to_correct_higher_derivatives =
-    higher_derivatives_need_correcting(mapping,
-                                       mapping_data,
-                                       quadrature.size(),
-                                       fe_data.update_each);
-
-  // transform gradients and higher derivatives. there is nothing to do
-  // for values since we already emplaced them into output_data when
-  // we were in get_data()
-  if (fe_data.update_each & update_gradients &&
-      cell_similarity != CellSimilarity::translation)
-    for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-      mapping.transform(make_array_view(fe_data.shape_gradients, k),
-                        mapping_covariant,
-                        mapping_internal,
-                        make_array_view(output_data.shape_gradients, k));
-
-  if (fe_data.update_each & update_hessians &&
-      cell_similarity != CellSimilarity::translation)
-    {
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        mapping.transform(make_array_view(fe_data.shape_hessians, k),
-                          mapping_covariant_gradient,
-                          mapping_internal,
-                          make_array_view(output_data.shape_hessians, k));
-
-      if (need_to_correct_higher_derivatives)
-        correct_hessians(output_data, mapping_data, quadrature.size());
-    }
-
-  if (fe_data.update_each & update_3rd_derivatives &&
-      cell_similarity != CellSimilarity::translation)
-    {
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k),
-                          mapping_covariant_hessian,
-                          mapping_internal,
-                          make_array_view(output_data.shape_3rd_derivatives,
-                                          k));
-
-      if (need_to_correct_higher_derivatives)
-        correct_third_derivatives(output_data, mapping_data, quadrature.size());
-    }
-}
-
-
-
-template <>
-void
-FE_Poly<TensorProductPolynomials<2>, 2, 3>::fill_fe_values(
-  const Triangulation<2, 3>::cell_iterator &,
-  const CellSimilarity::Similarity       cell_similarity,
-  const Quadrature<2> &                  quadrature,
-  const Mapping<2, 3> &                  mapping,
-  const Mapping<2, 3>::InternalDataBase &mapping_internal,
-  const dealii::internal::FEValuesImplementation::MappingRelatedData<2, 3>
-    &                                          mapping_data,
-  const FiniteElement<2, 3>::InternalDataBase &fe_internal,
-  dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 3>
-    &output_data) const
-{
-  // convert data object to internal data for this class. fails with an
-  // exception if that is not possible
-  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
-         ExcInternalError());
-  const InternalData &fe_data =
-    static_cast<const InternalData &>(fe_internal); // NOLINT
-
-  const bool need_to_correct_higher_derivatives =
-    higher_derivatives_need_correcting(mapping,
-                                       mapping_data,
-                                       quadrature.size(),
-                                       fe_data.update_each);
-
-  // transform gradients and higher derivatives. there is nothing to do
-  // for values since we already emplaced them into output_data when
-  // we were in get_data()
-  if (fe_data.update_each & update_gradients &&
-      cell_similarity != CellSimilarity::translation)
-    for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-      mapping.transform(make_array_view(fe_data.shape_gradients, k),
-                        mapping_covariant,
-                        mapping_internal,
-                        make_array_view(output_data.shape_gradients, k));
-
-  if (fe_data.update_each & update_hessians &&
-      cell_similarity != CellSimilarity::translation)
-    {
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        mapping.transform(make_array_view(fe_data.shape_hessians, k),
-                          mapping_covariant_gradient,
-                          mapping_internal,
-                          make_array_view(output_data.shape_hessians, k));
-
-      if (need_to_correct_higher_derivatives)
-        correct_hessians(output_data, mapping_data, quadrature.size());
-    }
-
-  if (fe_data.update_each & update_3rd_derivatives &&
-      cell_similarity != CellSimilarity::translation)
-    {
-      for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
-        mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k),
-                          mapping_covariant_hessian,
-                          mapping_internal,
-                          make_array_view(output_data.shape_3rd_derivatives,
-                                          k));
-
-      if (need_to_correct_higher_derivatives)
-        correct_third_derivatives(output_data, mapping_data, quadrature.size());
-    }
-}
 
 
 #include "fe_poly.inst"
diff --git a/source/fe/fe_poly.inst.in b/source/fe/fe_poly.inst.in
index 486fbd92dd..dcaa87fc0d 100644
--- a/source/fe/fe_poly.inst.in
+++ b/source/fe/fe_poly.inst.in
@@ -18,28 +18,6 @@
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
-    template class FE_Poly<TensorProductPolynomials<deal_II_dimension>,
-                           deal_II_dimension,
-                           deal_II_space_dimension>;
-    template class FE_Poly<TensorProductPolynomialsConst<deal_II_dimension>,
-                           deal_II_dimension,
-                           deal_II_space_dimension>;
-    template class FE_Poly<TensorProductPolynomialsBubbles<deal_II_dimension>,
-                           deal_II_dimension,
-                           deal_II_space_dimension>;
-    template class FE_Poly<
-      TensorProductPolynomials<deal_II_dimension,
-                               Polynomials::PiecewisePolynomial<double>>,
-      deal_II_dimension,
-      deal_II_space_dimension>;
-    template class FE_Poly<PolynomialSpace<deal_II_dimension>,
-                           deal_II_dimension,
-                           deal_II_space_dimension>;
-    template class FE_Poly<PolynomialsP<deal_II_dimension>,
-                           deal_II_dimension,
-                           deal_II_space_dimension>;
-    template class FE_Poly<PolynomialsRannacherTurek<deal_II_dimension>,
-                           deal_II_dimension,
-                           deal_II_space_dimension>;
+    template class FE_Poly<deal_II_dimension, deal_II_space_dimension>;
 #endif
   }
diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc
index 006f9e1c4e..db825ba22d 100644
--- a/source/fe/fe_q.cc
+++ b/source/fe/fe_q.cc
@@ -99,8 +99,10 @@ FE_Q<dim, spacedim>::get_name() const
   std::vector<double> points(this->degree + 1);
 
   // Decode the support points in one coordinate direction.
+  TensorProductPolynomials<dim> *poly_space_derived_ptr =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
   std::vector<unsigned int> lexicographic =
-    this->poly_space.get_numbering_inverse();
+    poly_space_derived_ptr->get_numbering_inverse();
   for (unsigned int j = 0; j <= this->degree; j++)
     points[j] = this->unit_support_points[lexicographic[j]][0];
 
diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc
index 4ab6de7034..dc183f5151 100644
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -180,12 +180,14 @@ struct FE_Q_Base<PolynomialType, xdim, xspacedim>::Implementation
 
     // use that the element evaluates to 1 at index 0 and along the line at
     // zero
+    TensorProductPolynomials<dim> *poly_space_derived_ptr =
+      dynamic_cast<TensorProductPolynomials<dim> *>(fe.poly_space.get());
     const std::vector<unsigned int> &index_map_inverse =
-      fe.poly_space.get_numbering_inverse();
+      poly_space_derived_ptr->get_numbering_inverse();
     const std::vector<unsigned int> face_index_map =
       FETools::lexicographic_to_hierarchic_numbering<dim - 1>(q_deg);
     Assert(std::abs(
-             fe.poly_space.compute_value(index_map_inverse[0], Point<dim>()) -
+             fe.poly_space->compute_value(index_map_inverse[0], Point<dim>()) -
              1.) < 1e-14,
            ExcInternalError());
 
@@ -195,7 +197,7 @@ struct FE_Q_Base<PolynomialType, xdim, xspacedim>::Implementation
           Point<dim> p;
           p[0] = constraint_points[i](0);
           fe.interface_constraints(i, face_index_map[j]) =
-            fe.poly_space.compute_value(index_map_inverse[j], p);
+            fe.poly_space->compute_value(index_map_inverse[j], p);
 
           // if the value is small up to round-off, then simply set it to zero
           // to avoid unwanted fill-in of the constraint matrices (which would
@@ -315,12 +317,14 @@ struct FE_Q_Base<PolynomialType, xdim, xspacedim>::Implementation
 
     // use that the element evaluates to 1 at index 0 and along the line at
     // zero
+    TensorProductPolynomials<dim> *poly_space_derived_ptr =
+      dynamic_cast<TensorProductPolynomials<dim> *>(fe.poly_space.get());
     const std::vector<unsigned int> &index_map_inverse =
-      fe.poly_space.get_numbering_inverse();
+      poly_space_derived_ptr->get_numbering_inverse();
     const std::vector<unsigned int> face_index_map =
       FETools::lexicographic_to_hierarchic_numbering<dim - 1>(q_deg);
     Assert(std::abs(
-             fe.poly_space.compute_value(index_map_inverse[0], Point<dim>()) -
+             fe.poly_space->compute_value(index_map_inverse[0], Point<dim>()) -
              1.) < 1e-14,
            ExcInternalError());
 
@@ -385,8 +389,8 @@ struct FE_Q_Base<PolynomialType, xdim, xspacedim>::Implementation
               indices[1] * (q_deg + 1) + indices[0];
 
             fe.interface_constraints(i, face_index_map[j]) =
-              fe.poly_space.compute_value(index_map_inverse[new_index],
-                                          constraint_point);
+              fe.poly_space->compute_value(index_map_inverse[new_index],
+                                           constraint_point);
 
             // if the value is small up to round-off, then simply set it to
             // zero to avoid unwanted fill-in of the constraint matrices
@@ -407,7 +411,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::FE_Q_Base(
   const PolynomialType &        poly_space,
   const FiniteElementData<dim> &fe_data,
   const std::vector<bool> &     restriction_is_additive_flags)
-  : FE_Poly<PolynomialType, dim, spacedim>(
+  : FE_Poly<dim, spacedim>(
       poly_space,
       fe_data,
       restriction_is_additive_flags,
@@ -445,7 +449,9 @@ FE_Q_Base<PolynomialType, dim, spacedim>::initialize(
       FETools::hierarchic_to_lexicographic_numbering<dim>(q_degree);
     for (unsigned int i = q_dofs_per_cell; i < this->dofs_per_cell; ++i)
       renumber.push_back(i);
-    this->poly_space.set_numbering(renumber);
+    TensorProductPolynomials<dim> *poly_space_derived_ptr =
+      dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+    poly_space_derived_ptr->set_numbering(renumber);
   }
 
   // Finally fill in support points on cell and face and initialize
@@ -502,12 +508,12 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_interpolation_matrix(
 
           // FE_Q element evaluates to 1 in unit support point and to zero in
           // all other points by construction
-          Assert(std::abs(this->poly_space.compute_value(j, p) - 1.) < 1e-13,
+          Assert(std::abs(this->poly_space->compute_value(j, p) - 1.) < 1e-13,
                  ExcInternalError());
 
           for (unsigned int i = 0; i < source_q_dofs_per_cell; ++i)
             interpolation_matrix(j, i) =
-              source_fe->poly_space.compute_value(i, p);
+              source_fe->poly_space->compute_value(i, p);
         }
 
       // for FE_Q_DG0, add one last row of identity
@@ -748,10 +754,15 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_line_dof_identities(
 
       std::vector<std::pair<unsigned int, unsigned int>> identities;
 
+      TensorProductPolynomials<dim> *poly_space_derived_ptr =
+        dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
       const std::vector<unsigned int> &index_map_inverse =
-        this->poly_space.get_numbering_inverse();
+        poly_space_derived_ptr->get_numbering_inverse();
+      TensorProductPolynomials<dim> *poly_space_derived_ptr_other =
+        dynamic_cast<TensorProductPolynomials<dim> *>(
+          fe_q_other->poly_space.get());
       const std::vector<unsigned int> &index_map_inverse_other =
-        fe_q_other->poly_space.get_numbering_inverse();
+        poly_space_derived_ptr_other->get_numbering_inverse();
 
       for (unsigned int i = 0; i < p - 1; ++i)
         for (unsigned int j = 0; j < q - 1; ++j)
@@ -810,10 +821,15 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_quad_dof_identities(
 
       std::vector<std::pair<unsigned int, unsigned int>> identities;
 
+      TensorProductPolynomials<dim> *poly_space_derived_ptr =
+        dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
       const std::vector<unsigned int> &index_map_inverse =
-        this->poly_space.get_numbering_inverse();
+        poly_space_derived_ptr->get_numbering_inverse();
+      TensorProductPolynomials<dim> *poly_space_derived_ptr_other =
+        dynamic_cast<TensorProductPolynomials<dim> *>(
+          fe_q_other->poly_space.get());
       const std::vector<unsigned int> &index_map_inverse_other =
-        fe_q_other->poly_space.get_numbering_inverse();
+        poly_space_derived_ptr_other->get_numbering_inverse();
 
       for (unsigned int i1 = 0; i1 < p - 1; ++i1)
         for (unsigned int i2 = 0; i2 < p - 1; ++i2)
@@ -870,8 +886,10 @@ void
 FE_Q_Base<PolynomialType, dim, spacedim>::initialize_unit_support_points(
   const std::vector<Point<1>> &points)
 {
+  TensorProductPolynomials<dim> *poly_space_derived_ptr =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
   const std::vector<unsigned int> &index_map_inverse =
-    this->poly_space.get_numbering_inverse();
+    poly_space_derived_ptr->get_numbering_inverse();
 
   // We can compute the support points by computing the tensor
   // product of the 1d set of points. We could do this by hand, but it's
@@ -1175,7 +1193,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_prolongation_matrix(
       // current numbering gives the identity operation
       for (unsigned int i = 0; i < q_dofs_per_cell; ++i)
         {
-          Assert(std::fabs(1. - this->poly_space.compute_value(
+          Assert(std::fabs(1. - this->poly_space->compute_value(
                                   i, this->unit_support_points[i])) < eps,
                  ExcInternalError("The Lagrange polynomial does not evaluate "
                                   "to one or zero in a nodal point. "
@@ -1185,7 +1203,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_prolongation_matrix(
                                   "prevents the sum to be one."));
           for (unsigned int j = 0; j < q_dofs_per_cell; ++j)
             if (j != i)
-              Assert(std::fabs(this->poly_space.compute_value(
+              Assert(std::fabs(this->poly_space->compute_value(
                        i, this->unit_support_points[j])) < eps,
                      ExcInternalError(
                        "The Lagrange polynomial does not evaluate "
@@ -1204,8 +1222,11 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_prolongation_matrix(
       const unsigned int            dofs1d = q_degree + 1;
       std::vector<Table<2, double>> subcell_evaluations(
         dim, Table<2, double>(dofs1d, dofs1d));
+
+      TensorProductPolynomials<dim> *poly_space_derived_ptr =
+        dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
       const std::vector<unsigned int> &index_map_inverse =
-        this->poly_space.get_numbering_inverse();
+        poly_space_derived_ptr->get_numbering_inverse();
 
       // helper value: step size how to walk through diagonal and how many
       // points we have left apart from the first dimension
@@ -1238,7 +1259,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_prolongation_matrix(
                 Point<dim> point;
                 point[0] = p_cell[d];
                 const double cell_value =
-                  this->poly_space.compute_value(index_map_inverse[i], point);
+                  this->poly_space->compute_value(index_map_inverse[i], point);
 
                 // cut off values that are too small. note that we have here
                 // Lagrange interpolation functions, so they should be zero at
@@ -1280,7 +1301,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_prolongation_matrix(
                   subcell_evaluations[d](j_indices[d - 1], i_indices[d - 1]);
 
               // innermost sum where we actually compute. the same as
-              // prolongate(j,i) = this->poly_space.compute_value (i, p_cell)
+              // prolongate(j,i) = this->poly_space->compute_value (i, p_cell)
               for (unsigned int jj = 0; jj < dofs1d; ++jj)
                 {
                   const unsigned int j_ind = index_map_inverse[j + jj];
@@ -1379,9 +1400,12 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_restriction_matrix(
       // assumption that whenever a row makes a non-zero contribution to the
       // mother's residual, the correct value is interpolated.
 
-      const double                     eps = 1e-15 * q_degree * dim;
+      const double eps = 1e-15 * q_degree * dim;
+
+      TensorProductPolynomials<dim> *poly_space_derived_ptr =
+        dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
       const std::vector<unsigned int> &index_map_inverse =
-        this->poly_space.get_numbering_inverse();
+        poly_space_derived_ptr->get_numbering_inverse();
 
       const unsigned int          dofs1d = q_degree + 1;
       std::vector<Tensor<1, dim>> evaluations1d(dofs1d);
@@ -1411,8 +1435,8 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_restriction_matrix(
                     Point<dim> point;
                     point[0] = p_subcell[d];
                     evaluations1d[j][d] =
-                      this->poly_space.compute_value(index_map_inverse[j],
-                                                     point);
+                      this->poly_space->compute_value(index_map_inverse[j],
+                                                      point);
                   }
               unsigned int j_indices[dim];
               internal::FE_Q_Base::zero_indices<dim>(j_indices);
diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc
index d77f389ad5..41e7021282 100644
--- a/source/fe/fe_q_hierarchical.cc
+++ b/source/fe/fe_q_hierarchical.cc
@@ -58,8 +58,9 @@ namespace internal
 
 template <int dim>
 FE_Q_Hierarchical<dim>::FE_Q_Hierarchical(const unsigned int degree)
-  : FE_Poly<TensorProductPolynomials<dim>, dim>(
-      Polynomials::Hierarchical::generate_complete_basis(degree),
+  : FE_Poly<dim>(
+      TensorProductPolynomials<dim>(
+        Polynomials::Hierarchical::generate_complete_basis(degree)),
       FiniteElementData<dim>(get_dpo_vector(degree),
                              1,
                              degree,
@@ -72,7 +73,9 @@ FE_Q_Hierarchical<dim>::FE_Q_Hierarchical(const unsigned int degree)
         std::vector<bool>(1, true)))
   , face_renumber(face_fe_q_hierarchical_to_hierarchic_numbering(degree))
 {
-  this->poly_space.set_numbering(
+  TensorProductPolynomials<dim> *poly_space_derived_ptr =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+  poly_space_derived_ptr->set_numbering(
     hierarchic_to_fe_q_hierarchical_numbering(*this));
 
   // The matrix @p{dofs_cell} contains the
@@ -673,7 +676,10 @@ FE_Q_Hierarchical<dim>::initialize_embedding_and_restriction(
   unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
 
   const unsigned int dofs_1d = 2 * this->dofs_per_vertex + this->dofs_per_line;
-  const std::vector<unsigned int> &renumber = this->poly_space.get_numbering();
+  TensorProductPolynomials<dim> *poly_space_derived_ptr =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+  const std::vector<unsigned int> &renumber =
+    poly_space_derived_ptr->get_numbering();
 
   for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; ++c)
     {
@@ -810,8 +816,10 @@ FE_Q_Hierarchical<dim>::initialize_generalized_support_points()
 
   this->generalized_support_points.resize(n);
 
+  TensorProductPolynomials<dim> *poly_space_derived_ptr =
+    dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
   const std::vector<unsigned int> &index_map_inverse =
-    this->poly_space.get_numbering_inverse();
+    poly_space_derived_ptr->get_numbering_inverse();
 
   Point<dim> p;
   // the method of numbering allows
diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc
index 24d4a5714b..dded3b04f5 100644
--- a/source/fe/fe_rannacher_turek.cc
+++ b/source/fe/fe_rannacher_turek.cc
@@ -32,14 +32,13 @@ template <int dim>
 FE_RannacherTurek<dim>::FE_RannacherTurek(
   const unsigned int order,
   const unsigned int n_face_support_points)
-  : FE_Poly<PolynomialsRannacherTurek<dim>, dim>(
-      PolynomialsRannacherTurek<dim>(),
-      FiniteElementData<dim>(this->get_dpo_vector(),
-                             1,
-                             2,
-                             FiniteElementData<dim>::L2),
-      std::vector<bool>(4, false), // restriction not implemented
-      std::vector<ComponentMask>(4, std::vector<bool>(1, true)))
+  : FE_Poly<dim>(PolynomialsRannacherTurek<dim>(),
+                 FiniteElementData<dim>(this->get_dpo_vector(),
+                                        1,
+                                        2,
+                                        FiniteElementData<dim>::L2),
+                 std::vector<bool>(4, false), // restriction not implemented
+                 std::vector<ComponentMask>(4, std::vector<bool>(1, true)))
   , order(order)
   , n_face_support_points(n_face_support_points)
 {
-- 
2.39.5