Change TensorProductPolynomialsConst to no longer derive from TensorProductPolynomials

diff --git a/include/deal.II/base/tensor_product_polynomials_const.h b/include/deal.II/base/tensor_product_polynomials_const.h
index a9b1e131f80d5e5ce94b0dc96a80b80176ada21b..0f8cb93d0c302c8447e88e5b8e4779715e8c8930 100644 (file)
--- a/include/deal.II/base/tensor_product_polynomials_const.h
+++ b/include/deal.II/base/tensor_product_polynomials_const.h
@@ -45,7 +45,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Timo Heister, 2012
  */
 template <int dim>
-class TensorProductPolynomialsConst : public TensorProductPolynomials<dim>
+class TensorProductPolynomialsConst
 {
 public:
   /**
@@ -62,6 +62,32 @@ public:
   template <class Pol>
   TensorProductPolynomialsConst(const std::vector<Pol> &pols);
 
+  /**
+   * Print the list of <tt>tensor_polys</tt> indices to <tt>out</tt>.
+   */
+  void
+  output_indices(std::ostream &out) const;
+
+  /**
+   * Set the ordering of the polynomials. Requires
+   * <tt>renumber.size()==tensor_polys.n()</tt>.  Stores a copy of
+   * <tt>renumber</tt>.
+   */
+  void
+  set_numbering(const std::vector<unsigned int> &renumber);
+
+  /**
+   * Give read access to the renumber vector.
+   */
+  const std::vector<unsigned int> &
+  get_numbering() const;
+
+  /**
+   * Give read access to the inverse renumber vector.
+   */
+  const std::vector<unsigned int> &
+  get_numbering_inverse() const;
+
   /**
    * Compute the value and the first and second derivatives of each tensor
    * product polynomial at <tt>unit_point</tt>.
@@ -151,6 +177,22 @@ public:
    */
   unsigned int
   n() const;
+
+private:
+  /**
+   * The TensorProductPolynomials object
+   */
+  TensorProductPolynomials<dim> tensor_polys;
+
+  /**
+   * Index map for reordering the polynomials.
+   */
+  std::vector<unsigned int> index_map;
+
+  /**
+   * Index map for reordering the polynomials.
+   */
+  std::vector<unsigned int> index_map_inverse;
 };
 
 /** @} */
@@ -164,12 +206,10 @@ template <int dim>
 template <class Pol>
 inline TensorProductPolynomialsConst<dim>::TensorProductPolynomialsConst(
   const std::vector<Pol> &pols)
-  : TensorProductPolynomials<dim>(pols)
-{
-  // append index for renumbering
-  this->index_map.push_back(TensorProductPolynomials<dim>::n());
-  this->index_map_inverse.push_back(TensorProductPolynomials<dim>::n());
-}
+  : tensor_polys(pols)
+  , index_map(tensor_polys.n() + 1)
+  , index_map_inverse(tensor_polys.n() + 1)
+{}
 
 
 
@@ -177,7 +217,24 @@ template <int dim>
 inline unsigned int
 TensorProductPolynomialsConst<dim>::n() const
 {
-  return TensorProductPolynomials<dim>::n() + 1;
+  return tensor_polys.n() + 1;
+}
+
+
+
+template <int dim>
+inline const std::vector<unsigned int> &
+TensorProductPolynomialsConst<dim>::get_numbering() const
+{
+  return index_map;
+}
+
+
+template <int dim>
+inline const std::vector<unsigned int> &
+TensorProductPolynomialsConst<dim>::get_numbering_inverse() const
+{
+  return index_map_inverse;
 }
 
 
@@ -196,13 +253,12 @@ TensorProductPolynomialsConst<dim>::compute_derivative(
   const unsigned int i,
   const Point<dim> & p) const
 {
-  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
+  const unsigned int max_indices = tensor_polys.n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
   if (i < max_indices)
-    return this
-      ->TensorProductPolynomials<dim>::template compute_derivative<order>(i, p);
+    return tensor_polys.template compute_derivative<order>(i, p);
   else
     // this is for the constant function
     return Tensor<order, dim>();

diff --git a/source/base/tensor_product_polynomials_const.cc b/source/base/tensor_product_polynomials_const.cc
index 54286ba84cf35c4d4a4fcc432c2af862cd48917f..c9b988073d4b74ff5fb6ce768afb886de57b8ad8 100644 (file)
--- a/source/base/tensor_product_polynomials_const.cc
+++ b/source/base/tensor_product_polynomials_const.cc
@@ -24,17 +24,55 @@ DEAL_II_NAMESPACE_OPEN
 /* ------------------- TensorProductPolynomialsConst -------------- */
 
 
+
+template <int dim>
+void
+TensorProductPolynomialsConst<dim>::output_indices(std::ostream &out) const
+{
+  unsigned int ix[dim];
+  for (unsigned int i = 0; i < tensor_polys.n(); ++i)
+    {
+      tensor_polys.compute_index(i, ix);
+      out << i << "\t";
+      for (unsigned int d = 0; d < dim; ++d)
+        out << ix[d] << " ";
+      out << std::endl;
+    }
+}
+
+
+
+template <int dim>
+void
+TensorProductPolynomialsConst<dim>::set_numbering(
+  const std::vector<unsigned int> &renumber)
+{
+  Assert(renumber.size() == index_map.size(),
+         ExcDimensionMismatch(renumber.size(), index_map.size()));
+
+  index_map = renumber;
+  for (unsigned int i = 0; i < index_map.size(); ++i)
+    index_map_inverse[index_map[i]] = i;
+
+  std::vector<unsigned int> renumber_base;
+  for (unsigned int i = 0; i < tensor_polys.n(); ++i)
+    renumber_base.push_back(renumber[i]);
+
+  tensor_polys.set_numbering(renumber_base);
+}
+
+
 template <int dim>
 double
 TensorProductPolynomialsConst<dim>::compute_value(const unsigned int i,
                                                   const Point<dim> & p) const
 {
-  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
+  const unsigned int max_indices = tensor_polys.n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
   if (i < max_indices)
-    return this->TensorProductPolynomials<dim>::compute_value(i, p);
+    return tensor_polys.compute_value(i, p);
   else
     // this is for the constant function
     return 1.;
@@ -57,12 +95,12 @@ Tensor<1, dim>
 TensorProductPolynomialsConst<dim>::compute_grad(const unsigned int i,
                                                  const Point<dim> & p) const
 {
-  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
+  const unsigned int max_indices = tensor_polys.n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
   if (i < max_indices)
-    return this->TensorProductPolynomials<dim>::compute_grad(i, p);
+    return tensor_polys.compute_grad(i, p);
   else
     // this is for the constant function
     return Tensor<1, dim>();
@@ -73,12 +111,12 @@ Tensor<2, dim>
 TensorProductPolynomialsConst<dim>::compute_grad_grad(const unsigned int i,
                                                       const Point<dim> &p) const
 {
-  const unsigned int max_indices = TensorProductPolynomials<dim>::n();
+  const unsigned int max_indices = tensor_polys.n();
   Assert(i <= max_indices, ExcInternalError());
 
   // treat the regular basis functions
   if (i < max_indices)
-    return this->TensorProductPolynomials<dim>::compute_grad_grad(i, p);
+    return tensor_polys.compute_grad_grad(i, p);
   else
     // this is for the constant function
     return Tensor<2, dim>();
@@ -94,30 +132,21 @@ TensorProductPolynomialsConst<dim>::evaluate(
   std::vector<Tensor<3, dim>> &third_derivatives,
   std::vector<Tensor<4, dim>> &fourth_derivatives) const
 {
-  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 ||
+  Assert(values.size() == tensor_polys.n() + 1 || values.size() == 0,
+         ExcDimensionMismatch2(values.size(), tensor_polys.n() + 1, 0));
+  Assert(grads.size() == tensor_polys.n() + 1 || grads.size() == 0,
+         ExcDimensionMismatch2(grads.size(), tensor_polys.n() + 1, 0));
+  Assert(grad_grads.size() == tensor_polys.n() + 1 || grad_grads.size() == 0,
+         ExcDimensionMismatch2(grad_grads.size(), tensor_polys.n() + 1, 0));
+  Assert(third_derivatives.size() == tensor_polys.n() + 1 ||
            third_derivatives.size() == 0,
          ExcDimensionMismatch2(third_derivatives.size(),
-                               TensorProductPolynomials<dim>::n() + 1,
+                               tensor_polys.n() + 1,
                                0));
-  Assert(fourth_derivatives.size() == TensorProductPolynomials<dim>::n() + 1 ||
+  Assert(fourth_derivatives.size() == tensor_polys.n() + 1 ||
            fourth_derivatives.size() == 0,
          ExcDimensionMismatch2(fourth_derivatives.size(),
-                               TensorProductPolynomials<dim>::n() + 1,
+                               tensor_polys.n() + 1,
                                0));
 
   // remove slot for const value, go into the base class compute method and
@@ -141,16 +170,16 @@ TensorProductPolynomialsConst<dim>::evaluate(
     }
   if (third_derivatives.empty() == false)
     {
-      third_derivatives.resize(TensorProductPolynomials<dim>::n());
+      third_derivatives.resize(tensor_polys.n());
       do_3rd_derivatives = true;
     }
   if (fourth_derivatives.empty() == false)
     {
-      fourth_derivatives.resize(TensorProductPolynomials<dim>::n());
+      fourth_derivatives.resize(tensor_polys.n());
       do_4th_derivatives = true;
     }
 
-  this->TensorProductPolynomials<dim>::evaluate(
+  tensor_polys.evaluate(
     p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
 
   // for dgq node: values =1, grads=0, grads_grads=0, third_derivatives=0,