Fix numbering

diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h
index fea8899e25af36a1a2814e1e6e37b2baaff64bde..2f2e1776b7558086a303e98e8288efeb16f2ade8 100644 (file)
--- a/include/deal.II/fe/fe_poly.templates.h
+++ b/include/deal.II/fe/fe_poly.templates.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/polynomial_space.h>
+#include <deal.II/base/polynomials_piecewise.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 #include <deal.II/base/tensor_product_polynomials_bubbles.h>
@@ -579,13 +580,23 @@ FE_Poly<dim, spacedim>::get_poly_space_numbering() const
 {
   auto *const space_tensor_prod =
     dynamic_cast<TensorProductPolynomials<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> *>(this->poly_space.get());
+  auto *const space_tensor_prod_piecewise = dynamic_cast<
+    TensorProductPolynomials<dim, Polynomials::PiecewisePolynomial<double>> *>(
+    this->poly_space.get());
+  if (space_tensor_prod_piecewise != nullptr)
+    return space_tensor_prod_piecewise->get_numbering();
 
+  auto *const space_tensor_prod_bubbles =
+    dynamic_cast<TensorProductPolynomialsBubbles<dim> *>(
+      this->poly_space.get());
+  if (space_tensor_prod_bubbles != nullptr)
+    return space_tensor_prod_bubbles->get_numbering();
+
+  auto *const space_tensor_prod_const =
+    dynamic_cast<TensorProductPolynomialsConst<dim> *>(this->poly_space.get());
   if (space_tensor_prod_const != nullptr)
     return space_tensor_prod_const->get_numbering();
 

diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc
index f54154595d5992cce08fc2e2f348f1925042a7cc..ad26b3049daa0538ba3f2e1d1e4813f764ccdde5 100644 (file)
--- a/source/fe/fe_dgq.cc
+++ b/source/fe/fe_dgq.cc
@@ -835,11 +835,11 @@ FE_DGQArbitraryNodes<dim, spacedim>::get_name() const
   bool                equidistant = true;
   std::vector<double> points(this->degree + 1);
 
-  auto *const polynomial_space_p =
+  auto *const polynomial_space =
     dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
-  Assert(polynomial_space_p != nullptr, ExcInternalError());
+  Assert(polynomial_space != nullptr, ExcInternalError());
   std::vector<unsigned int> lexicographic =
-    polynomial_space_p->get_numbering_inverse();
+    polynomial_space->get_numbering_inverse();
   for (unsigned int j = 0; j <= this->degree; j++)
     points[j] = this->unit_support_points[lexicographic[j]][0];
 
@@ -952,11 +952,11 @@ FE_DGQArbitraryNodes<dim, spacedim>::clone() const
 {
   // Construct a dummy quadrature formula containing the FE's nodes:
   std::vector<Point<1>> qpoints(this->degree + 1);
-  auto *const           polynomial_space_p =
+  auto *const           polynomial_space =
     dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
-  Assert(polynomial_space_p != nullptr, ExcInternalError());
+  Assert(polynomial_space != nullptr, ExcInternalError());
   std::vector<unsigned int> lexicographic =
-    polynomial_space_p->get_numbering_inverse();
+    polynomial_space->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_q_base.cc b/source/fe/fe_q_base.cc
index dc183f51510103e444c9fb6ff51e4118f07a6767..1dff72f2b480204b33f2f18629b4e4bbd4c925fa 100644 (file)
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -180,10 +180,8 @@ 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 =
-      poly_space_derived_ptr->get_numbering_inverse();
+      fe.get_poly_space_numbering_inverse();
     const std::vector<unsigned int> face_index_map =
       FETools::lexicographic_to_hierarchic_numbering<dim - 1>(q_deg);
     Assert(std::abs(
@@ -317,10 +315,8 @@ 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 =
-      poly_space_derived_ptr->get_numbering_inverse();
+      fe.get_poly_space_numbering_inverse();
     const std::vector<unsigned int> face_index_map =
       FETools::lexicographic_to_hierarchic_numbering<dim - 1>(q_deg);
     Assert(std::abs(
@@ -444,15 +440,44 @@ FE_Q_Base<PolynomialType, dim, spacedim>::initialize(
            q_dofs_per_cell + dim == this->dofs_per_cell,
          ExcInternalError());
 
-  {
+  [this, q_dofs_per_cell]() {
     std::vector<unsigned int> renumber =
       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);
-    TensorProductPolynomials<dim> *poly_space_derived_ptr =
+    auto *tensor_poly_space_ptr =
       dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
-    poly_space_derived_ptr->set_numbering(renumber);
-  }
+    if (tensor_poly_space_ptr != nullptr)
+      {
+        tensor_poly_space_ptr->set_numbering(renumber);
+        return;
+      }
+    auto *tensor_piecewise_poly_space_ptr = dynamic_cast<
+      TensorProductPolynomials<dim, Polynomials::PiecewisePolynomial<double>>
+        *>(this->poly_space.get());
+    if (tensor_piecewise_poly_space_ptr != nullptr)
+      {
+        tensor_piecewise_poly_space_ptr->set_numbering(renumber);
+        return;
+      }
+    auto *tensor_bubbles_poly_space_ptr =
+      dynamic_cast<TensorProductPolynomialsBubbles<dim> *>(
+        this->poly_space.get());
+    if (tensor_bubbles_poly_space_ptr != nullptr)
+      {
+        tensor_bubbles_poly_space_ptr->set_numbering(renumber);
+        return;
+      }
+    auto *tensor_const_poly_space_ptr =
+      dynamic_cast<TensorProductPolynomialsConst<dim> *>(
+        this->poly_space.get());
+    if (tensor_const_poly_space_ptr != nullptr)
+      {
+        tensor_const_poly_space_ptr->set_numbering(renumber);
+        return;
+      }
+    Assert(false, ExcNotImplemented());
+  }();
 
   // Finally fill in support points on cell and face and initialize
   // constraints. All of this can happen in parallel
@@ -754,15 +779,10 @@ 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 =
-        poly_space_derived_ptr->get_numbering_inverse();
-      TensorProductPolynomials<dim> *poly_space_derived_ptr_other =
-        dynamic_cast<TensorProductPolynomials<dim> *>(
-          fe_q_other->poly_space.get());
+        this->get_poly_space_numbering_inverse();
       const std::vector<unsigned int> &index_map_inverse_other =
-        poly_space_derived_ptr_other->get_numbering_inverse();
+        fe_q_other->get_poly_space_numbering_inverse();
 
       for (unsigned int i = 0; i < p - 1; ++i)
         for (unsigned int j = 0; j < q - 1; ++j)
@@ -821,15 +841,10 @@ 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 =
-        poly_space_derived_ptr->get_numbering_inverse();
-      TensorProductPolynomials<dim> *poly_space_derived_ptr_other =
-        dynamic_cast<TensorProductPolynomials<dim> *>(
-          fe_q_other->poly_space.get());
+        this->get_poly_space_numbering_inverse();
       const std::vector<unsigned int> &index_map_inverse_other =
-        poly_space_derived_ptr_other->get_numbering_inverse();
+        fe_q_other->get_poly_space_numbering_inverse();
 
       for (unsigned int i1 = 0; i1 < p - 1; ++i1)
         for (unsigned int i2 = 0; i2 < p - 1; ++i2)
@@ -886,10 +901,8 @@ 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 =
-    poly_space_derived_ptr->get_numbering_inverse();
+    this->get_poly_space_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
@@ -1223,10 +1236,8 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_prolongation_matrix(
       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 =
-        poly_space_derived_ptr->get_numbering_inverse();
+        this->get_poly_space_numbering_inverse();
 
       // helper value: step size how to walk through diagonal and how many
       // points we have left apart from the first dimension
@@ -1400,12 +1411,9 @@ 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;
-
-      TensorProductPolynomials<dim> *poly_space_derived_ptr =
-        dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
+      const double                     eps = 1e-15 * q_degree * dim;
       const std::vector<unsigned int> &index_map_inverse =
-        poly_space_derived_ptr->get_numbering_inverse();
+        this->get_poly_space_numbering_inverse();
 
       const unsigned int          dofs1d = q_degree + 1;
       std::vector<Tensor<1, dim>> evaluations1d(dofs1d);