Change TensorProductPolynomialsBubbles to no longer derive from TensorProductPolynomials

diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h
index 9c39dec516942ed5dafaa68046c37907f5290baf..3fd35d075ca0ca32c5532b4341f74517e6d23581 100644 (file)
--- a/include/deal.II/base/tensor_product_polynomials.h
+++ b/include/deal.II/base/tensor_product_polynomials.h
@@ -29,6 +29,9 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+template <int dim>
+class TensorProductPolynomialsBubbles;
+
 /**
  * @addtogroup Polynomials
  * @{
@@ -227,6 +230,12 @@ protected:
   void
   compute_index(const unsigned int i,
                 unsigned int (&indices)[(dim > 0 ? dim : 1)]) const;
+
+  /**
+   * TensorProductPolynomialsBubbles has a TensorProductPolynomials class
+   * so we declare it as a friend class.
+   */
+  friend class TensorProductPolynomialsBubbles<dim>;
 };
 
 

diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h
index 7a0f3dd0187e5fec7f5b73f9be835952026e27eb..af66095a6d57aef4776d0dc0bb68c5443f9f1e24 100644 (file)
--- a/include/deal.II/base/tensor_product_polynomials_bubbles.h
+++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h
@@ -30,7 +30,6 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-
 /**
  * @addtogroup Polynomials
  * @{
@@ -45,9 +44,15 @@ DEAL_II_NAMESPACE_OPEN
  * @author Daniel Arndt, 2015
  */
 template <int dim>
-class TensorProductPolynomialsBubbles : public TensorProductPolynomials<dim>
+class TensorProductPolynomialsBubbles
 {
 public:
+  /**
+   * Access to the dimension of this object, for checking and automatic
+   * setting of dimension in other classes.
+   */
+  static const unsigned int dimension = dim;
+
   /**
    * Constructor. <tt>pols</tt> is a vector of objects that should be derived
    * or otherwise convertible to one-dimensional polynomial objects. It will
@@ -56,6 +61,32 @@ public:
   template <class Pol>
   TensorProductPolynomialsBubbles(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>.
@@ -145,6 +176,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;
 };
 
 /** @} */
@@ -158,29 +205,31 @@ template <int dim>
 template <class Pol>
 inline TensorProductPolynomialsBubbles<dim>::TensorProductPolynomialsBubbles(
   const std::vector<Pol> &pols)
-  : TensorProductPolynomials<dim>(pols)
+  : tensor_polys(pols)
+  , index_map(tensor_polys.n() +
+              ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
+  , index_map_inverse(tensor_polys.n() +
+                      ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
 {
-  const unsigned int q_degree  = this->polynomials.size() - 1;
+  const unsigned int q_degree  = tensor_polys.polynomials.size() - 1;
   const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
   // append index for renumbering
-  for (unsigned int i = 0; i < n_bubbles; ++i)
+  for (unsigned int i = 0; i < tensor_polys.n() + n_bubbles; ++i)
     {
-      this->index_map.push_back(i + this->n_tensor_pols);
-      this->index_map_inverse.push_back(i + this->n_tensor_pols);
+      index_map[i]         = i;
+      index_map_inverse[i] = i;
     }
 }
 
 
-
 template <int dim>
 inline unsigned int
 TensorProductPolynomialsBubbles<dim>::n() const
 {
-  return this->n_tensor_pols + dim;
+  return tensor_polys.n() + dim;
 }
 
 
-
 template <>
 inline unsigned int
 TensorProductPolynomialsBubbles<0>::n() const
@@ -188,6 +237,23 @@ TensorProductPolynomialsBubbles<0>::n() const
   return numbers::invalid_unsigned_int;
 }
 
+
+template <int dim>
+inline const std::vector<unsigned int> &
+TensorProductPolynomialsBubbles<dim>::get_numbering() const
+{
+  return index_map;
+}
+
+
+template <int dim>
+inline const std::vector<unsigned int> &
+TensorProductPolynomialsBubbles<dim>::get_numbering_inverse() const
+{
+  return index_map_inverse;
+}
+
+
 template <int dim>
 template <int order>
 Tensor<order, dim>
@@ -195,18 +261,17 @@ TensorProductPolynomialsBubbles<dim>::compute_derivative(
   const unsigned int i,
   const Point<dim> & p) const
 {
-  const unsigned int q_degree      = this->polynomials.size() - 1;
-  const unsigned int max_q_indices = this->n_tensor_pols;
+  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
   if (i < max_q_indices)
-    return this
-      ->TensorProductPolynomials<dim>::template compute_derivative<order>(i, p);
+    return tensor_polys.template compute_derivative<order>(i, p);
 
-  const unsigned int comp = i - this->n_tensor_pols;
+  const unsigned int comp = i - tensor_polys.n();
 
   Tensor<order, dim> derivative;
   switch (order)

diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc
index afd88bfb827462945fc7ca94629d860374bc2f4f..dd231f8304870fa028bd3dd2d89486556b1c3236 100644 (file)
--- a/source/base/tensor_product_polynomials_bubbles.cc
+++ b/source/base/tensor_product_polynomials_bubbles.cc
@@ -24,22 +24,60 @@ DEAL_II_NAMESPACE_OPEN
 /* ------------------- TensorProductPolynomialsBubbles -------------- */
 
 
+
+template <int dim>
+void
+TensorProductPolynomialsBubbles<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
+TensorProductPolynomialsBubbles<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
 TensorProductPolynomialsBubbles<dim>::compute_value(const unsigned int i,
                                                     const Point<dim> & p) const
 {
-  const unsigned int q_degree      = this->polynomials.size() - 1;
-  const unsigned int max_q_indices = this->n_tensor_pols;
+  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
   if (i < max_q_indices)
-    return this->TensorProductPolynomials<dim>::compute_value(i, p);
+    return tensor_polys.compute_value(i, p);
 
-  const unsigned int comp = i - this->n_tensor_pols;
+  const unsigned int comp = i - tensor_polys.n();
 
   // compute \prod_{i=1}^d 4*(1-x_i^2)(p)
   double value = 1.;
@@ -63,22 +101,23 @@ TensorProductPolynomialsBubbles<0>::compute_value(const unsigned int,
 }
 
 
+
 template <int dim>
 Tensor<1, dim>
 TensorProductPolynomialsBubbles<dim>::compute_grad(const unsigned int i,
                                                    const Point<dim> & p) const
 {
-  const unsigned int q_degree      = this->polynomials.size() - 1;
-  const unsigned int max_q_indices = this->n_tensor_pols;
+  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
   if (i < max_q_indices)
-    return this->TensorProductPolynomials<dim>::compute_grad(i, p);
+    return tensor_polys.compute_grad(i, p);
 
-  const unsigned int comp = i - this->n_tensor_pols;
+  const unsigned int comp = i - tensor_polys.n();
   Tensor<1, dim>     grad;
 
   for (unsigned int d = 0; d < dim; ++d)
@@ -116,17 +155,17 @@ TensorProductPolynomialsBubbles<dim>::compute_grad_grad(
   const unsigned int i,
   const Point<dim> & p) const
 {
-  const unsigned int q_degree      = this->polynomials.size() - 1;
-  const unsigned int max_q_indices = this->n_tensor_pols;
+  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
   if (i < max_q_indices)
-    return this->TensorProductPolynomials<dim>::compute_grad_grad(i, p);
+    return tensor_polys.compute_grad_grad(i, p);
 
-  const unsigned int comp = i - this->n_tensor_pols;
+  const unsigned int comp = i - tensor_polys.n();
 
   double v[dim + 1][3];
   {
@@ -213,6 +252,8 @@ TensorProductPolynomialsBubbles<dim>::compute_grad_grad(
   return grad_grad;
 }
 
+
+
 template <int dim>
 void
 TensorProductPolynomialsBubbles<dim>::evaluate(
@@ -223,8 +264,8 @@ TensorProductPolynomialsBubbles<dim>::evaluate(
   std::vector<Tensor<3, dim>> &third_derivatives,
   std::vector<Tensor<4, dim>> &fourth_derivatives) const
 {
-  const unsigned int q_degree      = this->polynomials.size() - 1;
-  const unsigned int max_q_indices = this->n_tensor_pols;
+  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);
   Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
@@ -249,36 +290,34 @@ TensorProductPolynomialsBubbles<dim>::evaluate(
   bool do_3rd_derivatives = false, do_4th_derivatives = false;
   if (values.empty() == false)
     {
-      values.resize(this->n_tensor_pols);
+      values.resize(tensor_polys.n());
       do_values = true;
     }
   if (grads.empty() == false)
     {
-      grads.resize(this->n_tensor_pols);
+      grads.resize(tensor_polys.n());
       do_grads = true;
     }
   if (grad_grads.empty() == false)
     {
-      grad_grads.resize(this->n_tensor_pols);
+      grad_grads.resize(tensor_polys.n());
       do_grad_grads = true;
     }
   if (third_derivatives.empty() == false)
     {
-      third_derivatives.resize(this->n_tensor_pols);
+      third_derivatives.resize(tensor_polys.n());
       do_3rd_derivatives = true;
     }
   if (fourth_derivatives.empty() == false)
     {
-      fourth_derivatives.resize(this->n_tensor_pols);
+      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 (unsigned int i = this->n_tensor_pols;
-       i < this->n_tensor_pols + n_bubbles;
-       ++i)
+  for (unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i)
     {
       if (do_values)
         values.push_back(compute_value(i, p));