From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 31 Jul 2020 21:55:57 +0000 (+0200)
Subject: Convert all face-related data structures of FiniteElement to vectors
X-Git-Tag: v9.3.0-rc1~1230^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=88cd461ff35f5a81c261c1af903d285db1a3e1bd;p=dealii.git

Convert all face-related data structures of FiniteElement to vectors
---

diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h
index 586a3784ff..8a24356dbc 100644
--- a/include/deal.II/fe/fe.h
+++ b/include/deal.II/fe/fe.h
@@ -2436,7 +2436,7 @@ protected:
    * get_unit_face_support_points() function for a discussion of what
    * contributes a face support point.
    */
-  std::vector<Point<dim - 1>> unit_face_support_points;
+  std::vector<std::vector<Point<dim - 1>>> unit_face_support_points;
 
   /**
    * Support points used for interpolation functions of non-Lagrangian
@@ -2448,7 +2448,7 @@ protected:
    * Face support points used for interpolation functions of non-Lagrangian
    * elements.
    */
-  std::vector<Point<dim - 1>> generalized_face_support_points;
+  std::vector<std::vector<Point<dim - 1>>> generalized_face_support_points;
 
   /**
    * For faces with non-standard face_orientation in 3D, the dofs on faces
@@ -2465,7 +2465,7 @@ protected:
    * no permutation at all. Derived finite element classes have to
    * fill this Table with the correct values.
    */
-  Table<2, int> adjust_quad_dof_index_for_face_orientation_table;
+  std::vector<Table<2, int>> adjust_quad_dof_index_for_face_orientation_table;
 
   /**
    * For lines with non-standard line_orientation in 3D, the dofs on lines
@@ -2496,7 +2496,7 @@ protected:
    * information thus makes only sense if a shape function is non-zero in only
    * one component.
    */
-  std::vector<std::pair<unsigned int, unsigned int>>
+  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
     face_system_to_component_table;
 
   /**
@@ -2521,7 +2521,8 @@ protected:
   /**
    * Likewise for the indices on faces.
    */
-  std::vector<std::pair<std::pair<unsigned int, unsigned int>, unsigned int>>
+  std::vector<
+    std::vector<std::pair<std::pair<unsigned int, unsigned int>, unsigned int>>>
     face_system_to_base_table;
 
   /**
@@ -3122,7 +3123,7 @@ inline std::pair<unsigned int, unsigned int>
 FiniteElement<dim, spacedim>::face_system_to_component_index(
   const unsigned int index) const
 {
-  AssertIndexRange(index, face_system_to_component_table.size());
+  AssertIndexRange(index, face_system_to_component_table[0].size());
 
   // in debug mode, check whether the
   // function is primitive, since
@@ -3141,7 +3142,7 @@ FiniteElement<dim, spacedim>::face_system_to_component_index(
          (typename FiniteElement<dim, spacedim>::ExcShapeFunctionNotPrimitive(
            index)));
 
-  return face_system_to_component_table[index];
+  return face_system_to_component_table[0][index];
 }
 
 
@@ -3162,8 +3163,8 @@ inline std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
 FiniteElement<dim, spacedim>::face_system_to_base_index(
   const unsigned int index) const
 {
-  AssertIndexRange(index, face_system_to_base_table.size());
-  return face_system_to_base_table[index];
+  AssertIndexRange(index, face_system_to_base_table[0].size());
+  return face_system_to_base_table[0][index];
 }
 
 
diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index bb90aaf84c..e46a6ef4c7 100644
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -59,14 +59,9 @@ FiniteElement<dim, spacedim>::FiniteElement(
   const std::vector<bool> &         r_i_a_f,
   const std::vector<ComponentMask> &nonzero_c)
   : FiniteElementData<dim>(fe_data)
-  , adjust_quad_dof_index_for_face_orientation_table(dim == 3 ?
-                                                       this->n_dofs_per_quad() :
-                                                       0,
-                                                     dim == 3 ? 8 : 0)
   , adjust_line_dof_index_for_line_orientation_table(
       dim == 3 ? this->n_dofs_per_line() : 0)
   , system_to_base_table(this->n_dofs_per_cell())
-  , face_system_to_base_table(this->n_dofs_per_face())
   , component_to_base_table(this->components,
                             std::make_pair(std::make_pair(0U, 0U), 0U))
   ,
@@ -109,17 +104,23 @@ FiniteElement<dim, spacedim>::FiniteElement(
   if (this->is_primitive())
     {
       system_to_component_table.resize(this->n_dofs_per_cell());
-      face_system_to_component_table.resize(this->n_dofs_per_face());
       for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
         system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
+
+      face_system_to_component_table.resize(1);
+      face_system_to_component_table[0].resize(this->n_dofs_per_face());
       for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
-        face_system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
+        face_system_to_component_table[0][j] =
+          std::pair<unsigned, unsigned>(0, j);
     }
 
   for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
     system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
+
+  face_system_to_base_table.resize(1);
+  face_system_to_base_table[0].resize(this->n_dofs_per_face());
   for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
-    face_system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
+    face_system_to_base_table[0][j] = std::make_pair(std::make_pair(0U, 0U), j);
 
   // Fill with default value; may be changed by constructor of derived class.
   base_to_block_indices.reinit(1, 1);
@@ -140,7 +141,17 @@ FiniteElement<dim, spacedim>::FiniteElement(
                                   FullMatrix<double>());
     }
 
-  adjust_quad_dof_index_for_face_orientation_table.fill(0);
+  if (dim == 3)
+    {
+      adjust_quad_dof_index_for_face_orientation_table.resize(1);
+
+      adjust_quad_dof_index_for_face_orientation_table[0] =
+        Table<2, int>(this->n_dofs_per_quad(), 8);
+      adjust_quad_dof_index_for_face_orientation_table[0].fill(0);
+    }
+
+  unit_face_support_points.resize(1);
+  generalized_face_support_points.resize(1);
 }
 
 
@@ -653,10 +664,10 @@ FiniteElement<dim, spacedim>::adjust_quad_dof_index_for_face_orientation(
   // the function should also not have been
   // called
   AssertIndexRange(index, this->n_dofs_per_quad());
-  Assert(adjust_quad_dof_index_for_face_orientation_table.n_elements() ==
+  Assert(adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
            8 * this->n_dofs_per_quad(),
          ExcInternalError());
-  return index + adjust_quad_dof_index_for_face_orientation_table(
+  return index + adjust_quad_dof_index_for_face_orientation_table[0](
                    index, 4 * face_orientation + 2 * face_flip + face_rotation);
 }
 
@@ -1057,10 +1068,10 @@ FiniteElement<dim, spacedim>::get_unit_face_support_points() const
   // support points, but only if
   // there are as many as there are
   // degrees of freedom on a face
-  Assert((unit_face_support_points.size() == 0) ||
-           (unit_face_support_points.size() == this->n_dofs_per_face()),
+  Assert((unit_face_support_points[0].size() == 0) ||
+           (unit_face_support_points[0].size() == this->n_dofs_per_face()),
          ExcInternalError());
-  return unit_face_support_points;
+  return unit_face_support_points[0];
 }
 
 
@@ -1069,7 +1080,7 @@ template <int dim, int spacedim>
 bool
 FiniteElement<dim, spacedim>::has_face_support_points() const
 {
-  return (unit_face_support_points.size() != 0);
+  return (unit_face_support_points[0].size() != 0);
 }
 
 
@@ -1080,9 +1091,9 @@ FiniteElement<dim, spacedim>::unit_face_support_point(
   const unsigned int index) const
 {
   AssertIndexRange(index, this->n_dofs_per_face());
-  Assert(unit_face_support_points.size() == this->n_dofs_per_face(),
+  Assert(unit_face_support_points[0].size() == this->n_dofs_per_face(),
          ExcFEHasNoSupportPoints());
-  return unit_face_support_points[index];
+  return unit_face_support_points[0][index];
 }
 
 
diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc
index 83732630e5..c55552234b 100644
--- a/source/fe/fe_abf.cc
+++ b/source/fe/fe_abf.cc
@@ -160,7 +160,7 @@ FE_ABF<dim>::initialize_support_points(const unsigned int deg)
 
   this->generalized_support_points.resize(
     GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
-  this->generalized_face_support_points.resize(n_face_points);
+  this->generalized_face_support_points[0].resize(n_face_points);
 
 
   // These might be required when the faces contribution is computed
@@ -194,7 +194,7 @@ FE_ABF<dim>::initialize_support_points(const unsigned int deg)
 
       for (unsigned int k = 0; k < n_face_points; ++k)
         {
-          this->generalized_face_support_points[k] = face_points.point(k);
+          this->generalized_face_support_points[0][k] = face_points.point(k);
           // Compute its quadrature
           // contribution for each
           // moment.
diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc
index a85ede4162..dd1e4b469e 100644
--- a/source/fe/fe_bdm.cc
+++ b/source/fe/fe_bdm.cc
@@ -333,9 +333,9 @@ FE_BDM<dim>::initialize_support_points(const unsigned int deg)
   QGauss<dim - 1> face_points(deg + 1);
 
   // Copy the quadrature formula to the face points.
-  this->generalized_face_support_points.resize(face_points.size());
+  this->generalized_face_support_points[0].resize(face_points.size());
   for (unsigned int k = 0; k < face_points.size(); ++k)
-    this->generalized_face_support_points[k] = face_points.point(k);
+    this->generalized_face_support_points[0][k] = face_points.point(k);
 
   // In the interior, we only test with polynomials of degree up to
   // deg-2, thus we use deg points. Note that deg>=1 and the lowest
diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc
index 3a37f8fd2b..012f4af037 100644
--- a/source/fe/fe_enriched.cc
+++ b/source/fe/fe_enriched.cc
@@ -440,7 +440,7 @@ FE_Enriched<dim, spacedim>::initialize(
     // If the system is not primitive, these have not been initialized by
     // FiniteElement
     this->system_to_component_table.resize(this->n_dofs_per_cell());
-    this->face_system_to_component_table.resize(this->n_dofs_per_face());
+    this->face_system_to_component_table[0].resize(this->n_dofs_per_face());
 
     FETools::Compositing::build_cell_tables(this->system_to_base_table,
                                             this->system_to_component_table,
@@ -449,8 +449,8 @@ FE_Enriched<dim, spacedim>::initialize(
                                             false);
 
     FETools::Compositing::build_face_tables(
-      this->face_system_to_base_table,
-      this->face_system_to_component_table,
+      this->face_system_to_base_table[0],
+      this->face_system_to_component_table[0],
       *this,
       false);
   }
@@ -468,8 +468,8 @@ FE_Enriched<dim, spacedim>::initialize(
   // this FE sits on the boundary or not. Thus for moment just copy support
   // points from fe system:
   {
-    this->unit_support_points      = fe_system->unit_support_points;
-    this->unit_face_support_points = fe_system->unit_face_support_points;
+    this->unit_support_points         = fe_system->unit_support_points;
+    this->unit_face_support_points[0] = fe_system->unit_face_support_points[0];
   }
 
   // take adjust_quad_dof_index_for_face_orientation_table from FESystem:
diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc
index 8eefdcbbb6..fb4343c0f4 100644
--- a/source/fe/fe_face.cc
+++ b/source/fe/fe_face.cc
@@ -61,12 +61,12 @@ FE_FaceQ<dim, spacedim>::FE_FaceQ(const unsigned int degree)
 {
   // initialize unit face support points
   const unsigned int codim = dim - 1;
-  this->unit_face_support_points.resize(
+  this->unit_face_support_points[0].resize(
     Utilities::fixed_power<codim>(this->degree + 1));
 
   if (this->degree == 0)
     for (unsigned int d = 0; d < codim; ++d)
-      this->unit_face_support_points[0][d] = 0.5;
+      this->unit_face_support_points[0][0][d] = 0.5;
   else
     {
       std::vector<Point<1>> points =
@@ -85,16 +85,16 @@ FE_FaceQ<dim, spacedim>::FE_FaceQ(const unsigned int degree)
               if (codim > 2)
                 p(2) = points[iz][0];
 
-              this->unit_face_support_points[k++] = p;
+              this->unit_face_support_points[0][k++] = p;
             }
-      AssertDimension(k, this->unit_face_support_points.size());
+      AssertDimension(k, this->unit_face_support_points[0].size());
     }
 
   // initialize unit support points (this makes it possible to assign initial
   // values to FE_FaceQ)
   this->unit_support_points.resize(GeometryInfo<dim>::faces_per_cell *
-                                   this->unit_face_support_points.size());
-  const unsigned int n_face_dofs = this->unit_face_support_points.size();
+                                   this->unit_face_support_points[0].size());
+  const unsigned int n_face_dofs = this->unit_face_support_points[0].size();
   for (unsigned int i = 0; i < n_face_dofs; ++i)
     for (unsigned int d = 0; d < dim; ++d)
       {
@@ -106,9 +106,9 @@ FE_FaceQ<dim, spacedim>::FE_FaceQ(const unsigned int degree)
               if (dim == 3 && d == 1)
                 renumber = i / (degree + 1) + (degree + 1) * (i % (degree + 1));
               this->unit_support_points[n_face_dofs * 2 * d + i][e] =
-                this->unit_face_support_points[renumber][c];
+                this->unit_face_support_points[0][renumber][c];
               this->unit_support_points[n_face_dofs * (2 * d + 1) + i][e] =
-                this->unit_face_support_points[renumber][c];
+                this->unit_face_support_points[0][renumber][c];
               this->unit_support_points[n_face_dofs * (2 * d + 1) + i][d] = 1;
               ++c;
             }
@@ -514,7 +514,7 @@ FE_FaceQ<1, spacedim>::FE_FaceQ(const unsigned int degree)
       std::vector<bool>(1, true),
       std::vector<ComponentMask>(1, ComponentMask(1, true)))
 {
-  this->unit_face_support_points.resize(1);
+  this->unit_face_support_points[0].resize(1);
 
   // initialize unit support points (this makes it possible to assign initial
   // values to FE_FaceQ)
diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index 9645e6a0c9..f5da3d679c 100644
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -273,11 +273,11 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order)
     QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
                                           reference_edge_quadrature);
 
-  this->generalized_face_support_points.resize(n_edge_points);
+  this->generalized_face_support_points[0].resize(n_edge_points);
 
   // Create face support points.
   for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
-    this->generalized_face_support_points[q_point] =
+    this->generalized_face_support_points[0][q_point] =
       reference_edge_quadrature.point(q_point);
 
   if (order > 0)
@@ -310,7 +310,7 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order)
             boundary_weights(q_point, i) =
               reference_edge_quadrature.weight(q_point) *
               lobatto_polynomials_grad[i + 1].value(
-                this->generalized_face_support_points[q_point](0));
+                this->generalized_face_support_points[0][q_point](0));
         }
 
       for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
@@ -380,8 +380,8 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
 
       boundary_weights.reinit(n_edge_points + n_face_points,
                               2 * (order + 1) * order);
-      this->generalized_face_support_points.resize(4 * n_edge_points +
-                                                   n_face_points);
+      this->generalized_face_support_points[0].resize(4 * n_edge_points +
+                                                      n_face_points);
       this->generalized_support_points.resize(n_boundary_points +
                                               n_interior_points);
 
@@ -391,8 +391,8 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
           for (unsigned int line = 0;
                line < GeometryInfo<dim - 1>::lines_per_cell;
                ++line)
-            this->generalized_face_support_points[line * n_edge_points +
-                                                  q_point] =
+            this->generalized_face_support_points[0][line * n_edge_points +
+                                                     q_point] =
               edge_quadrature.point(
                 QProjector<dim - 1>::DataSetDescriptor::face(
                   ReferenceCell::get_hypercube(dim - 1),
@@ -421,13 +421,14 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
             boundary_weights(q_point, i) =
               reference_edge_quadrature.weight(q_point) *
               lobatto_polynomials_grad[i + 1].value(
-                this->generalized_face_support_points[q_point](1));
+                this->generalized_face_support_points[0][q_point](1));
         }
 
       // Create support points on faces.
       for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
         {
-          this->generalized_face_support_points[q_point + 4 * n_edge_points] =
+          this
+            ->generalized_face_support_points[0][q_point + 4 * n_edge_points] =
             reference_face_quadrature.point(q_point);
 
           for (unsigned int i = 0; i <= order; ++i)
@@ -436,23 +437,23 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
                 boundary_weights(q_point + n_edge_points, 2 * (i * order + j)) =
                   reference_face_quadrature.weight(q_point) *
                   lobatto_polynomials_grad[i].value(
-                    this->generalized_face_support_points[q_point +
-                                                          4 * n_edge_points](
+                    this->generalized_face_support_points[0][q_point +
+                                                             4 * n_edge_points](
                       0)) *
                   lobatto_polynomials[j + 2].value(
-                    this->generalized_face_support_points[q_point +
-                                                          4 * n_edge_points](
+                    this->generalized_face_support_points[0][q_point +
+                                                             4 * n_edge_points](
                       1));
                 boundary_weights(q_point + n_edge_points,
                                  2 * (i * order + j) + 1) =
                   reference_face_quadrature.weight(q_point) *
                   lobatto_polynomials_grad[i].value(
-                    this->generalized_face_support_points[q_point +
-                                                          4 * n_edge_points](
+                    this->generalized_face_support_points[0][q_point +
+                                                             4 * n_edge_points](
                       1)) *
                   lobatto_polynomials[j + 2].value(
-                    this->generalized_face_support_points[q_point +
-                                                          4 * n_edge_points](
+                    this->generalized_face_support_points[0][q_point +
+                                                             4 * n_edge_points](
                       0));
               }
         }
@@ -485,7 +486,7 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
 
   else
     {
-      this->generalized_face_support_points.resize(4 * n_edge_points);
+      this->generalized_face_support_points[0].resize(4 * n_edge_points);
       this->generalized_support_points.resize(
         GeometryInfo<dim>::lines_per_cell * n_edge_points);
 
@@ -494,8 +495,8 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
           for (unsigned int line = 0;
                line < GeometryInfo<dim - 1>::lines_per_cell;
                ++line)
-            this->generalized_face_support_points[line * n_edge_points +
-                                                  q_point] =
+            this->generalized_face_support_points[0][line * n_edge_points +
+                                                     q_point] =
               edge_quadrature.point(
                 QProjector<dim - 1>::DataSetDescriptor::face(
                   ReferenceCell::get_hypercube(dim - 1),
@@ -3182,7 +3183,7 @@ FE_Nedelec<dim>::convert_generalized_support_point_values_to_dof_values(
                     system_matrix(i, j) +=
                       boundary_weights(q_point, j) *
                       lobatto_polynomials_grad[i + 1].value(
-                        this->generalized_face_support_points[q_point](0));
+                        this->generalized_face_support_points[0][q_point](0));
 
               FullMatrix<double> system_matrix_inv(this->degree - 1,
                                                    this->degree - 1);
@@ -3461,7 +3462,7 @@ FE_Nedelec<dim>::convert_generalized_support_point_values_to_dof_values(
                     system_matrix(i, j) +=
                       boundary_weights(q_point, j) *
                       lobatto_polynomials_grad[i + 1].value(
-                        this->generalized_face_support_points[q_point](1));
+                        this->generalized_face_support_points[0][q_point](1));
 
               FullMatrix<double> system_matrix_inv(this->degree - 1,
                                                    this->degree - 1);
@@ -3540,10 +3541,10 @@ FE_Nedelec<dim>::convert_generalized_support_point_values_to_dof_values(
                                            2 * (k * (this->degree - 1) + l)) *
                           legendre_polynomials[i].value(
                             this->generalized_face_support_points
-                              [q_point + 4 * n_edge_points](0)) *
+                              [0][q_point + 4 * n_edge_points](0)) *
                           lobatto_polynomials[j + 2].value(
                             this->generalized_face_support_points
-                              [q_point + 4 * n_edge_points](1));
+                              [0][q_point + 4 * n_edge_points](1));
 
               system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
               system_matrix_inv.invert(system_matrix);
diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc
index e86950d1c7..8d40cbce7c 100644
--- a/source/fe/fe_p1nc.cc
+++ b/source/fe/fe_p1nc.cc
@@ -28,9 +28,9 @@ FE_P1NC::FE_P1NC()
       std::vector<ComponentMask>(4, ComponentMask(1, true)))
 {
   // face support points: 2 end vertices
-  unit_face_support_points.resize(2);
-  unit_face_support_points[0][0] = 0.0;
-  unit_face_support_points[1][0] = 1.0;
+  unit_face_support_points[0].resize(2);
+  unit_face_support_points[0][0][0] = 0.0;
+  unit_face_support_points[0][1][0] = 1.0;
 
   // initialize constraints matrix
   initialize_constraints();
diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc
index d6664d6ea9..9528d3958e 100644
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -942,7 +942,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::initialize_unit_face_support_points(
   if (dim == 1)
     return;
 
-  this->unit_face_support_points.resize(
+  this->unit_face_support_points[0].resize(
     Utilities::fixed_power<dim - 1>(q_degree + 1));
 
   // find renumbering of faces and assign from values of quadrature
@@ -958,9 +958,9 @@ FE_Q_Base<PolynomialType, dim, spacedim>::initialize_unit_face_support_points(
 
   // The only thing we have to do is reorder the points from tensor
   // product order to the order in which we enumerate DoFs on cells
-  this->unit_face_support_points.resize(support_quadrature.size());
+  this->unit_face_support_points[0].resize(support_quadrature.size());
   for (unsigned int k = 0; k < support_quadrature.size(); ++k)
-    this->unit_face_support_points[face_index_map[k]] =
+    this->unit_face_support_points[0][face_index_map[k]] =
       support_quadrature.point(k);
 }
 
@@ -975,8 +975,8 @@ FE_Q_Base<PolynomialType, dim, spacedim>::
   if (dim < 3)
     return;
 
-  Assert(this->adjust_quad_dof_index_for_face_orientation_table.n_elements() ==
-           8 * this->n_dofs_per_quad(),
+  Assert(this->adjust_quad_dof_index_for_face_orientation_table[0]
+             .n_elements() == 8 * this->n_dofs_per_quad(),
          ExcInternalError());
 
   const unsigned int n = q_degree - 1;
@@ -1007,27 +1007,27 @@ FE_Q_Base<PolynomialType, dim, spacedim>::
       unsigned int i = local % n, j = local / n;
 
       // face_orientation=false, face_flip=false, face_rotation=false
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 0) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 0) =
         j + i * n - local;
       // face_orientation=false, face_flip=false, face_rotation=true
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 1) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 1) =
         i + (n - 1 - j) * n - local;
       // face_orientation=false, face_flip=true,  face_rotation=false
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 2) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 2) =
         (n - 1 - j) + (n - 1 - i) * n - local;
       // face_orientation=false, face_flip=true,  face_rotation=true
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 3) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 3) =
         (n - 1 - i) + j * n - local;
       // face_orientation=true,  face_flip=false, face_rotation=false
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 4) = 0;
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 4) = 0;
       // face_orientation=true,  face_flip=false, face_rotation=true
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 5) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 5) =
         j + (n - 1 - i) * n - local;
       // face_orientation=true,  face_flip=true,  face_rotation=false
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 6) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 6) =
         (n - 1 - i) + (n - 1 - j) * n - local;
       // face_orientation=true,  face_flip=true,  face_rotation=true
-      this->adjust_quad_dof_index_for_face_orientation_table(local, 7) =
+      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 7) =
         (n - 1 - j) + i * n - local;
     }
 
diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc
index 964aae3ba3..099c58900e 100644
--- a/source/fe/fe_q_hierarchical.cc
+++ b/source/fe/fe_q_hierarchical.cc
@@ -1903,7 +1903,7 @@ FE_Q_Hierarchical<dim>::initialize_generalized_face_support_points()
   for (unsigned int i = 1; i < codim; ++i)
     n *= this->degree + 1;
 
-  this->generalized_face_support_points.resize(n);
+  this->generalized_face_support_points[0].resize(n);
 
   Point<codim> p;
 
@@ -1936,7 +1936,7 @@ FE_Q_Hierarchical<dim>::initialize_generalized_face_support_points()
               else
                 p(2) = .5;
             }
-          this->generalized_face_support_points[face_renumber[k++]] = p;
+          this->generalized_face_support_points[0][face_renumber[k++]] = p;
         }
 }
 
diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc
index 310806e5bc..5b4334ae6c 100644
--- a/source/fe/fe_raviart_thomas.cc
+++ b/source/fe/fe_raviart_thomas.cc
@@ -160,7 +160,7 @@ FE_RaviartThomas<dim>::initialize_support_points(const unsigned int deg)
 
   this->generalized_support_points.resize(
     GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
-  this->generalized_face_support_points.resize(n_face_points);
+  this->generalized_face_support_points[0].resize(n_face_points);
 
   // Number of the point being entered
   unsigned int current = 0;
@@ -178,7 +178,7 @@ FE_RaviartThomas<dim>::initialize_support_points(const unsigned int deg)
 
       for (unsigned int k = 0; k < n_face_points; ++k)
         {
-          this->generalized_face_support_points[k] = face_points.point(k);
+          this->generalized_face_support_points[0][k] = face_points.point(k);
           // Compute its quadrature
           // contribution for each
           // moment.
diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc
index 0ddeb83502..f633cc23b2 100644
--- a/source/fe/fe_raviart_thomas_nodal.cc
+++ b/source/fe/fe_raviart_thomas_nodal.cc
@@ -148,7 +148,7 @@ void
 FE_RaviartThomasNodal<dim>::initialize_support_points(const unsigned int deg)
 {
   this->generalized_support_points.resize(this->n_dofs_per_cell());
-  this->generalized_face_support_points.resize(this->n_dofs_per_face());
+  this->generalized_face_support_points[0].resize(this->n_dofs_per_face());
 
   // Number of the point being entered
   unsigned int current = 0;
@@ -164,7 +164,7 @@ FE_RaviartThomasNodal<dim>::initialize_support_points(const unsigned int deg)
       QGauss<dim - 1> face_points(deg + 1);
       Assert(face_points.size() == this->n_dofs_per_face(), ExcInternalError());
       for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
-        this->generalized_face_support_points[k] = face_points.point(k);
+        this->generalized_face_support_points[0][k] = face_points.point(k);
       Quadrature<dim> faces =
         QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
                                               face_points);
@@ -614,7 +614,7 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
   // which returns the support
   // points on the face.
   Quadrature<dim - 1> quad_face_support(
-    source_fe.generalized_face_support_points);
+    source_fe.generalized_face_support_points[0]);
 
   // Rule of thumb for FP accuracy,
   // that can be expected for a
@@ -721,7 +721,7 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
   // which returns the support
   // points on the face.
   Quadrature<dim - 1> quad_face_support(
-    source_fe.generalized_face_support_points);
+    source_fe.generalized_face_support_points[0]);
 
   // Rule of thumb for FP accuracy,
   // that can be expected for a
diff --git a/source/fe/fe_rt_bubbles.cc b/source/fe/fe_rt_bubbles.cc
index b7c66bdf63..9f0b2025db 100644
--- a/source/fe/fe_rt_bubbles.cc
+++ b/source/fe/fe_rt_bubbles.cc
@@ -132,7 +132,7 @@ void
 FE_RT_Bubbles<dim>::initialize_support_points(const unsigned int deg)
 {
   this->generalized_support_points.resize(this->n_dofs_per_cell());
-  this->generalized_face_support_points.resize(this->n_dofs_per_face());
+  this->generalized_face_support_points[0].resize(this->n_dofs_per_face());
 
   // Index of the point being entered
   unsigned int current = 0;
@@ -145,7 +145,7 @@ FE_RT_Bubbles<dim>::initialize_support_points(const unsigned int deg)
       QGaussLobatto<dim - 1> face_points(deg + 1);
       Assert(face_points.size() == this->n_dofs_per_face(), ExcInternalError());
       for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
-        this->generalized_face_support_points[k] = face_points.point(k);
+        this->generalized_face_support_points[0][k] = face_points.point(k);
       Quadrature<dim> faces =
         QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
                                               face_points);
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index ef668c0171..f0392cccaa 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -1422,7 +1422,7 @@ FESystem<dim, spacedim>::build_interface_constraints()
         // data type, first value in pair is (base element,instance of base
         // element), second is index within this instance
         const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
-          n_index = this->face_system_to_base_table[n];
+          n_index = this->face_system_to_base_table[0][n];
 
         // likewise for the m index. this is more complicated due to the
         // strange ordering we have for the dofs on the refined faces.
@@ -1460,7 +1460,8 @@ FESystem<dim, spacedim>::build_interface_constraints()
                     // then use face_system_to_base_table
                     const unsigned int tmp1 =
                       2 * this->n_dofs_per_vertex() + index_in_line;
-                    m_index.first = this->face_system_to_base_table[tmp1].first;
+                    m_index.first =
+                      this->face_system_to_base_table[0][tmp1].first;
 
                     // what we are still missing is the index of m within the
                     // base elements interface_constraints table
@@ -1471,12 +1472,12 @@ FESystem<dim, spacedim>::build_interface_constraints()
                     // dof, we can construct the rest: tmp2 will denote the
                     // index of this shape function among the line shape
                     // functions:
-                    Assert(this->face_system_to_base_table[tmp1].second >=
+                    Assert(this->face_system_to_base_table[0][tmp1].second >=
                              2 * base_element(m_index.first.first)
                                    .n_dofs_per_vertex(),
                            ExcInternalError());
                     const unsigned int tmp2 =
-                      this->face_system_to_base_table[tmp1].second -
+                      this->face_system_to_base_table[0][tmp1].second -
                       2 * base_element(m_index.first.first).n_dofs_per_vertex();
                     Assert(tmp2 < base_element(m_index.first.first)
                                     .n_dofs_per_line(),
@@ -1515,14 +1516,15 @@ FESystem<dim, spacedim>::build_interface_constraints()
 
                     const unsigned int tmp1 =
                       4 * this->n_dofs_per_vertex() + index_in_line;
-                    m_index.first = this->face_system_to_base_table[tmp1].first;
+                    m_index.first =
+                      this->face_system_to_base_table[0][tmp1].first;
 
-                    Assert(this->face_system_to_base_table[tmp1].second >=
+                    Assert(this->face_system_to_base_table[0][tmp1].second >=
                              4 * base_element(m_index.first.first)
                                    .n_dofs_per_vertex(),
                            ExcInternalError());
                     const unsigned int tmp2 =
-                      this->face_system_to_base_table[tmp1].second -
+                      this->face_system_to_base_table[0][tmp1].second -
                       4 * base_element(m_index.first.first).n_dofs_per_vertex();
                     Assert(tmp2 < base_element(m_index.first.first)
                                     .n_dofs_per_line(),
@@ -1553,18 +1555,19 @@ FESystem<dim, spacedim>::build_interface_constraints()
                     const unsigned int tmp1 = 4 * this->n_dofs_per_vertex() +
                                               4 * this->n_dofs_per_line() +
                                               index_in_quad;
-                    Assert(tmp1 < this->face_system_to_base_table.size(),
+                    Assert(tmp1 < this->face_system_to_base_table[0].size(),
                            ExcInternalError());
-                    m_index.first = this->face_system_to_base_table[tmp1].first;
+                    m_index.first =
+                      this->face_system_to_base_table[0][tmp1].first;
 
-                    Assert(this->face_system_to_base_table[tmp1].second >=
+                    Assert(this->face_system_to_base_table[0][tmp1].second >=
                              4 * base_element(m_index.first.first)
                                    .n_dofs_per_vertex() +
                                4 * base_element(m_index.first.first)
                                      .n_dofs_per_line(),
                            ExcInternalError());
                     const unsigned int tmp2 =
-                      this->face_system_to_base_table[tmp1].second -
+                      this->face_system_to_base_table[0][tmp1].second -
                       4 *
                         base_element(m_index.first.first).n_dofs_per_vertex() -
                       4 * base_element(m_index.first.first).n_dofs_per_line();
@@ -1644,7 +1647,9 @@ FESystem<dim, spacedim>::initialize(
     // If the system is not primitive, these have not been initialized by
     // FiniteElement
     this->system_to_component_table.resize(this->n_dofs_per_cell());
-    this->face_system_to_component_table.resize(this->n_dofs_per_face());
+
+    this->face_system_to_component_table.resize(1);
+    this->face_system_to_component_table[0].resize(this->n_dofs_per_face());
 
     FETools::Compositing::build_cell_tables(this->system_to_base_table,
                                             this->system_to_component_table,
@@ -1652,8 +1657,8 @@ FESystem<dim, spacedim>::initialize(
                                             *this);
 
     FETools::Compositing::build_face_tables(
-      this->face_system_to_base_table,
-      this->face_system_to_component_table,
+      this->face_system_to_base_table[0],
+      this->face_system_to_component_table[0],
       *this);
   }
 
@@ -1712,28 +1717,28 @@ FESystem<dim, spacedim>::initialize(
         if (!base_element(base_el).has_support_points() &&
             (base_element(base_el).n_dofs_per_face() > 0))
           {
-            this->unit_face_support_points.resize(0);
+            this->unit_face_support_points[0].resize(0);
             return;
           }
 
 
       // generate unit face support points from unit support points of sub
       // elements
-      this->unit_face_support_points.resize(this->n_dofs_per_face());
+      this->unit_face_support_points[0].resize(this->n_dofs_per_face());
 
       for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
         {
           const unsigned int base_i =
-            this->face_system_to_base_table[i].first.first;
+            this->face_system_to_base_table[0][i].first.first;
           const unsigned int index_in_base =
-            this->face_system_to_base_table[i].second;
+            this->face_system_to_base_table[0][i].second;
 
           Assert(index_in_base <
-                   base_element(base_i).unit_face_support_points.size(),
+                   base_element(base_i).unit_face_support_points[0].size(),
                  ExcInternalError());
 
-          this->unit_face_support_points[i] =
-            base_element(base_i).unit_face_support_points[index_in_base];
+          this->unit_face_support_points[0][i] =
+            base_element(base_i).unit_face_support_points[0][index_in_base];
         }
     });
 
@@ -1810,7 +1815,7 @@ FESystem<dim, spacedim>::initialize(
     init_tasks += Threads::new_task([&]() {
       // the array into which we want to write should have the correct size
       // already.
-      Assert(this->adjust_quad_dof_index_for_face_orientation_table
+      Assert(this->adjust_quad_dof_index_for_face_orientation_table[0]
                  .n_elements() == 8 * this->n_dofs_per_quad(),
              ExcInternalError());
 
@@ -1821,12 +1826,12 @@ FESystem<dim, spacedim>::initialize(
         {
           const Table<2, int> &temp =
             this->base_element(b)
-              .adjust_quad_dof_index_for_face_orientation_table;
+              .adjust_quad_dof_index_for_face_orientation_table[0];
           for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
             {
               for (unsigned int i = 0; i < temp.size(0); ++i)
                 for (unsigned int j = 0; j < 8; ++j)
-                  this->adjust_quad_dof_index_for_face_orientation_table(
+                  this->adjust_quad_dof_index_for_face_orientation_table[0](
                     index + i, j) = temp(i, j);
               index += temp.size(0);
             }
@@ -2343,13 +2348,14 @@ Point<dim - 1>
 FESystem<dim, spacedim>::unit_face_support_point(const unsigned int index) const
 {
   AssertIndexRange(index, this->n_dofs_per_face());
-  Assert((this->unit_face_support_points.size() == this->n_dofs_per_face()) ||
-           (this->unit_face_support_points.size() == 0),
+  Assert((this->unit_face_support_points[0].size() ==
+          this->n_dofs_per_face()) ||
+           (this->unit_face_support_points[0].size() == 0),
          (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints()));
 
   // let's see whether we have the information pre-computed
-  if (this->unit_face_support_points.size() != 0)
-    return this->unit_face_support_points[index];
+  if (this->unit_face_support_points[0].size() != 0)
+    return this->unit_face_support_points[0][index];
   else
     // no. ask the base element whether it would like to provide this
     // information
diff --git a/source/fe/fe_trace.cc b/source/fe/fe_trace.cc
index 577c3ae778..0f990292a5 100644
--- a/source/fe/fe_trace.cc
+++ b/source/fe/fe_trace.cc
@@ -53,7 +53,7 @@ FE_TraceQ<dim, spacedim>::FE_TraceQ(const unsigned int degree)
     FETools::hierarchic_to_lexicographic_numbering<dim - 1>(degree));
 
   // Initialize face support points
-  this->unit_face_support_points = fe_q.get_unit_face_support_points();
+  this->unit_face_support_points[0] = fe_q.get_unit_face_support_points();
 
   // initialize unit support points (this makes it possible to assign initial
   // values to FE_TraceQ). Note that we simply take the points of fe_q but
diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc
index 75af427581..a3545d4e8d 100644
--- a/source/simplex/fe_lib.cc
+++ b/source/simplex/fe_lib.cc
@@ -120,8 +120,8 @@ namespace Simplex
             this->unit_support_points.emplace_back(0.0, 0.0);
 
             // TODO
-            this->unit_face_support_points.emplace_back(0.0);
-            this->unit_face_support_points.emplace_back(1.0);
+            this->unit_face_support_points[0].emplace_back(0.0);
+            this->unit_face_support_points[0].emplace_back(1.0);
           }
         else if (degree == 2)
           {
@@ -133,9 +133,9 @@ namespace Simplex
             this->unit_support_points.emplace_back(0.5, 0.0);
 
             // TODO
-            this->unit_face_support_points.emplace_back(0.0);
-            this->unit_face_support_points.emplace_back(1.0);
-            this->unit_face_support_points.emplace_back(0.5);
+            this->unit_face_support_points[0].emplace_back(0.0);
+            this->unit_face_support_points[0].emplace_back(1.0);
+            this->unit_face_support_points[0].emplace_back(0.5);
           }
         else
           {
@@ -152,9 +152,9 @@ namespace Simplex
             this->unit_support_points.emplace_back(0.0, 0.0, 1.0);
 
             // TODO
-            this->unit_face_support_points.emplace_back(1.0, 0.0);
-            this->unit_face_support_points.emplace_back(0.0, 1.0);
-            this->unit_face_support_points.emplace_back(0.0, 0.0);
+            this->unit_face_support_points[0].emplace_back(1.0, 0.0);
+            this->unit_face_support_points[0].emplace_back(0.0, 1.0);
+            this->unit_face_support_points[0].emplace_back(0.0, 0.0);
           }
         else if (degree == 2)
           {
@@ -170,12 +170,12 @@ namespace Simplex
             this->unit_support_points.emplace_back(0.0, 0.5, 0.5);
 
             // TODO
-            this->unit_face_support_points.emplace_back(1.0, 0.0);
-            this->unit_face_support_points.emplace_back(0.0, 1.0);
-            this->unit_face_support_points.emplace_back(0.0, 0.0);
-            this->unit_face_support_points.emplace_back(0.5, 0.5);
-            this->unit_face_support_points.emplace_back(0.0, 0.5);
-            this->unit_face_support_points.emplace_back(0.5, 0.0);
+            this->unit_face_support_points[0].emplace_back(1.0, 0.0);
+            this->unit_face_support_points[0].emplace_back(0.0, 1.0);
+            this->unit_face_support_points[0].emplace_back(0.0, 0.0);
+            this->unit_face_support_points[0].emplace_back(0.5, 0.5);
+            this->unit_face_support_points[0].emplace_back(0.0, 0.5);
+            this->unit_face_support_points[0].emplace_back(0.5, 0.0);
           }
         else
           {